Field sparsening for the construction of the correlation functions in lattice QCD

Two field-sparsening methods, namely the sparse-grid method and the random field selection method, are used in this paper for the construction of the 2-point and 3-point correlation functions in lattice QCD. We argue that, due to the high correlation among the lattice correlators at different field points associated with source, current, and sink locations, one can save a lot of computational time by performing the summation over a subset of the lattice sites. Furthermore, with this strategy, one only needs to store a small fraction of the full quark propagators. It is found that the number of field points can be reduced by a factor of $\sim$100 for the point-source operator and a factor of $\sim$1000 for the Gaussian-smeared operator, while the uncertainties of the correlators only increase by $\sim$15\%. Therefore, with a modest cost of the computational resources, one can approach the precision of the all-to-all correlators using the field-sparsening methods.


I. INTRODUCTION
Lattice QCD provides a non-perturbative approach to the numerical solution of Quantum Chromodynamics (QCD), which is believed to be the basic theory of strong interactions among quarks and gluons. With the development of the cutting-edge supercomputers, the new algorithms and the advanced methodologies, lattice QCD is now playing an increasingly important role in the understanding of low-energy QCD.
A typical way to gain a better efficiency in a numerical lattice QCD calculation is to reduce the redundant costs. Here are some examples.
• To guarantee an accurate generation of a quark propagator, the residual criterion in the conjugate gradient (CG) inversion for the quark propagator is usually set to 10 −8 or even smaller. When realizing that most of the CG iterations can be saved, the all mode averaging (AMA) technique [1] is proposed, where the residual criterion is raised to a level of ∼ 10 −4 . In this way the computational time is significantly reduced while the physical quantities can still be obtained with no bias by performing a correction, which compensates the systematic effects from the approximated propagators by adding the difference between some samples of the precise correlators and the approximated ones.
• In the calculation of hadronic light-by-light contributions to the muon anomalous magnetic moment, a four-point hadronic function with vector-currents is required. To construct such a four-current correlator, one expects a spacetime summation over the locations of at least three of the four currents which is very challenging numerically.
Realizing that in the connected diagram, when the locations of two vector currents are separated with a spacetime distance r, the hadronic function falls exponentially with the increase of r, an importance sampling is introduced to evaluate the stochastic sum over r efficiently [2]. Therefore, in the important regions where r 1 fm, the summation is run in a complete way while in the other regions r > 1 fm, the contributions are calculated with a probability of p(r) ∝ 1/r 3.5 . In this way, much less computational resources are spent in the very long distance region, where the lattice correlation functions mainly contribute noise rather than the signal.
• In many cases, it is appealing to utilize the translational invariance and construct the correlation function using the all-to-all propagators [3]. As a result, the information over the whole spacetime volume is summed and one expects to gain a good precision for the correlation function. On the other hand, generating all-to-all propagators is quite time consuming. Since the correlation functions from the same configuration is highly correlated, one can achieve nearly the same precision by averaging over part of the source and sink locations. Such kind of techniques, called field-sparsening methods here, are the main focus of this paper.
The field-sparsening techniques have been studied by Detmold and Murphy in Ref. [4], where a sparse-grid technique is introduced. An earlier application of the sparse grid can be traced back to a study from χQCD collaboration in 2010 [5]. In this work, in addition to the sparse-grid approach, we developed another field-sparsening approach, which we will refer to as the random field selection method. We will study and compare these two methods in some detail. The paper is organized as follows. We start with Sect. II by introducing the two field-sparsening techniques mentioned above. In Sect. III we discuss the results of the 2-point functions for pion and proton with both point-source and Gaussian-smeared-source operators. The advantages and disadvantages for each field-sparsening method are also discussed. In Sect. IV we employ a simple model to analyze the sources of the uncertainties for the random field selection method. In Sect. V we extend the study to the 3-point function, where the proton axial charge g A is used as an example to demonstrate the efficiency of the field sparsening methods.

