A new class of variance reduction techniques using lattice symmetries
Abstract
We present a general class of unbiased improved estimators for physical observables in lattice gauge theory computations which significantly reduces statistical errors at modest computational cost. The error reduction techniques, referred to as covariant approximation averaging, utilize approximations which are covariant under lattice symmetry transformations. We observed cost reductions from the new method compared to the traditional one, for fixed statistical error, of 16 times for the nucleon mass at MeV (DomainWall quark) and 2.620 times for the hadronic vacuum polarization at MeV (Asqtad quark). These cost reductions should improve with decreasing quark mass and increasing lattice sizes.
pacs:
11.15.Ha,12.38.Gc,07.05.TpAs nonperturbative computations using lattice gauge theory are applied to a wider range of physically interesting observables, it is increasingly important to find numerical strategies that provide precise results. In Monte Carlo simulations our reach to important physics is still often limited by statistical uncertainties. Examples include hadronic contributions to the muon’s anomalous magnetic moment Aubin et al. (2012), nucleon form factors and structure functions Yamazaki et al. (2009), including nucleon electric dipole moments Shintani et al. (2005); Berruto et al. (2006); Shintani et al. (2007, 2008), hadron matrix elements relevant to flavor physics (, amplitudes) Blum et al. (2012), and multihadron state physics Savage (2012), to name only a few.
As a generalization of lowmode averaging (LMA) Giusti et al. (2004); DeGrand and Schaefer (2004), we present a class of unbiased statistical error reduction techniques, utilizing approximations that are covariant under lattice symmetry transformations. LMA has worked well in cases where low eigenmodes of the Dirac operator dominate Neff et al. (2001): low energy constants in the regime Giusti et al. (2003, 2004); DeGrand and Schaefer (2005); Luscher (2007); Fukaya et al. (2007), pseudoscalar meson masses and decay constants Giusti and Necco (2006); Noaki et al. (2008); Li et al. (2010), an so on. With a modest increase in computational cost, the generalized method can reduce statistical errors by an order of magnitude, or more, even in cases where LMA fails.
Unlike LMA, we account for all modes of the Dirac operator, averaging over (most of) the lattice volume, with modest additional computational cost. The alltoall methods Bali et al. (2005); Foley et al. (2005) implement this stochastically for the higher modes, while treating the lowmodes exactly. For expectation values invariant under translations, statistics effectively increase by averaging over the whole lattice. The alltoall method is advantageous when the stochastic noise introduced in the target observable is comparable to, or smaller than, the gauge field fluctuations of the ensemble Bali et al. (2010a), which typically holds only for many random source vectors per measurement. The error reduction techniques presented here, which do not rely on stochastic noise, are potentially more effective, provided an inexpensive approximation can be found for the desired observable.
In lattice gauge theory simulations an ensemble of gauge field configurations is generated randomly, according to the Boltzmann weight, , where is the latticeregularized action. The expectation value of a primary, covariant observable, ,
(1) 
is estimated as the ensemble average, over a large number of configurations, . Here, we primarily consider observables made of fermion propagators computed on the background gauge configuration .
By exploiting lattice symmetry transformations , that transform , a general class of variance reduction techniques is introduced. First construct an approximation to which must fulfill the following conditions,
 appx1

should fluctuate closely with ,
, and , where .  appx2

the cost to compute is smaller than ’s, .
 appx3

