Frequency-splitting estimators for single-propagator traces
Leonardo Giusti, Tim Harris, Alessandro Nada, Stefan Schaefer
FFrequency-splitting estimators for single-propagatortraces
Leonardo Giusti a , Tim Harris ∗ a , Alessandro Nada b , Stefan Schaefer b a Dipartimento di Fisica, Università di Milano-Bicocca, and INFN, Sezione di Milano-BicoccaPiazza della Scienza 3, I-20126 Milano, Italy b John von Neumann Institute for Computing, DESYPlatanenallee 6, D-15738 Zeuthen, GermanyE-mail: [email protected] , [email protected] , [email protected] , [email protected] In these proceedings we address the computation of quark-line disconnected diagrams in latticeQCD. The evaluation of these diagrams is required for many phenomenologically interesting ob-servables, but suffers from large statistical errors due to the vacuum and random-noise contribu-tions to their variances. Motivated by a theoretical analysis of the variances, we introduce a newfamily of stochastic estimators of single-propagator traces built upon a frequency splitting com-bined with a hopping expansion of the quark propagator, and test their efficiency in two-flavourQCD with pions as light as 190 MeV. The use of these estimators reduces the cost of the com-putation by one to two orders of magnitude over standard estimators depending on the fermionbilinear. As a concrete application, we show the impact of these findings on the computation ofthe hadronic vacuum polarization contribution to the muon anomalous magnetic moment.DESY 20-010 ∗ Speaker. c (cid:13) Copyright owned by the author(s) under the terms of the Creative CommonsAttribution-NonCommercial-NoDerivatives 4.0 International License (CC BY-NC-ND 4.0). https://pos.sissa.it/ a r X i v : . [ h e p - l a t ] J a n requency-splitting estimators Tim Harris
1. Introduction
Quark-line disconnected Wick contractions are ubiquitous in lattice QCD, and those involv-ing single-propagator traces contribute to many interesting observables, such as Standard Modelprocesses like K → ππ , strong and electromagnetic isospin-breaking corrections and isosingletspectroscopy. These observables are particularly computationally challenging as they can havelarge vacuum contributions to their variance, as well as large random noise contributions from theauxiliary fields introduced to evaluate them stochastically [1, 2]. The former can be amelioratedby multi-level integration [3, 4], whereas in the following, we suppress the latter contribution byintroducing a new family of reduced-variance stochastic estimators which are constructed using afrequency splitting and a hopping expansion applied at a large quark mass. These estimators canbe applied in standard Monte Carlo simulations, and here we report on numerical tests with N f = a )-improved Wilson fermions, where we observe between one and two orders of magntitude re-duction in the variance (or cost), depending on the fermion bilinear considered. Further details andany unexplained notation can be found in ref. [5].
2. Variances of single-propagator traces
It is often sufficient to study the variance of individual disconnected Wick contractions, e.g. W and W (assumed to be real), whose product comprises a larger correlation function C W W ( x , x ) = (cid:68)(cid:104) W ( x ) − (cid:104) W ( x ) (cid:105) (cid:105)(cid:104) W ( x ) − (cid:104) W ( x ) (cid:105) (cid:105)(cid:69) . (2.1)For large | x − x | the variance factorizes into the product of the variances of the individual con-tractions, σ C W W ( x , x ) ≈ σ C W ( x ) · σ C W ( x ) + . . . , where σ C Wi ( x i ) = (cid:68)(cid:104) W i ( x i ) − (cid:104) W i ( x i ) (cid:105) (cid:105) (cid:69) , (2.2)and the ellipsis stands for exponentially suppressed terms.A simple example is the disconnected contribution to a two-point function of fermion bilinears,whose disconnected components are the single-propagator traces¯ t Γ , r ( x ) = − a Γ aL ∑ xxx tr (cid:2) Γ D − m r ( x , x ) (cid:3) , (2.3)where D m r is the massive Dirac operator with bare-quark mass m r , a is the lattice spacing and L isthe lattice volume. The factor a Γ = − i for Γ = γ µ and a Γ = Γ = I , γ , γ µ γ , σ µν , is chosen sothat the Wick contraction is real. The gauge variance can be defined in terms of local operators σ t Γ , r = a Γ L ∑ x a (cid:104) O Γ , rr ( , x ) O Γ , r (cid:48) r (cid:48) ( ) (cid:105) c , (2.4)where c stands for connected correlation function, O Γ , rs ( x ) = ¯ ψ r ( x ) Γ ψ s ( x ) , and m r (cid:48) = m r . It isevident that the gauge variance is itself a disconnected contraction, and so begins only at order g orhigher in perturbation theory and may be expected to be suppressed. Nevertheless, by the operatorproduct expansion and power-counting the variance has a cubic divergence in the continuum limit.1 requency-splitting estimators Tim Harrisid L / a κ MDU N cfg M π [MeV] M π L F7 48 0 . ( )
268 4 . . . Table 1:
Overview of the ensembles and statistics presented in this study and their simulation and physicsparameters. − − − − − σ · ( L a ) N s SPT jk A k V k Figure 1:
Variances of the standard random noise estimators defined in eq. (2.6) as a function of N s for theensemble F7. The symbols S , P , T jk , A k and V k stand for Γ = I , γ , σ jk , γ k γ and γ k respectively. The dashedlines indicate the gauge noise contributions to the variances computed in sec. 4. In practice, it is unfeasible to compute the single-propagator trace exactly, to which end weintroduce N s independent auxiliary fields η i ( x ) , whose components must have unit variance andzero mean, which we choose to be drawn from a Gaussian distribution, and a stochastic estimator¯ τ Γ , r ( x ) = − aL N s N s ∑ i = ∑ xxx Re (cid:104) a Γ η † i ( x ) Γ { D − m r η i } ( x ) (cid:105) . (2.5)The variance of the stochastic estimator receives contributions from the fluctuations of the auxiliaryfields σ τ Γ , r = σ t Γ , r − L N s (cid:26) a Γ ∑ x a (cid:104) O Γ , rr (cid:48) ( , x ) O Γ , r (cid:48) r ( ) (cid:105) + a ∑ x a (cid:104) P rr (cid:48) ( x ) P r (cid:48) r ( ) (cid:105) (cid:27) , (2.6)where P rs = O γ , rs . The second and third terms in eq. (2.6) are connected diagrams which occur attree level in perturbation theory, which suggests that they will be larger than the gauge varianceunless N s is large. Note that the third term, which is chirally-enhanced, is independent of Γ .In order to investigate the relative magnitude of the contributions to the variance, we employedthe ensembles listed in table 1, with N f = a )-improved Wilsonfermions. In fig. 1 we plot the variance as a function of N s for the scalar and pseudoscalar densities,and the tensor, axial and vector currents for the F7 ensemble. The linear behaviour in N − s illustratesthat the random noise variance dominates by orders of magnitude for small N s , and is practicallyindependent of the bilinear, as expected from the arguments outlined above. For large N s , thegauge variance (dashed lines) is saturated, and is orders of magnitude larger for the densities thanthe currents. It is therefore highly desirable to examine variance reduction methods in particularfor the currents in order to reach the gauge noise.2 requency-splitting estimators Tim Harris
3. Estimators for differences of single-propagator traces
In this section, we investigate the difference of single-propagator traces ¯ t Γ , rs = ¯ t Γ , r − ¯ t Γ , s withtwo different bare quark masses, m r < m s , which contains the infrared contributions to the single-propagator trace. Such differences arise, for example, from the disconnected contraction of theelectromagnetic current in the isosymmetric theory, as well as constituting a component of ourimproved frequency-splitting estimator in sec. 4.We examine two estimators for the differences of single-propagator traces defined by ¯ θ Γ , rs ( x ) = − ( m s − m r ) aL N s N s ∑ i = ∑ xxx Re (cid:104) a Γ η † i ( x ) Γ { D − m r D − m s η i } ( x ) (cid:105) , (3.1)¯ τ Γ , rs ( x ) = − ( m s − m r ) aL N s N s ∑ i = ∑ xxx Re (cid:104) a Γ { η † i D − m r } ( x ) Γ { D − m s η i } ( x ) (cid:105) , (3.2)which we denote as the standard [6] and split-even estimators [5] in the following. Note that thecost of the two estimators is the same for a given N s . Analogously to the previous section, theirvariances can be defined in terms of local operators σ θ Γ , rs = σ t Γ , rs − ( m s − m r ) L N s (cid:40) a Γ ∑ y , y , y a (cid:104) S rs ( y ) O Γ , ss (cid:48) ( , y ) S s (cid:48) r (cid:48) ( y ) O Γ , r (cid:48) r ( ) (cid:105) + a ∑ y , y , y a (cid:104) S rs ( y ) P ss (cid:48) ( y ) S s (cid:48) r (cid:48) ( y ) P r (cid:48) r ( ) (cid:105) (cid:41) , (3.3) σ τ Γ , rs = σ t Γ , rs − a Γ ( m s − m r ) L N s ∑ y , y , y a (cid:110)(cid:10) S rs ( y ) O Γ , ss (cid:48) ( , y ) S s (cid:48) r (cid:48) ( y ) O Γ , r (cid:48) r ( ) (cid:11) + (cid:10) P rr (cid:48) ( y ) O Γ , r (cid:48) s (cid:48) ( , y ) P s (cid:48) s ( y ) O Γ , sr ( ) (cid:11)(cid:111) (3.4)where S rs = O I , rs and the gauge variance σ t Γ , rs is the fully-connected analogue of the first four-pointfunction in eq. (3.3), just as for the single-propagator trace. The two estimators differ only betweenthe second four-point functions in eq. (3.3) and eq. (3.4). For the standard estimator, the (cid:104) SPSP (cid:105) term is independent of Γ and integrated over one more time coordinate compared to the (cid:104) POPO (cid:105) term in the split-even estimator, which furthermore depends on Γ .In fig. 2, we compare the variances of the standard (filled) and split-even (open) estimatorsfor the pseudoscalar and vector channels for the difference of light- and strange-quark propagatorswith bare-quark masses am q , r = . am q , s = . N s ∼ O ( ) for the pseudoscalar channel, and O(100) in the vectorchannel.Interestingly, the large gain using the split-even estimator reported here has been confirmedin ref. [7]. The partial cancellation of stochastic noise between the light- and strange-quark tracesis already present in the baseline standard estimator, which is not the origin of the significant gainas suggested there. The large reduction in the variance is well explained instead by the preceding Here, we have used the identity D − m r − D − m s = ( m s − m r ) D − m r D − m s to rewrite the difference as a product. requency-splitting estimators Tim Harris − − − − − σ · ( L a ) N s P standardsplit-even 10 − − − − − σ · ( L a ) N s V k standardsplit-even Figure 2:
Variances of the standard θ Γ , rs (filled symbols) and the split-even τ Γ , rs (open symbols) estimatorsfor the pseudoscalar density (left) and the vector current (right) for ensemble F7. -1012345 0 5 10 15 20 25 30 x /aa C rsV V × split-even -200-1000100200300400500600700 0 10 20 30 40 a µ ( x c u t ) × − x cut0 /a (ud) conn.(uds) disc. × − Figure 3:
Left: the disconnected contribution to the electromagnetic current correlator using the split-even estimator from N cfg = N s =
512 noise sources for F7, and (right)the corresponding contribution to a µ from distances up to x cut0 . A striking plateau in x cut0 is visible in theconnected case, but no evident for the disconnected contribution and moreover, the error grows quickly dueto the vacuum contributions to the gauge noise, which can be tackled with multi-level integration. formulæ for the variances of the two estimators. This analysis also explains the origin of the empir-ical gains observed for the one-end trick for the pseudoscalar density in twisted-mass QCD [8], inwhich case it is even possible to show the estimator has a strictly smaller variance than the standardone [5].As an application, we consider the disconnected contribution to the electromagnetic currentcorrelator for an isodoublet of light quarks and a valence strange quark, C rsVV ( x ) ∼ (cid:104) ¯ t γ k , rs ( x + y ) ¯ t γ k , rs ( y ) (cid:105) . The disconnected contribution is shown in fig. 3 (left) using N cfg = a µ , and in fig. 3 (right) we show the resulting disconnected contribution (timesa factor 100, for visibility) computed using the split-even estimator, along with the light-quark4 requency-splitting estimators Tim Harris connected contribution, from distances up to x cut0 using the time-momentum representation [9]. Inparticular, this variance-reduction technique represents an important computational advance forprecision computations of hadronic contributions to muon g −
2, and as suggested in ref. [5], theuse of this estimator can significantly speed up many other computations such as of hadronic matrixelements and their strong and electromagnetic isospin-breaking corrections, see ref. [10] for a firstapplication to isospin-breaking corrections.
4. Frequency-splitting estimators
In order to construct a reduced-variance estimator for the single-propagator trace, we combinethe split-even estimator and for the large quark-mass contribution use the order-2 n hopping expan-sion of the Dirac operator [11, 12], D − m r = M n , m r + D − m r H nm r , where H m r denotes the hopping partof the Dirac operator and M n , m r the first 2 n − t Γ , r = ¯ t M Γ , r + ¯ t R Γ , r . The first term can be computed exactlywith 24 n applications of M n , m r onto probing vectors, while an efficient stochastic estimator forthe remainder of the hopping expansion is given by¯ τ R Γ , r ( x ) = − aL N s N s ∑ i = ∑ x Re (cid:110) a Γ (cid:2) η † i H nm r (cid:3) ( x ) Γ (cid:2) D − m r H nm r η i (cid:3) ( x ) (cid:111) . (4.1)We therefore define the frequency-splitting estimator for the target quark mass m by using thesplit-even estimator for K − m r k < m r k + and applying the hopping decompositonat the largest quark mass m r K which controls the ultraviolet fluctuations,¯ τ fs Γ , r ( x ) = K − ∑ k = ¯ τ Γ , rkrk + ( x ) + ¯ t M Γ , r K ( x ) + ¯ τ R Γ , r K ( x ) , (4.2)In fig. 4, we investigate two frequency-splitting estimators for single-propagator traces at the sea-quark mass for the G8 (left) and F7 (right) ensemble respectively. For G8 we use K = am r K = .
1, while for F7 we take K = am r K = .
3, and n = . In both cases, we see around two orders of magnitude reduction in the random-noisecontributions to the variance, but due to the estimated increase of about 3 . −
30 depending on the mass.
5. Conclusions
In these proceedings, we have investigated improved split-even estimators for differences ofsingle-propagator traces, which can be applied to efficiently evaluate the disconnected contribu-tions arising from the electromagnetic current, so that the leading contribution to the variancesarises from fluctuations of the gauge field. The cost is reduced by between one to two ordersof magnitude depending on the bilinear. Furthermore, these can be used to construct frequency-splitting estimators for single-propagator traces by combining them with the hopping expansionfor large quark masses, which reduces the cost by more than an order of magnitude for the vectorcurrent close to the physical point. These techniques are compatible with other variance-reductiontechniques, such as low-mode averaging [13, 14] and dilution [15]. See ref. [5] for the relative N s used in each component. requency-splitting estimators Tim Harris − − − − σ · ( L a ) N s V k standardFS1 10 − − − − σ · ( L a ) N s V k standardFS2 Figure 4:
Variances of the frequency-splitting estimators (open symbols) for the vector current, comparedwith the standard random-noise estimators (filled symbols) on the G8 ensemble (left) and F7 ensemble(right). The cost of one iteration of the FS estimator is about 3.3 and 6 times the standard one on F7 and G8respectively.
Acknowledgments
Simulations have been performed on the PC clusters Marconi at CINECA (CINECA-INFN and CINECA-Bicocca agreements) and Wilson at Milano-Bicocca. We are grateful to our colleagueswithin the CLS initiative for sharing the ensembles of gauge configurations with two dynamical flavours.L.G. and T. H. acknowledge partial support by the INFN project “High performance data network”.
References [1] K. Bitar et al.
Nucl. Phys.
B313 (1989), pp. 348–376.[2] C. Michael and J. Peisa.
Phys. Rev.
D58 (1998), p. 034506.[3] M. Cè et al.
Phys. Rev.
D93.9 (2016), p. 094507.[4] M. Cè et al.
Phys. Rev.
D95.3 (2017), p. 034503.[5] L. Giusti et al.
Eur. Phys. J.
C79.7 (2019), p. 586.[6] V. Gülpers et al.
PoS
LATTICE2014 (2014), p. 128.[7] H. Wittig et al.
FCCP2019 Anacapri, Capri Island, Italy, August 29-31, 2019 . 2019.[8] P. Boucaud et al.
Comput. Phys. Commun.
179 (2008), pp. 695–715.[9] D. Bernecker and H. B. Meyer.
Eur. Phys. J.
A47 (2011), p. 148.[10] A. Risch and H. Wittig.
Lattice 2019 Wuhan, Hubei, China, June 16-22, 2019 . 2019.[11] G. S. Bali et al.
Comput. Phys. Commun.
181 (2010), pp. 1570–1583.[12] V. Gülpers et al.
Phys. Rev.
D89.9 (2014), p. 094503.[13] L. Giusti et al.
JHEP
04 (2004), p. 013.[14] T. A. DeGrand and S. Schaefer.
Comput. Phys. Commun.
159 (2004), pp. 185–191.[15] J. Foley et al.
Comput. Phys. Commun.
172 (2005), pp. 145–162.172 (2005), pp. 145–162.