II. FIELD SPARSENING METHODS
In the lattice QCD calculation of a generic n-point function, with O i (x i ) being the interpolating operators at temporal-spatial point x i = (x 0 i , x i ) and Λ full the full set of spatial lattice points, one needs to perform the volume summation (n − 1) times. This results in a computational cost of (L 3 ) n−1 in the quark contraction, with L being the spatial lattice size. If one wants to gain a better precision by making another spatial-volume average over the locations of x n , then the cost becomes (L 3 ) n . The typical size of L 3 is about 10 4 -10 6 for practical lattice QCD simulations. Given a relatively large lattice, the complexity of O(L 6 ) usually exceeds the capability of the current lattice QCD calculations. In many cases, the techniques of all-to-all propagators [3] or sequential-source propagators [6] are used to reduce the computational costs. On the other hand, there are some limitation for the usage of these propagators. For example, the all-to-all propagators work more efficiently in the mesonic sector than in the baryonic sector. For the sequential- For simplicity, we start with the 2-point correlation function as an example to introduce the field-sparsening techniques. The standard 2-point correlation function with zero spatial momentum insertion is written as where the subscript "full" indicates that the summation of x runs over all the sink location points. By using field sparsening, one can replace the summation x∈Λ full by L 3 N Λ x∈Λ , where Λ is a subset of Λ full , which contains only N Λ location points. Due to the high correlation in the lattice data, we expect that the replacement does not increase the noise much but reduces the propagator storage and contraction time for modest size N Λ . In Eq. (2) only one source location ( x 0 , t 0 ) is used. In principle one can use multiple source locations and write the correlation functions as where the source spatial location takes the value from the set Λ 0 and source time slice takes the value from Λ t . The size of the set Λ 0 and Λ t is given by N Λ 0 and N Λt , respectively.
In this work we will compare two different sets for Λ (Λ 0 ), namely the sparse-grid method and the random field selection method, and determine the optimal values for N Λ and N Λ 0 in each case.

A. Sparse-grid method
Following the sparse-grid method introduced in Ref. [4], the set Λ is chosen as, where k is an integer factor of L. By using this setup, a L 3 -point spatial lattice is reduced to a (L/k) 3 -point one. For all the time slices t 0 and t 0 + t used in Eq. (2), the same sparse grid set is implemented. One can also extend the type of the sparse grid to type II : {(n 1 , n 2 , n 3 ) 0 ≤ n i < L; n i = 0(mod k); n 1 + n 2 + n 3 = 0(mod 2k)}, type III : {(n 1 , n 2 , n 3 ) 0 ≤ n i < L; n i = 0(mod k); n i + n j = 0(mod 2k), i = j}.
With the above definitions, we have N Λ = L 3 k 3 , L 3 2k 3 , L 3 4k 3 for type I, II and III, respectively. In our numerical study, we use a lattice gauge ensemble with L = 24 and pick up 15 values for N Λ . For convenience, these 15 cases are labelled by an integer denoted as N th , running from 0 to 14, and the corresponding values for N Λ are given by the following list, It means that we have N Λ = 24 3 for N th = 0 and N Λ = 1 for N th = 14.
Note that the sets of type I, II, III always include the location point of (n 1 , n 2 , n 3 ) = (0, 0, 0). We therefore can consider it as a reference point. To reduce the correlation in the lattice calculation, for each configuration, one can shift the reference point randomly with Another choice of Λ is called random field selection, with Λ chosen as The random numbers for n i vary when the time slice t and configuration trajectory change.
In principle, we can use any value for N Λ . To favor a comparison with the sparse-grid method, we use the same choices of N Λ as that in Eq. (6).