is covariant under a lattice symmetry transformation, , (in the examples below, a stronger condition holds: is covariant on each configuration, rather than on average, ).
Note and refers to the approximations before and after applying a symmetry transformation .
Using and one can define an improved observable
(2)  
where an average over symmetry transformations in is taken.
For appx1, the statistical error of is
(3) 
which can be made smaller than the original (err) by a judicious choice of . The fluctuation from , the first term in (3), is suppressed due to , while the second term is reduced by without too much additional cost as required by appx2 (correlations among , and have been ignored, which is a good approximation for noisy observables or large volume). Due to covariance, appx3, it is easy to prove the ensemble averages of (primary observables) , and are all equal, so the improved estimator (2) is unbiased,
The idea of exploiting covariance Giusti et al. (2004); DeGrand and Schaefer (2004) to improve statistical errors has a wider range of applicability than LMA, so in general we call it covariant approximation averaging (CAA). Several comments on CAA follow. From Eq. 3 the accuracy of the approximation (appx1) should be precise enough so that the statistical error from is below, say, onehalf of the desired final precision. Too accurate an approximation wastes resources. In , most of the statistical fluctuation is carried by , which is reduced by averaging over measurements with smaller cost (appx2). Balance between these opposing parts of the method allows CAA to reduce statistical errors significantly while keeping the computational cost low.
In the framework of CAA the best choice of approximation depends on the target observables and lattice parameters such as quark mass and volume. In principle, any set of lattice symmetries, , can be used in CAA. We limit ourselves to the case of translation symmetries in the following examples.
The first example is LMA. In LMA eigensystems of the Hermitian Dirac operator are obtained for the part of the spectrum closest to zero,
(4)  
(5) 
which is then used to construct, through spectral decomposition, the lowmode approximation of the fermion propagator,
(6)  
(7) 
is the total dimension of the Dirac matrix. The recipe for LMA in terms of the CAA master Eq. (2) is shown in left column of Table 1.
LMA algorithm  AMA algorithm  
1:  Compute lowmodes of  1:  if 
Compute lowmode of  
2:  Set source and invariant inital guess  
3:  Compute exact and precisely (use deflation if exits)  
4:  Repeat for in (6)  4:  Repeat for in (8) 
and  and using  
deflated CG (if )  
5:  5:  ;  
6:  Set shifted source and invariant inital guess  
7:  Average  7:  Average 
over to get  over to get  
8: 
Although LMA is particularly good for observables dominated by lowmodes, such as the single pion state for lighter fermion masses, LMA does not work so well for heavier hadrons or when the quark mass is heavier Giusti and Necco (2006); Li et al. (2010) (see also Bali et al. (2010b) for dependence on parity of states and (non)Hermiticity of Dirac operators). This is due to the truncation of the sum in (6), i.e., for .
One could improve the above by constructing a polynomial for and using it to obtain a better (allmode) approximation of the propagator above :
(8)  
(9) 
where is a polynomial of degree , From (8) and (9), one computes the approximate propagator using in the subspace orthogonal to the eigenvectors below ,
(10) 
with number of lowmodes . In analogy to LMA, we refer to the above as allmode averaging (AMA). A recipe similar to LMA is shown in the right column of Table 1.
As emphasized in Giusti et al. (2004) approximate eigenvectors can be used in LMA (and AMA) to reduce the cost of this part of the calculation. We have not done that as we find the cost of computing them exactly is not too burdensome and is partly recouped in the deflation of the Dirac operator.
Among many different ways Luscher (1994); Borici and de Forcrand (1995); Kamleh and Peardon (2012) to obtain , one of the easiest is to use the polynomial implicitly generated by an iterative linear solver such as conjugate gradient (CG). For example (8) can be implemented as a CG solution using the lowmode approximation applied to the source vector (the coefficients of depend on ) as the starting vector, , which is nothing but a deflated CG with iteration number set to the degree of the polynomial, . One can either fix (number of iterations) or the CG residual vector stopping criterion. Either satisfies the covariance condition (appx3). This particular construction of is called the truncated solver method (TSM) Bali et al. (2010a). The difference with AMA is that TSM is applied in Bali et al. (2010a) to a random source, and the unbiased result is guaranteed by stochasticity while AMA relies on covariance, so it does not need the random source.
In Li et al. (2010) lowmodes are utilized with noise to compute manytoall hadron correlation functions for variance reduction. One may also choose in (9), not to use eigenvectors at all. This may be effective for heavier quark masses, but for lighter quarks one needs a larger degree polynomial for an accurate approximation and is likely more costeffective.
To ensure unbiasness one should check, on a few configurations, the covariance of the particular implementation of the approximation by computing the approximation explicitly on a transformed gauge field to compare with the original gauge field to see that they are equivalent to numerical precision.
To compare the LMA and AMA methods, we use the 2+1 flavor DomainWall fermion (DWF) ensemble generated by the RBC/UKQCD collaboration Aoki et al. (2011) with lattice size , extra dimension size , and Iwasaki gauge action (, or GeV). The lowmodes of the Hermitian DWF Dirac operator are obtained using a 4Devenoddpreconditioned, shifted Lanczos algorithm Neff et al. (2001) with accuracy . The eigenmodes are used for LMA as in Eqs. (6) and (7), to deflate the CG, and to evaluate the lowmode parts of both and , and similarly for AMA as in Eqs. (8) and (9). In this paper we compute 180 lowmodes for light quark mass and 400 lowmodes for .
We adopt translational symmetry on the lattice as and take propagator source locations, starting from the origin, separated by 12 lattice units in space and 16 in time, and the total set of translations numbers . For AMA, the stopping condition of the “sloppy CG” for our approximation is while it is in Yamazaki et al. (2009). Note that when using an evenodd preconditioned Dirac operator, LMA and AMA guarantee unbiased estimators for translations by an even number of sites (appx3). We have explicitly checked this in our calculations.
Table 2 lists the relative statistical errors for various hadronic twopoint correlation functions computed using LMA, AMA, and the original CG method, for . All were obtained with the same Gaussian smeared sources and point (Gaussian) sinks for pseudoscalar and vector (Nucleon) used in Yamazaki et al. (2009). At short distance , there is no improvement between the original and LMA cases, except in the pseudoscalar (PS) channel. This is because the contribution of higher modes is still important in the shortdistance region. Although for LMA could be taken as large as the lattice size with modest cost, we set since larger is not effective due to correlations between nearby gauge fields in our examples. On the other hand, AMA dramatically reduces the errors (more than 46) for all channels (and different momenta) and for all distances. In this example the variance reduction by AMA comes almost entirely from the second term in Eq. (3) since is very close to one ( for ), even though the residual stopping criterion used for is loose (). For LMA at short distance so the error from is significant. We also confirm that for the PS channel both LMA and AMA yield improvement, with even in the short distance region, as suggested previously for LMA using overlap fermions Noaki et al. (2008); Li et al. (2010). For is somewhat smaller (), so the contribution from is more significant. Only 180 lowmodes were used for .
Hadron  Original [%]  LMA [%]  AMA [%]  