III. 2-POINT FUNCTION
In this section we will present the results for 2-point correlation functions. The calculation is performed using a gauge ensemble with 2 + 1 + 1 clover-improved Wilson twisted mass quarks, generated by the ETM Collaboration [7]. The lattice volume is 24 3 × 48 with a pion mass m π ≈ 350 MeV and a lattice spacing a ≈ 0.093 fm. In total 91 configurations are utilized in this analysis.
We use Gaussian-smeared-source propagators to construct the smeared-source smeared- with α = π for pion and α = p for proton. We use 24 time slices for t 0 . For each time slice, we use 8 random source locations for x 0 . Thus we have N Λt = 24 and N Λ 0 = 8. For the sink location x, it can be summed over the full set Λ full or field-sparsening set Λ. and C α (t + 1) as inputs. For the two field-sparsening methods, the data points are slightly shifted horizontally to favor a comparison with that from the full set.
From Fig. 1 a clear enhancement of the excited-state contamination is found in the sparsegrid method at N Λ = 2 3 . This is due to the mixing of the hadron states with high momenta.
Let us take the sparse-grid set type I as an example. The summation over Λ can be written  with Therefore the higher-momentum modes with m = 0 mix with the zero-momentum mode. As a consequence, the excited-state contaminations increase as N Λ decreases. The situation becomes even more problematic when one targets on the calculation of the correlation functions with large momentum transfer. In this case it is possible that a state with an assigned momentum mixes with the states carrying smaller momenta and consequently the correlation functions are distorted by the low-lying states.
For the random-field selection method, there is no enhancement of the excited-state contribution. On the other hand, the statistical errors become larger due to the random noise arising from the field sparsening. Note that the correlation among C α (t) at various t is weakened by the random-field selection. As a result, the uncertainty for the effective mass can be further reduced if one performs a fit over a temporal window.
In Fig. 2 The efficiency of the field sparsening likely depends on the type of interpolating operators used in the lattice calculation. As the Gaussian-smeared-source propagator is more correlated than the point-source propagator, we expect that the field sparsening works more efficiently for the former case. To confirm this conjecture we calculate the effective masses for the pion and proton using both Gaussian-smeared-source and point-source propagators. To save the cost, we keep N Λt at 24 while reducing N Λ 0 to 1 for both smeared-source smearedsink (s-s) and point-source point-sink (p-p) correlation functions. The effective masses from the correlated fit are shown in the left panels of Fig. 3, while the ratios between the statistical uncertainty and the effective masses m full α for α = π and p are shown in the corresponding right panels. With the same statistics, the results of p-p correlators are noisier than that  of the s-s ones. The field-sparsening method works less efficiently in the p-p corelators as we expect. Nevertheless, we find that at N th = 5 (N Λ = 108) the uncertainty of the p-p effective mass increases only 15% for the pion and 18% for the proton. This implies that one can reduce the field points by a factor of 128 in trade with a small increase of the statistical error. Two order of magnitude reduction in the propagator storage and I/O data transfer would make a typical lattice calculation much easier.

IV. NOISE FROM RANDOM FIELD SELECTION METHOD
When using the random field selection method, the correlation function receives two types of the noise. The first is the gauge noise, δ gauge , and the second is the noise from the selection of the random field points, δ rand . The increase of the gauge noise is an inevitable price when the number of the field points is reduced in the summation. Since the lattice data at the different points are highly correlated, we expect the gauge noise only increases mildly. We therefore focus on the estimation of δ rand .
We write the time dependence of the 2-point correlation function as where the sign of . = reminds us that the excited-state contributions are neglected at sufficiently large t. The energy E satisfies the dispersion relation E = m 2 + p 2 with m the hadron's mass.
To fully understand the impact on the uncertainty of the correlation function from the random field selection, one needs to determine the weight functionÃ( p). Considering the fact that we have used the Gaussian-smeared source and sink in our calculation, we assume that the weight function in the coordinate space is given by a Gaussian distribution Then the weight function in the momentum space is given by the Fourier transformation of  formula to fit the pion lattice correlator with three free parameters, namely A 0 , σ and m.
We obtain m = 0.165(1) which is consistent with the effective mass m full π calculated before. Furthermore, we can obtain the result for C(r, t) by averaging the lattice data of C( x, t) in a range of r ≤ | x| < r + 1. In Fig. 4 we show the C(r, t) at t = 10 together with the fitting curve. The good agreement suggests that the functional form given in Eq. (10) describes the lattice data well at large distances.
As a next step, we use the parameters A 0 , σ and m and Eqs. (10) and (12) to construct the correlation functionĈ( x, t). As the gauge noise is eliminated inĈ( x, t), we can isolate the noise δ rand by replacing the lattice correlator C( x, t) withĈ( x, t) in the random field selection method. In Fig. 5 we compare the uncertainties for the three different types of 2-point functions at t = 10, 15 and 20: The blue bar indicates the uncertainty of the correlation function with random field selection method with N Λ = 1. The orange bar shows the uncertainty for the data with N Λ = Λ full while the green bar indicates the size of δ rand only. It is noticed that the size of δ rand is much smaller than the gauge noise and thus can be safely neglected. We have thus reached a conclusion that the precision of the pion 2-point functions using the random field selection method is equally good as that using the sparse-grid method. The random field selection is theoretically cleaner as there is no enhancement of the excited-state contamination.

V. 3-POINT FUNCTION
In this section, we extend the random field selection method to the 3-point function, where we calculate the proton axial charge g A as an example which is one of the most fundamental quantities in nuclear physics. A global average of the lattice results of g A can be found in the latest FLAG review [8]. In this study we mainly focus on the efficiency of the random field selection method rather than the full control of various systematic effects for g A . We start with the 3-point and 2-point functions as and where O p is the Gaussian-smeared operator for the proton, A 3 is the axial vector current with the polarization in the z direction and P is the spin projection operator.
In order to compare the field-sparsening results with the full-size ones, we use the sequential-source propagators, which start from the source (t 0 , x 0 ), go through the current insertion point (t 0 + t i , x) and end at the sink location (t 0 + t s , x s ). We use 24 time slices for t 0 . For each t 0 , one random source is used and the time t i and t s are fixed as t i = 5 and This is very consistent with the observation in the effective masses of the 2-point functions.
Thus we can conclude that the field-sparsening method seems to work equally well for both 2-point and 3-point functions.
As a next step, we want to determine the optimal values of N . The data points with various N Λ 0 -N Λs pairs are plotted using the different symbols. In the right panel of Fig. 7, the statistical uncertainty of C 3pt /C 2pt as a function of N Λ 0 = N Λs are shown. We find that with N Λ 0 changing from 1 to 4, the statistical error drops relatively fast. From N Λ 0 = 4 to 8, the error decreases slower. Due to the high correlation, the change of the uncertainty is very mild from N Λ 0 = 8 to 32. It is unnecessary to move on to larger N Λ 0 as one can expect that the precision at N Λ 0 = N Λs = 8 is close to the best precision using the all-to-all setup (N Λ 0 = N Λs = L 3 ). In practise, we can choose N Λ 0 = 4 or 8 and invest the additional computational resources in accumulating data from more gauge configurations. In Fig. 7, the ratio of C 3pt /C 2pt at t s = 10 is shown. As a next step we add the lattice results at another two values of t s (t s = 8 and 12) and vary t i − t s /2 in a range of The corresponding results are shown in Fig. 8. At large time separation, the time dependence of C 3pt /C 2pt can be approximated by a two-state form as where ∆ is the energy difference between the excited state and the ground state. The  Figure 9. The proton axial charge g A as a function of N th . Here N th is related to the number of field points at the location of current insertion. We use N Λs = N Λ 0 = 32 and N Λt = 24.
As a last step, we calculate g A with different value of N Λ , which is the number of the field points at the location of current insertion (t 0 + t i , x i ). In Fig. 9 we show g A as a function of N th . By increasing N th from 0 to 5 (or equavilently reducing N Λ from L 3 to 108), we find that the uncertainty of g A only increases by 10%. Therefore, by using the field sparsening, one can reduce the size of the Gaussian-smeared source point-sink propagators by a factor of 128 at the expense of 10% increase in the statistical error of the correlation function.

VI. CONCLUSION
In this work we make an exploratory study on the field-sparsening methods. The observables under investigation include the pion and proton 2-point correlation functions and the proton axial charge g A involving the 3-point functions. For the sparse-grid method, the results are not affected by the noise from field sparsening but receive additional excited-state contamination from the higher-momenta states. For the random field selection method, the situation is just the opposite. There is no enhancement of the excited-state contamination, but the correlation functions are affected by the noise from random selection. Fortunately, we confirm, in both numerical lattice results and a model analysis mimicking pion correlator, that the noise from the random selection can be safely neglected.
In the calculation, we construct the correlation function using both Gaussian-smeared operator and the point-source operator. We find that Gaussian-smeared correlators do have higher correlation than the point-source ones. At the expense of a ∼15% increase in the statistical error, we can reduce the number of field points by a factor of ∼100 for the point-source operator and a factor of ∼1000 for the Gaussian-smeared operator. This is a surprisingly significant reduction in the size of the propagator. It also makes the storage of propagators much easier and saves the time for both I/O and quark contraction. Another interesting observation is that by using the field sparsening methods, one can approach the precision of the all-to-all correlators with the modest cost of the computational resources.
Due to its high efficiency, we can foresee a vast application prospect of the field-sparsening methods proposed in this paper.