N  4  6.9  5.0  1.5 
8  9.2  3.2  1.9  
12  23  4.8  3.5  
PS  4  4.5  0.98  0.86 
12  4.9  0.91  0.86  
28  5.0  1.3  1.3  
V  4  3.9  2.9  0.6 
8  5.2  2.1  1.1  
12  12  3.4  2.3 
Figure 1 shows the nucleon effective mass using LMA and AMA for the data in Tab. 2, and Tab. 3 compares these to an earlier high statistics study of nucleon structure functions Yamazaki et al. (2009). The rightmost panel in Fig. 1 shows significant improvement of the effective mass plateau for AMA. Using the same fitting range, the precision of the nucleon mass attained with AMA is smaller by more than a factor of 1.5 compared to the high statistics study Yamazaki et al. (2009) where 3728 and 1424 measurements were made for and 0.01, respectively. The improved statistics make it easier to choose the fit range based on , as seen in Fig. 1. LMA for nucleon masses was examined in Giusti and Necco (2006).
,  

sink  fit range  LMA  AMA  High stat.  
0.005  pt  812  1.1391(145)  1.1413(61)  1.1561(104) 
0.005  gauss  612  1.1305(143)  1.1420(58)  1.1481(100) 
0.01  pt  915  1.2446(164)  1.2363(59)  1.2101(89) 
0.01  gauss  715  1.2240(148)  1.2268(60)  1.2169(93) 
Most of the cost of AMA comes from the lowmode and sloppy CG parts of the approximation (deflating the Dirac operator significantly reduces the cost of computing ), and the larger , the lesser the relative cost of the former. The various costs for AMA in our examples are broken down in Table 4 and compared to the high statistics study Yamazaki et al. (2009). In the example using Gaussian sinks, AMA is roughly 16 and 5 times less expensive for roughly the same statistical error, for and 0.01, respectively. LMA is significantly less effective, 3.6 and 2.3 times less expensive. As increases, AMA improves statistics with relatively little extra cost. For instance, for AMA costs an additional 114, in units of the original propagator. The advantage of AMA clearly grows with increasing lattice size and decreasing quark mass. The cost of calculating the correlation functions in this example is negligible, but this may not be the case for more complicated observables. Although disk space and CPU time for eigenvector I/O can be nonnegligible, we ignore these as the costs strongly depend on the implementation details (, we could (de)compress eigenvectors) and the features of the I/O systems used.
Another impressive example of AMA is shown in Fig. 2, which depicts the hadronic vacuum polarization (HVP) from Aubin et al. (2012) and using AMA for roughly the same amount of computational resource (20 configurations, 1400 lowmodes with accuracy , , and sloppy CG stopping residual criterion compared to in Aubin et al. (2012)). The pion mass is MeV and lattice size . The HVP contribution to the muon’s anomalous magnetic moment is sensitive to the low region Aubin et al. (2012), so constraining the HVP in this region is crucial to precisely extract the anomaly. In this test case (which was not optimized), to achieve the same errors on the HVP in the range 01 GeV as the original calculation required about 2.620 times less computer time. Interestingly, LMA actually increases the error in this case by about because the lowmodes do not saturate the WardTakahashi identity. The stopping criterion for can not be too low for the same reason, though our choice may have been too conservative. The costs are summarized in Tab. 4. We note that in this case the cost of constructing the low mode part of the propagator is roughly equivalent to the sloppy CG cost, and that here again the contraction costs are negligible.
LM  Tot.  scaled cost  
, 400 LM  gauss  pt  
AMA  110  1  213  18  91+23  350  0.063  0.065 
LMA  110  1  213  18  23  254  0.279  0.265 
Ref. Yamazaki et al. (2009)  932  4    3728    3728^{1}^{1}1 In Yamazaki et al. (2009) a doubled source was used to reduce this cost by two.  1  1 
, 180 LM  
AMA  158  1  297  74  300+22  693  0.203  0.214 
LMA  158  1  297  74  22  393  0.699  0.937 
Ref. Yamazaki et al. (2009)  356  4    1424    1424  1  1 
HVP  , 1400 LM  max  min  
AMA  20  1  96  11  504+420  1031  0.387  0.050 
LMA  20  1  96  11  420  527  10.3  3.56 
Ref. Aubin et al. (2012)  292  2    584    584  1  1 
In this letter a new class of unbiased error reduction techniques is introduced, using approximations that are covariant under lattice symmetries. This is a generalization of lowmode averaging which reduces the statistical error for observables that are not dominated by lowmodes. We have shown through several numerical examples that allmode averaging is a powerful example of CAA, performing better than LMA and works well even in cases where LMA fails. In the examples given here, AMA reduced the cost by factors up to , over conventional computations, and these factors will only increase for larger lattice sizes and smaller quark masses. The method has great potential for investigations of difficult but important physics problems where statistical fluctuations still dominate the total uncertainty, like the nucleon electric dipole moment or hadronic contributions to the muon anomalous magnetic moment. Since CAA works without introducing any statistical bias (so long as condition appx3 holds), there are many possibilities that also satisfy appx1 and appx2: One can construct using different lattice fermions and parameters (mass, (for DWF), boundary conditions and so on). can be measured on a larger number of gauge configurations, which is potentially advantageous for observables dominated by gauge noise such as disconnected diagrams. One may also consider other types of approximations such as the hopping parameter expansion used in Bali et al. (2010a), or approximations at the level of hadronic Green’s functions.
Acknowledgements.
Numerical calculations were performed using the RICC at RIKEN and the Ds cluster at FNAL. We thank Sinya Aoki, Rudy Arthur, Gunnar Bali, Peter Boyle, Norman Christ, Thomas DeGrand, Leonardo Giusti, Shoji Hashimoto, Tomomi Ishikawa, Chulwoo Jung, Takashi Kaneko, Christoph Lehner, Meifeng Lin, Stefan Schaefer, Ruth Van de Water, Oliver Witzel, Takeshi Yamazaki, Jianglei Yu and other members of JLQCD,RBC,UKQCD collaborations for valuable discussions and comments. CPS QCD library Columbia Physics System code developed and maintained by members of RBC/UKQCD collaborations http://qcdoc.phys.columbia.edu/cps.html () and other softwares (QMP,QIO) are used, which are supported by USQCD and USDOE SciDAC program. This work was supported by the Japanese Ministry of Education GrantinAid, Nos. 22540301 (TI), 23105714 (ES), 23105715 (TI) and U.S. DOE grants DEAC0298CH10886 (TI) and DEFG0292ER40716 (TB). We also thank BNL, the RIKEN BNL Research Center, and USQCD for providing resources necessary for completion of this work.References
 Aubin et al. (2012) C. Aubin, T. Blum, M. Golterman, and S. Peris, (2012), arXiv:1205.3695 [heplat] .
 Yamazaki et al. (2009) T. Yamazaki, Y. Aoki, T. Blum, H.W. Lin, S. Ohta, et al., Phys.Rev. D79, 114505 (2009), arXiv:0904.2039 [heplat] .
 Shintani et al. (2005) E. Shintani, S. Aoki, N. Ishizuka, K. Kanaya, Y. Kikukawa, et al., Phys.Rev. D72, 014504 (2005), arXiv:heplat/0505022 [heplat] .
 Berruto et al. (2006) F. Berruto, T. Blum, K. Orginos, and A. Soni, Phys. Rev. D73, 054509 (2006), arXiv:heplat/0512004 .
 Shintani et al. (2007) E. Shintani, S. Aoki, N. Ishizuka, K. Kanaya, Y. Kikukawa, et al., Phys.Rev. D75, 034507 (2007), arXiv:heplat/0611032 [heplat] .
 Shintani et al. (2008) E. Shintani, S. Aoki, and Y. Kuramashi, Phys. Rev. D78, 014503 (2008), arXiv:0803.0797 [heplat] .
 Blum et al. (2012) T. Blum, P. Boyle, N. Christ, N. Garron, E. Goode, et al., Phys.Rev.Lett. 108, 141601 (2012), arXiv:1111.1699 [heplat] .
 Savage (2012) M. J. Savage, Prog. Part. Nucl. Phys. 67, 140 (2012), arXiv:1110.5943 [nuclth] .
 Giusti et al. (2004) L. Giusti, P. Hernandez, M. Laine, P. Weisz, and H. Wittig, JHEP 04, 013 (2004), arXiv:heplat/0402002 .
 DeGrand and Schaefer (2004) T. A. DeGrand and S. Schaefer, Comput. Phys. Commun. 159, 185 (2004), arXiv:heplat/0401011 .
 Neff et al. (2001) H. Neff, N. Eicker, T. Lippert, J. W. Negele, and K. Schilling, Phys. Rev. D64, 114509 (2001), arXiv:heplat/0106016 .
 Giusti et al. (2003) L. Giusti, C. Hoelbling, M. Luscher, and H. Wittig, Comput. Phys. Commun. 153, 31 (2003), arXiv:heplat/0212012 .
 DeGrand and Schaefer (2005) T. A. DeGrand and S. Schaefer, Phys. Rev. D72, 054503 (2005), arXiv:heplat/0506021 .
 Luscher (2007) M. Luscher, JHEP 07, 081 (2007), arXiv:0706.2298 [heplat] .
 Fukaya et al. (2007) H. Fukaya et al. (JLQCD), Phys. Rev. Lett. 98, 172001 (2007), arXiv:heplat/0702003 .
 Giusti and Necco (2006) L. Giusti and S. Necco, PoS LAT2005, 132 (2006), arXiv:heplat/0510011 .
 Noaki et al. (2008) J. Noaki et al. (JLQCD and TWQCD), Phys. Rev. Lett. 101, 202004 (2008), arXiv:0806.0894 [heplat] .
 Li et al. (2010) A. Li et al. (xQCD), Phys. Rev. D82, 114501 (2010), arXiv:1005.5424 [heplat] .
 Bali et al. (2005) G. S. Bali, H. Neff, T. Duessel, T. Lippert, and K. Schilling (SESAM), Phys. Rev. D71, 114513 (2005), arXiv:heplat/0505012 .
 Foley et al. (2005) J. Foley, K. Jimmy Juge, A. O’Cais, M. Peardon, S. M. Ryan, et al., Comput.Phys.Commun. 172, 145 (2005), arXiv:heplat/0505023 [heplat] .
 Bali et al. (2010a) G. S. Bali, S. Collins, and A. Schafer, Comput. Phys. Commun. 181, 1570 (2010a), arXiv:0910.3970 [heplat] .
 Bali et al. (2010b) G. Bali, L. Castagnini, and S. Collins, PoS LATTICE2010, 096 (2010b), arXiv:1011.1353 [heplat] .
 Luscher (1994) M. Luscher, Nucl.Phys. B418, 637 (1994), arXiv:heplat/9311007 [heplat] .
 Borici and de Forcrand (1995) A. Borici and P. de Forcrand, Nucl.Phys. B454, 645 (1995), arXiv:heplat/9505021 [heplat] .
 Kamleh and Peardon (2012) W. Kamleh and M. Peardon, Comput.Phys.Commun. 183, 1993 (2012), arXiv:1106.5625 [heplat] .
 Aoki et al. (2011) Y. Aoki et al. (RBC), Phys. Rev. D83, 074508 (2011), arXiv:1011.0892 [heplat] .
 (27) Columbia Physics System code developed and maintained by members of RBC/UKQCD collaborations http://qcdoc.phys.columbia.edu/cps.html, .