Several important observables, like the quark condensate and the Taylor
coefficients of the expansion of the QCD pressure with respect to the chemical
potential, are based on the trace of the inverse Dirac operator and of its
powers. Such traces are traditionally estimated with "noise vectors"
sandwiching the operator. We explore alternative approaches based on polynomial
approximations of the inverse Dirac operator.Comment: Eight pages, thirteen figures. Proceedings of the 35th International
  Symposium on Lattice Field Theor

de Forcrand, Philippe

Jaeger, Benjamin

EPJ Web of Conferences

English

arXiv

Several important observables, like the quark condensate and the Taylor coefficients of the expansion of the QCD pressure with respect to the chemical potential, are based on the trace of the inverse Dirac operator and of its powers. Such traces are traditionally estimated with “noise vectors” sandwiching the operator. We explore alternative approaches based on polynomial approximations of the inverse Dirac operator

Philippe de Forcrand

Benjamin Jäger

EDP Sciences OAI-PMH repository (1.2.0)

Alternatives to the stochastic “noise vector” approach

Several important observables, like the quark condensate and the Taylor coefficients of the expansion of the QCD pressure with respect to the chemical potential, are based on the trace of the inverse Dirac operator and of its powers. Such traces are traditionally estimated with "noise vectors" sandwiching the operator. We explore alternative approaches based on polynomial approximations of the inverse Dirac operator

Forcrand, Philippe de

Syddansk Universitets Forskerportal

       University of Southern DenmarkAlternatives to the stochastic "noise vector" approachForcrand, Philippe de; Jaeger, BenjaminPublished in:EPJ Web of ConferencesDOI:10.1051/epjconf/201817514022Publication date:2018Document version:Submitted manuscriptCitation for pulished version (APA):Forcrand, P. D., & Jaeger, B. (2018). Alternatives to the stochastic "noise vector" approach. EPJ Web ofConferences, 175, Article 14022. https://doi.org/10.1051/epjconf/201817514022Go to publication entry in University of Southern Denmark's Research PortalTerms of useThis work is brought to you by the University of Southern Denmark.Unless otherwise specified it has been shared according to the terms for self-archiving.If no other license is stated, these terms apply:            • You may download this work for personal use only.            • You may not further distribute the material or use it for any profit-making activity or commercial gain            • You may freely distribute the URL identifying this open access versionIf you believe that this document breaches copyright please contact us providing details and we will investigate your claim.Please direct all enquiries to puresupport@bib.sdu.dkDownload date: 21. Oct. 2025CERN-TH-2017-210CP3-Origins-2017-44 DNRF90Alternatives to the stochastic “noise vector” approachPhilippe de Forcrand1,2,? and Benjamin Jäger1,31Institute for Theoretical Physics, ETH, CH-8093 Zürich, Switzerland2CERN, Theory Division, CH-1211 Geneva, Switzerland3CP3-Origins & Danish Institute for Advanced Study, Department of Mathematics and Computer Science,University of Southern Denmark, 5230 Odense M, DenmarkAbstract. Several important observables, like the quark condensate and the Taylor coef-ficients of the expansion of the QCD pressure with respect to the chemical potential, arebased on the trace of the inverse Dirac operator and of its powers. Such traces are tradi-tionally estimated with "noise vectors" sandwiching the operator. We explore alternativeapproaches based on polynomial approximations of the inverse Dirac operator.1 IntroductionMeasuring observables on an ensemble of gauge field configurations can require computing resourcessimilar to those required to generate the configurations. However, comparatively little effort hasbeen devoted to optimizing the computing strategy for such measurements. An important categoryof fermionic observables can be expressed as the trace of some function of the Dirac operator, orsome other related operator. For example, the quark condensate ψ̄ψ is obtained from the trace ofthe inverse of the Dirac operator D/ . Generically, the trace for a given gauge field configuration isestimated by using “noise vectors” ηk, k = 1, .., n, viaTrA ≈1nn∑k=1η†k Aηk (1)with 〈ηi†k ηjk〉 = δi j. Here, i and j label components of a noise vector η, and 〈..〉 means averaging overnoise vectors (for a given gauge field). Clearly, there are two sources of noise here: the “stochasticnoise” coming from the limited number n of noise vectors used on a given gauge field configuration;and the “gauge noise” coming from the fluctuations caused by the variation of the gauge field fromconfiguration to configuration. The balance between these two types of noise depends on the observ-able, and on the parameters of the Monte Carlo simulation: in particular, fluctuations coming fromthe gauge field will severely increase in the vicinity of a phase transition.Here, we consider two situations of relevance at finite temperature:(i) the measurement of the chiral condensate and its first 4 moments near the finite-temperature de-confinement transition. The moments are necessary to construct the Binder cumulant (or kurtosis),whose value at the phase transition is known by universality;?e-mail: forcrand@phys.ethz.charXiv:1710.07305v1  [hep-lat]  19 Oct 2017(ii) the measurement of the Taylor coefficients of the expansion of the pressure in µ/T .In the latter situation, it is not unusual to take 1000 noise vectors or more per gauge configuration.Can one be more efficient?2 Related recent progressApplied mathematicians have proposed a novel approach to determine the spectrum of a Hermitianmatrix A [1]. Instead of obtaining the eigenvalues directly by a Krylov-based method, the idea is todetermine the “eigenvalue count”, namely the number µ(a, b) of eigenvalues in some interval [a, b].Upon varying a and b, the positions of the jumps in the eigenvalue count give the desired eigenvalues.Thus, one needs a sliding bandpass filter f (x, a, b) which is = 1 if x ∈ [a, b], 0 otherwise. In practice:- The function f is approximated by Chebyshev polynomials and applied to the matrix A. Then,the “eigenvalue count” µ(a, b) is Tr f (A, a, b).- The trace is approximated by using a small number of noise vectors following eq.(1).One great advantage is that, by varying the coefficients of the polynomial, a and b can be varied atwill without any additional matrix-vector multiplication, and thereby the complete spectrum of A canbe determined at low cost.Two lattice groups have tested this approach so far:• In [2], a degree-8000 Chebyshev approximation was constructed and applied to D/ †D/ , using 20noise vectors for the trace. This required 80000 Dirac-matrix-vector multiplications. These are verymoderate resources, considering that the lattice size was 643 × 96, and that the whole spectrum couldbe extracted. After averaging over 10 gauge configurations, a precise determination of 〈ψ̄ψ〉 could beobtained.• In [3], an overlap Dirac operator was used, on the same lattice size, with a Chebyshev approximationof the same degree 8000, and 1 noise vector per configuration only. The gauge ensemble consisted of50 configurations. Again, a precise determination of 〈ψ̄ψ〉 could be obtained.These promising first results were obtained at zero temperature. We are interested in applying asimilar approach at T ∼ Tc, the deconfinement temperature, where the “gauge noise” is expected to bemuch greater. And we want, in addition to ψ̄ψ =∑i λ−1i , higher negative moments of the eigenvalues,i.e.∑i λ−mi ,m = 1, 2, 3, 4, ...3 Application I: Binder cumulant for the finite-temperature transitionA suitable order parameter for the finite-T transition is the quark condensate ψ̄ψ. To probe the orderof the phase transition, it is customary to measure its Binder cumulant (or equivalently its kurtosis)B4(ψ̄ψ) ≡〈(δψ̄ψ)4〉〈(δψ̄ψ)2〉2, δψ̄ψ = ψ̄ψ − 〈ψ̄ψ〉 (2)Thus, the first 4 moments (ψ̄ψ)k, k = 1, .., 4 must be measured, or estimated in an unbiased way. Thelatter approach is normally chosen, using noise vectors. One must then decide how many noise vectorsper configuration to use: using too few will increase the “stochastic noise” on B4 and the associatedstatistical error; using too many is wasteful of computing resources.Here, we chose the most economical approach: 4 noise vectors η1, .., η4, forming unbiased estimators 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0  0.2  0.4  0.6  0.8  1  1.2  1.4Estimator pbpExact pbp 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0  0.2  0.4  0.6  0.8  1  1.2  1.4Estimator pbp2Exact pbp2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0  0.2  0.4  0.6  0.8  1  1.2  1.4Estimator pbp3Exact pbp3 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0  0.2  0.4  0.6  0.8  1  1.2  1.4Estimator pbp4Exact pbp4 0 100 200 300 400 500 600 0.4  0.6  0.8  1  1.2Distribution of psibarpsiExactNoiseFigure 1. Comparison of the first 4 powers of ψ̄ψ with their stochastic estimators eq.(3), on an ensemble of 63×4SU(3),N f = 4 configurations with staggered fermions, near the finite-temperature deconfinement transition. Thebottom figure compares the resulting distribution of the exactly measured ψ̄ψ with that of its estimator: thetwo-bump distribution signals the vicinity of the finite-T transition; the "stochastic noise" barely increases thevariance of the distribution.asψ̄ψ ←−14(η†1D/−1η1 + η†2D/−1η2 + . . .)(ψ̄ψ)2 ←−16(η†1D/−1η1 · η†2D/−1η2 + . . .)(ψ̄ψ)3 ←−14(η†1D/−1η1 · η†2D/−1η2 · η†3D/−1η3 + . . .)(ψ̄ψ)4 ←−(η†1D/−1η1 · η†2D/−1η2 + η†3D/−1η3 · η†4D/−1η4)(3)In Fig. 1, we compare the estimators above with the exact values obtained by diagonalization of D/ ,for an ensemble of 63 × 4 SU(3) configurations, where the parameters (N f = 4, am = 0.07) have beenchosen to be in the vicinity of a sharp temperature-driven crossover. Indeed, the exact distributionof ψ̄ψ in the last panel shows a clear two-peak distribution. What is striking is that the distributionof the stochastic estimator of ψ̄ψ, in the same figure, is hardly broader, even though the number ofnoise vectors has been kept to its minimum value 4. Thus, in this example, the “stochastic noise” isnegligible compared to the “gauge noise”. Whether or not this persists on larger lattices needs to bechecked.0.2 0.4 0.6 0.8 1.0Chebychevapprox.-1.0 -0.5 0.5 1.0200400600800μ(-1,x)Figure 2. Chebyshev bandpass filter (left) and eigenvalue count µ(−1, x) as a function of x, for the rescaled,staggered operator (iD/ ) on a 44 lattice.The novel method outlined in Sec. 2 can give us nearly exact results at a potentially reduced cost,compared with exact diagonalization. A Chebyshev bandpass filter is applied to (iD/ ), so that thenumber µ of eigenvalues in [a, b] can be approximated asµ(a, b) =1nn∑j=1pmax∑p=0gp(a, b) 〈〈η†jTp(iD/ ) η j〉〉 (4)where Tp is the Chebyshev polynomial of degree p, gp are the coefficients of the filter approximation,and 〈〈· · · 〉〉 represents the averaging over the gauge field. Fig. 2 illustrates the bandpass filter (left) andthe eigenvalue count µ(−1, x) versus x (right) 1, for one 44 configuration (D/ has been rescaled so thatits spectrum is contained in [−1,+1]). µ(−1, x) is a piecewise constant function, increasing by one asx passes each eigenvalue. We approximate the eigenvalues as solutions of µ(−1, x) being half-integer.Since all eigenvalues are obtained in this way, we can form arbitrary moments Tr(D/ −k), k = 1, 2, . . ..This is illustrated in Table 1 below. Note the remarkable accuracy, even for high moments.k exact approx. approx./exact2 −453.767 −451.148 0.9942294 3871.66 3844.25 0.9929210 −8.136 · 107 −8.091 · 107 0.994547Table 1. Tr(D/ −k)using 4 noise vectors and a degree-16k polynomial on a 44 lattice.4 Application II: Taylor expansion of free energyQCD at non-zero quark chemical potential µ represents a long-standing challenge for lattice simu-lations: the Dirac determinant becomes complex when µ , 0, which prevents its interpretation as asampling probability. Several strategies have been proposed to circumvent this “sign” (or phase) prob-lem. One of the first is that of Taylor-expanding the QCD free energy, or the pressure P, in powers of(µ/T ) around µ = 0:P(µ,T ) = P(µ = 0,T ) +∑k=1ck(T )(µT)k. (5)Charge conjugation symmetry implies that odd-order coefficients ck vanish. However, products of oddderivatives do contribute in the presence of multiple quark flavors, each with its chemical potential.Now, ck is simply the kth derivative of the free energy evaluated at µ = 0, so it can be measured in con-ventional, sign-problem free simulations. The µ-dependence appears only in the Dirac determinant,so that the generic expression for ck is the trace of a polynomial expression in D/ −1 and in derivatives∂mD/ /∂µm. The number of terms grows exponentially with the order k, with positive or negative signs.This makes the calculation of high-order coefficients ck extremely challenging: no calculation hasever yielded information at 10th order. Improvements in efficiency are needed.One such improvement has been proposed recently by Gavai and Sharma [5, 6]. The suggestionis to change the way the chemical potential is introduced on the lattice, and adopt a naive, linearprescription rather than the exponential prescription of Hasenfratz and Karsch [7]. This is rathercounter-intuitive, since it re-introduces a UV divergence in the pressure. Gavai and Sharma argue thatthis divergence can be removed by a numerical subtraction. The gain is that derivatives ∂mD/ /∂µmnow vanish for m > 1. This greatly simplifies the structure of the ck’s, which become the traces ofhomogeneous polynomials of degree k in A−1, with A = D/(∂D/∂µ)−1(note that A is neither hermitian,nor anti-hermitian – see Fig. 3).At this stage, our approach may become advantageous. Suppose we can determine the [complex]spectrum of A. Then, ck will be obtained as a sum of products of factors Tr (A−m), which can all bedetermined in one go!1Note the similarity with Lüscher’s “mode number” [4].-4 -2 2 4-4-2244x4³-4 -2 2 4-4-2244x6³-4 -2 2 4-4-2244x8³Figure 3. Complex spectrum of(D/(∂D/∂µ)−1), for lattices of increasing spatial size near the finite-T transition.Here, we take first steps in that promising direction. First, we generalize the notion of eigenvaluecount to the complex plane. The number of eigenvalues µ(Γ) inside a complex contour Γ can beobtained by a Cauchy integral [8]:µ(Γ) =12πi∮Γdz Tr[(z1 − A)−1]=12πi∮Γdz∑j1z − λ j(6)For Γ, we choose a circle of center c and radius r, and discretize the integral using the trapezoidal rule:µ(r, c) ∼1NN−1∑k=0r e2πiN k Tr[(zk1 − A)−1]with zk = c + r e2πiN k (7)where the trace is estimated using noise vectors. The so-obtained value of µ(r, c) will be nearly aninteger equal to the number of eigenvalues inside the circle. One can now discard the circles whichcontain zero eigenvalue, and divide by 2 the radius of non-empty circles. This procedure is repeatedrecursively, until each circle contains one and only one eigenvalue, and has a radius smaller than someprescribed tolerance ε. In this fashion, the complete spectrum of A is determined. Arbitrary momentsTrA−m can be evaluated.The crucial point is that the estimation of the trace in eq.(7) proceeds by solving linear systems(zk1 − A) x = η, using a multi-shift linear solver. Then the cost of solving many such systems isnot much greater than that of solving just one, because these systems differ only by the shift zk, andthe Krylov space built by the linear solver is shift-invariant. In addition, each iteration of the solverinvolves matrix-vector multiplication by A = D/(∂D/∂µ)−1. At first sight, the inversion of ∂D/∂µlooksdaunting. But this operator is completely local in space, so that the corresponding matrix is made ofsmall diagonal blocks, each of which is easy to invert.We illustrate our approach in Figs. 3 and 4. Fig. 3 shows the complex spectrum of A, for 3configurations obtained with the same action, near the finite-T transition, on lattices of increasingspatial sizes.Fig. 4 shows 3 successive stages of our recursive procedure: at each stage, circles empty of eigen-values are discarded, and non-empty circles are replaced by 4 circles with half the radius. The finaloutput is an estimate of all eigenvalues of A, from which successive moments Tr(D/ ∂D/∂µ−1)−mcan beestimated as well. Illustrative results are provided in Table 2 below.-4 -2 2 4-4-224-4 -2 2 4-4-224-4 -2 2 4-4-224Figure 4. Circular contours are recursively refined until each contour encloses one eigenvalue (or a few).m exact approx. ratio1 4.10306i 0.0352058 + 4.10474i 1.000412 142.497 142.607 − 0.00748385i 1.000773 10.8641i 0.04117 + 10.8765i 1.001144 −14.4443 −14.3171 − 0.0751784i 0.99119510 −133.716 −135.113 − 3.93177i 1.0104520 −79320.5 −81511.3 + 979.464i 1.02762Table 2. Tr (A−m) for 100 quadrature points and 100 noise vectors for ε = 10−3, on a 44 gauge configuration.5 OutlookWe have described a novel approach to obtain estimates of moments Tr (A−m) of a sparse matrix A,in the Hermitian and in the non-Hermitian case. Two applications were tested: measuring the Bindercumulant of the quark condensate (A = iD/ ), and the coefficients of the Taylor expansion of thepressure in µ/T (A = D/(∂D/∂µ)−1). In both cases, the regime of interest is near the finite-temperaturedeconfinement transition. In this regime, fluctuations of the gauge field are large, so that the “gaugenoise” may well exceed the “stochastic noise” caused by the limited number of noise vectors used toestimate Tr (A−m). Indeed, this is what happens in the first application: it is then sufficient to use 4noise vectors per configuration, at least for the small lattice sizes considered in this test.In the second application, we have generalized the eigenvalue count to the complex plane. Theviability of our approach rests on a massive usage of non-Hermitian multi-shift linear solvers. Greatsimplifications in the contour integrals are possible and are under study. However, while we obtainnearly exact results for the first 20 moments Tr (A−m), we have not measured the “gauge noise”,which presumably increases quickly with m. This may represent the ultimate, irreducible origin of thecomputing cost of the Taylor expansion of the pressure, as foreseen in [9].6 AcknowledgementsBJ was supported by the Schweizerischer Nationalfonds (SNF) under grant 200020-162515. PdFthanks the Kavli Institute for Theoretical Physics for hospitality. This research was supported in partby the National Science Foundation under Grant No. NSF PHY-1125915.⇒ ⇒References[1] E. Di Napoli, E. Polizzi, Y. Saad, Numerical Linear Algebra with Applications 23, 674 (2016)[2] Z. Fodor, K. Holland, J. Kuti, S. Mondal, D. Nogradi, C.H. Wong, PoS LATTICE2015, 310(2016), 1605.08091[3] G. Cossu, H. Fukaya, S. Hashimoto, T. Kaneko, J.I. Noaki, PTEP 2016, 093B06 (2016),1607.01099[4] L. Giusti, M. Lüscher, JHEP 03, 013 (2009), 0812.3638[5] R.V. Gavai, S. Sharma, Phys. Lett. B749, 8 (2015), 1406.0474[6] R.V. Gavai, S. Sharma, PoS CPOD2014, 015 (2015), 1503.03575[7] P. Hasenfratz, F. Karsch, Physics Letters B 125, 308 (1983)[8] Y. Futamura, H. Tadano, T. Sakurai, JSIAM Letters 2, 127 (2010)[9] P. de Forcrand, PoS LAT2009, 010 (2009), 1005.0539

Alternatives to the stochastic "noise vector" approach

       University of Southern DenmarkAlternatives to the stochastic "noise vector" approachForcrand, Philippe de; Jaeger, BenjaminPublished in:EPJ Web of ConferencesDOI:10.1051/epjconf/201817514022Publication date:2018Document version:Submitted manuscriptCitation for pulished version (APA):Forcrand, P. D., & Jaeger, B. (2018). Alternatives to the stochastic "noise vector" approach. EPJ Web ofConferences, 175, Article 14022. https://doi.org/10.1051/epjconf/201817514022Go to publication entry in University of Southern Denmark's Research PortalTerms of useThis work is brought to you by the University of Southern Denmark.Unless otherwise specified it has been shared according to the terms for self-archiving.If no other license is stated, these terms apply:            • You may download this work for personal use only.            • You may not further distribute the material or use it for any profit-making activity or commercial gain            • You may freely distribute the URL identifying this open access versionIf you believe that this document breaches copyright please contact us providing details and we will investigate your claim.Please direct all enquiries to puresupport@bib.sdu.dkDownload date: 23. Oct. 2025CERN-TH-2017-210CP3-Origins-2017-44 DNRF90Alternatives to the stochastic “noise vector” approachPhilippe de Forcrand1,2,? and Benjamin Jäger1,31Institute for Theoretical Physics, ETH, CH-8093 Zürich, Switzerland2CERN, Theory Division, CH-1211 Geneva, Switzerland3CP3-Origins & Danish Institute for Advanced Study, Department of Mathematics and Computer Science,University of Southern Denmark, 5230 Odense M, DenmarkAbstract. Several important observables, like the quark condensate and the Taylor coef-ficients of the expansion of the QCD pressure with respect to the chemical potential, arebased on the trace of the inverse Dirac operator and of its powers. Such traces are tradi-tionally estimated with "noise vectors" sandwiching the operator. We explore alternativeapproaches based on polynomial approximations of the inverse Dirac operator.1 IntroductionMeasuring observables on an ensemble of gauge field configurations can require computing resourcessimilar to those required to generate the configurations. However, comparatively little effort hasbeen devoted to optimizing the computing strategy for such measurements. An important categoryof fermionic observables can be expressed as the trace of some function of the Dirac operator, orsome other related operator. For example, the quark condensate ψ̄ψ is obtained from the trace ofthe inverse of the Dirac operator D/ . Generically, the trace for a given gauge field configuration isestimated by using “noise vectors” ηk, k = 1, .., n, viaTrA ≈1nn∑k=1η†k Aηk (1)with 〈ηi†k ηjk〉 = δi j. Here, i and j label components of a noise vector η, and 〈..〉 means averaging overnoise vectors (for a given gauge field). Clearly, there are two sources of noise here: the “stochasticnoise” coming from the limited number n of noise vectors used on a given gauge field configuration;and the “gauge noise” coming from the fluctuations caused by the variation of the gauge field fromconfiguration to configuration. The balance between these two types of noise depends on the observ-able, and on the parameters of the Monte Carlo simulation: in particular, fluctuations coming fromthe gauge field will severely increase in the vicinity of a phase transition.Here, we consider two situations of relevance at finite temperature:(i) the measurement of the chiral condensate and its first 4 moments near the finite-temperature de-confinement transition. The moments are necessary to construct the Binder cumulant (or kurtosis),whose value at the phase transition is known by universality;?e-mail: forcrand@phys.ethz.charXiv:1710.07305v1  [hep-lat]  19 Oct 2017(ii) the measurement of the Taylor coefficients of the expansion of the pressure in µ/T .In the latter situation, it is not unusual to take 1000 noise vectors or more per gauge configuration.Can one be more efficient?2 Related recent progressApplied mathematicians have proposed a novel approach to determine the spectrum of a Hermitianmatrix A [1]. Instead of obtaining the eigenvalues directly by a Krylov-based method, the idea is todetermine the “eigenvalue count”, namely the number µ(a, b) of eigenvalues in some interval [a, b].Upon varying a and b, the positions of the jumps in the eigenvalue count give the desired eigenvalues.Thus, one needs a sliding bandpass filter f (x, a, b) which is = 1 if x ∈ [a, b], 0 otherwise. In practice:- The function f is approximated by Chebyshev polynomials and applied to the matrix A. Then,the “eigenvalue count” µ(a, b) is Tr f (A, a, b).- The trace is approximated by using a small number of noise vectors following eq.(1).One great advantage is that, by varying the coefficients of the polynomial, a and b can be varied atwill without any additional matrix-vector multiplication, and thereby the complete spectrum of A canbe determined at low cost.Two lattice groups have tested this approach so far:• In [2], a degree-8000 Chebyshev approximation was constructed and applied to D/ †D/ , using 20noise vectors for the trace. This required 80000 Dirac-matrix-vector multiplications. These are verymoderate resources, considering that the lattice size was 643 × 96, and that the whole spectrum couldbe extracted. After averaging over 10 gauge configurations, a precise determination of 〈ψ̄ψ〉 could beobtained.• In [3], an overlap Dirac operator was used, on the same lattice size, with a Chebyshev approximationof the same degree 8000, and 1 noise vector per configuration only. The gauge ensemble consisted of50 configurations. Again, a precise determination of 〈ψ̄ψ〉 could be obtained.These promising first results were obtained at zero temperature. We are interested in applying asimilar approach at T ∼ Tc, the deconfinement temperature, where the “gauge noise” is expected to bemuch greater. And we want, in addition to ψ̄ψ =∑i λ−1i , higher negative moments of the eigenvalues,i.e.∑i λ−mi ,m = 1, 2, 3, 4, ...3 Application I: Binder cumulant for the finite-temperature transitionA suitable order parameter for the finite-T transition is the quark condensate ψ̄ψ. To probe the orderof the phase transition, it is customary to measure its Binder cumulant (or equivalently its kurtosis)B4(ψ̄ψ) ≡〈(δψ̄ψ)4〉〈(δψ̄ψ)2〉2, δψ̄ψ = ψ̄ψ − 〈ψ̄ψ〉 (2)Thus, the first 4 moments (ψ̄ψ)k, k = 1, .., 4 must be measured, or estimated in an unbiased way. Thelatter approach is normally chosen, using noise vectors. One must then decide how many noise vectorsper configuration to use: using too few will increase the “stochastic noise” on B4 and the associatedstatistical error; using too many is wasteful of computing resources.Here, we chose the most economical approach: 4 noise vectors η1, .., η4, forming unbiased estimators 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0  0.2  0.4  0.6  0.8  1  1.2  1.4Estimator pbpExact pbp 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0  0.2  0.4  0.6  0.8  1  1.2  1.4Estimator pbp2Exact pbp2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0  0.2  0.4  0.6  0.8  1  1.2  1.4Estimator pbp3Exact pbp3 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0  0.2  0.4  0.6  0.8  1  1.2  1.4Estimator pbp4Exact pbp4 0 100 200 300 400 500 600 0.4  0.6  0.8  1  1.2Distribution of psibarpsiExactNoiseFigure 1. Comparison of the first 4 powers of ψ̄ψ with their stochastic estimators eq.(3), on an ensemble of 63×4SU(3),N f = 4 configurations with staggered fermions, near the finite-temperature deconfinement transition. Thebottom figure compares the resulting distribution of the exactly measured ψ̄ψ with that of its estimator: thetwo-bump distribution signals the vicinity of the finite-T transition; the "stochastic noise" barely increases thevariance of the distribution.asψ̄ψ ←−14(η†1D/−1η1 + η†2D/−1η2 + . . .)(ψ̄ψ)2 ←−16(η†1D/−1η1 · η†2D/−1η2 + . . .)(ψ̄ψ)3 ←−14(η†1D/−1η1 · η†2D/−1η2 · η†3D/−1η3 + . . .)(ψ̄ψ)4 ←−(η†1D/−1η1 · η†2D/−1η2 + η†3D/−1η3 · η†4D/−1η4)(3)In Fig. 1, we compare the estimators above with the exact values obtained by diagonalization of D/ ,for an ensemble of 63 × 4 SU(3) configurations, where the parameters (N f = 4, am = 0.07) have beenchosen to be in the vicinity of a sharp temperature-driven crossover. Indeed, the exact distributionof ψ̄ψ in the last panel shows a clear two-peak distribution. What is striking is that the distributionof the stochastic estimator of ψ̄ψ, in the same figure, is hardly broader, even though the number ofnoise vectors has been kept to its minimum value 4. Thus, in this example, the “stochastic noise” isnegligible compared to the “gauge noise”. Whether or not this persists on larger lattices needs to bechecked.0.2 0.4 0.6 0.8 1.0Chebychevapprox.-1.0 -0.5 0.5 1.0200400600800μ(-1,x)Figure 2. Chebyshev bandpass filter (left) and eigenvalue count µ(−1, x) as a function of x, for the rescaled,staggered operator (iD/ ) on a 44 lattice.The novel method outlined in Sec. 2 can give us nearly exact results at a potentially reduced cost,compared with exact diagonalization. A Chebyshev bandpass filter is applied to (iD/ ), so that thenumber µ of eigenvalues in [a, b] can be approximated asµ(a, b) =1nn∑j=1pmax∑p=0gp(a, b) 〈〈η†jTp(iD/ ) η j〉〉 (4)where Tp is the Chebyshev polynomial of degree p, gp are the coefficients of the filter approximation,and 〈〈· · · 〉〉 represents the averaging over the gauge field. Fig. 2 illustrates the bandpass filter (left) andthe eigenvalue count µ(−1, x) versus x (right) 1, for one 44 configuration (D/ has been rescaled so thatits spectrum is contained in [−1,+1]). µ(−1, x) is a piecewise constant function, increasing by one asx passes each eigenvalue. We approximate the eigenvalues as solutions of µ(−1, x) being half-integer.Since all eigenvalues are obtained in this way, we can form arbitrary moments Tr(D/ −k), k = 1, 2, . . ..This is illustrated in Table 1 below. Note the remarkable accuracy, even for high moments.k exact approx. approx./exact2 −453.767 −451.148 0.9942294 3871.66 3844.25 0.9929210 −8.136 · 107 −8.091 · 107 0.994547Table 1. Tr(D/ −k)using 4 noise vectors and a degree-16k polynomial on a 44 lattice.4 Application II: Taylor expansion of free energyQCD at non-zero quark chemical potential µ represents a long-standing challenge for lattice simu-lations: the Dirac determinant becomes complex when µ , 0, which prevents its interpretation as asampling probability. Several strategies have been proposed to circumvent this “sign” (or phase) prob-lem. One of the first is that of Taylor-expanding the QCD free energy, or the pressure P, in powers of(µ/T ) around µ = 0:P(µ,T ) = P(µ = 0,T ) +∑k=1ck(T )(µT)k. (5)Charge conjugation symmetry implies that odd-order coefficients ck vanish. However, products of oddderivatives do contribute in the presence of multiple quark flavors, each with its chemical potential.Now, ck is simply the kth derivative of the free energy evaluated at µ = 0, so it can be measured in con-ventional, sign-problem free simulations. The µ-dependence appears only in the Dirac determinant,so that the generic expression for ck is the trace of a polynomial expression in D/ −1 and in derivatives∂mD/ /∂µm. The number of terms grows exponentially with the order k, with positive or negative signs.This makes the calculation of high-order coefficients ck extremely challenging: no calculation hasever yielded information at 10th order. Improvements in efficiency are needed.One such improvement has been proposed recently by Gavai and Sharma [5, 6]. The suggestionis to change the way the chemical potential is introduced on the lattice, and adopt a naive, linearprescription rather than the exponential prescription of Hasenfratz and Karsch [7]. This is rathercounter-intuitive, since it re-introduces a UV divergence in the pressure. Gavai and Sharma argue thatthis divergence can be removed by a numerical subtraction. The gain is that derivatives ∂mD/ /∂µmnow vanish for m > 1. This greatly simplifies the structure of the ck’s, which become the traces ofhomogeneous polynomials of degree k in A−1, with A = D/(∂D/∂µ)−1(note that A is neither hermitian,nor anti-hermitian – see Fig. 3).At this stage, our approach may become advantageous. Suppose we can determine the [complex]spectrum of A. Then, ck will be obtained as a sum of products of factors Tr (A−m), which can all bedetermined in one go!1Note the similarity with Lüscher’s “mode number” [4].-4 -2 2 4-4-2244x4³-4 -2 2 4-4-2244x6³-4 -2 2 4-4-2244x8³Figure 3. Complex spectrum of(D/(∂D/∂µ)−1), for lattices of increasing spatial size near the finite-T transition.Here, we take first steps in that promising direction. First, we generalize the notion of eigenvaluecount to the complex plane. The number of eigenvalues µ(Γ) inside a complex contour Γ can beobtained by a Cauchy integral [8]:µ(Γ) =12πi∮Γdz Tr[(z1 − A)−1]=12πi∮Γdz∑j1z − λ j(6)For Γ, we choose a circle of center c and radius r, and discretize the integral using the trapezoidal rule:µ(r, c) ∼1NN−1∑k=0r e2πiN k Tr[(zk1 − A)−1]with zk = c + r e2πiN k (7)where the trace is estimated using noise vectors. The so-obtained value of µ(r, c) will be nearly aninteger equal to the number of eigenvalues inside the circle. One can now discard the circles whichcontain zero eigenvalue, and divide by 2 the radius of non-empty circles. This procedure is repeatedrecursively, until each circle contains one and only one eigenvalue, and has a radius smaller than someprescribed tolerance ε. In this fashion, the complete spectrum of A is determined. Arbitrary momentsTrA−m can be evaluated.The crucial point is that the estimation of the trace in eq.(7) proceeds by solving linear systems(zk1 − A) x = η, using a multi-shift linear solver. Then the cost of solving many such systems isnot much greater than that of solving just one, because these systems differ only by the shift zk, andthe Krylov space built by the linear solver is shift-invariant. In addition, each iteration of the solverinvolves matrix-vector multiplication by A = D/(∂D/∂µ)−1. At first sight, the inversion of ∂D/∂µlooksdaunting. But this operator is completely local in space, so that the corresponding matrix is made ofsmall diagonal blocks, each of which is easy to invert.We illustrate our approach in Figs. 3 and 4. Fig. 3 shows the complex spectrum of A, for 3configurations obtained with the same action, near the finite-T transition, on lattices of increasingspatial sizes.Fig. 4 shows 3 successive stages of our recursive procedure: at each stage, circles empty of eigen-values are discarded, and non-empty circles are replaced by 4 circles with half the radius. The finaloutput is an estimate of all eigenvalues of A, from which successive moments Tr(D/ ∂D/∂µ−1)−mcan beestimated as well. Illustrative results are provided in Table 2 below.-4 -2 2 4-4-224-4 -2 2 4-4-224-4 -2 2 4-4-224Figure 4. Circular contours are recursively refined until each contour encloses one eigenvalue (or a few).m exact approx. ratio1 4.10306i 0.0352058 + 4.10474i 1.000412 142.497 142.607 − 0.00748385i 1.000773 10.8641i 0.04117 + 10.8765i 1.001144 −14.4443 −14.3171 − 0.0751784i 0.99119510 −133.716 −135.113 − 3.93177i 1.0104520 −79320.5 −81511.3 + 979.464i 1.02762Table 2. Tr (A−m) for 100 quadrature points and 100 noise vectors for ε = 10−3, on a 44 gauge configuration.5 OutlookWe have described a novel approach to obtain estimates of moments Tr (A−m) of a sparse matrix A,in the Hermitian and in the non-Hermitian case. Two applications were tested: measuring the Bindercumulant of the quark condensate (A = iD/ ), and the coefficients of the Taylor expansion of thepressure in µ/T (A = D/(∂D/∂µ)−1). In both cases, the regime of interest is near the finite-temperaturedeconfinement transition. In this regime, fluctuations of the gauge field are large, so that the “gaugenoise” may well exceed the “stochastic noise” caused by the limited number of noise vectors used toestimate Tr (A−m). Indeed, this is what happens in the first application: it is then sufficient to use 4noise vectors per configuration, at least for the small lattice sizes considered in this test.In the second application, we have generalized the eigenvalue count to the complex plane. Theviability of our approach rests on a massive usage of non-Hermitian multi-shift linear solvers. Greatsimplifications in the contour integrals are possible and are under study. However, while we obtainnearly exact results for the first 20 moments Tr (A−m), we have not measured the “gauge noise”,which presumably increases quickly with m. This may represent the ultimate, irreducible origin of thecomputing cost of the Taylor expansion of the pressure, as foreseen in [9].6 AcknowledgementsBJ was supported by the Schweizerischer Nationalfonds (SNF) under grant 200020-162515. PdFthanks the Kavli Institute for Theoretical Physics for hospitality. This research was supported in partby the National Science Foundation under Grant No. NSF PHY-1125915.⇒ ⇒References[1] E. Di Napoli, E. Polizzi, Y. Saad, Numerical Linear Algebra with Applications 23, 674 (2016)[2] Z. Fodor, K. Holland, J. Kuti, S. Mondal, D. Nogradi, C.H. Wong, PoS LATTICE2015, 310(2016), 1605.08091[3] G. Cossu, H. Fukaya, S. Hashimoto, T. Kaneko, J.I. Noaki, PTEP 2016, 093B06 (2016),1607.01099[4] L. Giusti, M. Lüscher, JHEP 03, 013 (2009), 0812.3638[5] R.V. Gavai, S. Sharma, Phys. Lett. B749, 8 (2015), 1406.0474[6] R.V. Gavai, S. Sharma, PoS CPOD2014, 015 (2015), 1503.03575[7] P. Hasenfratz, F. Karsch, Physics Letters B 125, 308 (1983)[8] Y. Futamura, H. Tadano, T. Sakurai, JSIAM Letters 2, 127 (2010)[9] P. de Forcrand, PoS LAT2009, 010 (2009), 1005.0539

Several important observables, like the quark condensate and the Taylor coefficients of the expansion of the QCD pressure with respect to the chemical potential, are based on the trace of the inverse Dirac operator and of its powers. Such traces are traditionally estimated with “noise vectors” sandwiching the operator. We explore alternative approaches based on polynomial approximations of the inverse Dirac operator.Several important observables, like the quark condensate and the Taylor coefficients of the expansion of the QCD pressure with respect to the chemical potential, are based on the trace of the inverse Dirac operator and of its powers. Such traces are traditionally estimated with "noise vectors" sandwiching the operator. We explore alternative approaches based on polynomial approximations of the inverse Dirac operator

Jäger, Benjamin

CERN Document Server

Several important observables, like the quark condensate and the Taylor coefficients of the expansion of the QCD pressure with respect to the chemical potential, are based on the trace of the inverse Dirac operator and of its powers. Such traces are traditionally estimated with “noise vectors” sandwiching the operator. We explore alternative approaches based on polynomial approximations of the inverse Dirac operator.ISSN:2100-014XISSN:2101-627

Repository for Publications and Research Data

Several important observables, like the quark condensate and the Taylor coefficients of the expansion of the QCD pressure with respect to the chemical potential, are based on the trace of the inverse Dirac operator and of its powers. Such traces are traditionally estimated with “noise vectors” sandwiching the operator. We explore alternative approaches based on polynomial approximations of the inverse Dirac operator.</jats:p

Crossref

       University of Southern DenmarkAlternatives to the stochastic "noise vector" approachForcrand, Philippe de; Jaeger, BenjaminPublished in:EPJ Web of ConferencesDOI:10.1051/epjconf/201817514022Publication date:2018Document version:Submitted manuscriptCitation for pulished version (APA):Forcrand, P. D., & Jaeger, B. (2018). Alternatives to the stochastic "noise vector" approach. EPJ Web ofConferences, 175, Article 14022. https://doi.org/10.1051/epjconf/201817514022Go to publication entry in University of Southern Denmark's Research PortalTerms of useThis work is brought to you by the University of Southern Denmark.Unless otherwise specified it has been shared according to the terms for self-archiving.If no other license is stated, these terms apply:            • You may download this work for personal use only.            • You may not further distribute the material or use it for any profit-making activity or commercial gain            • You may freely distribute the URL identifying this open access versionIf you believe that this document breaches copyright please contact us providing details and we will investigate your claim.Please direct all enquiries to puresupport@bib.sdu.dkDownload date: 25. Oct. 2025CERN-TH-2017-210CP3-Origins-2017-44 DNRF90Alternatives to the stochastic “noise vector” approachPhilippe de Forcrand1,2,? and Benjamin Jäger1,31Institute for Theoretical Physics, ETH, CH-8093 Zürich, Switzerland2CERN, Theory Division, CH-1211 Geneva, Switzerland3CP3-Origins & Danish Institute for Advanced Study, Department of Mathematics and Computer Science,University of Southern Denmark, 5230 Odense M, DenmarkAbstract. Several important observables, like the quark condensate and the Taylor coef-ficients of the expansion of the QCD pressure with respect to the chemical potential, arebased on the trace of the inverse Dirac operator and of its powers. Such traces are tradi-tionally estimated with "noise vectors" sandwiching the operator. We explore alternativeapproaches based on polynomial approximations of the inverse Dirac operator.1 IntroductionMeasuring observables on an ensemble of gauge field configurations can require computing resourcessimilar to those required to generate the configurations. However, comparatively little effort hasbeen devoted to optimizing the computing strategy for such measurements. An important categoryof fermionic observables can be expressed as the trace of some function of the Dirac operator, orsome other related operator. For example, the quark condensate ψ̄ψ is obtained from the trace ofthe inverse of the Dirac operator D/ . Generically, the trace for a given gauge field configuration isestimated by using “noise vectors” ηk, k = 1, .., n, viaTrA ≈1nn∑k=1η†k Aηk (1)with 〈ηi†k ηjk〉 = δi j. Here, i and j label components of a noise vector η, and 〈..〉 means averaging overnoise vectors (for a given gauge field). Clearly, there are two sources of noise here: the “stochasticnoise” coming from the limited number n of noise vectors used on a given gauge field configuration;and the “gauge noise” coming from the fluctuations caused by the variation of the gauge field fromconfiguration to configuration. The balance between these two types of noise depends on the observ-able, and on the parameters of the Monte Carlo simulation: in particular, fluctuations coming fromthe gauge field will severely increase in the vicinity of a phase transition.Here, we consider two situations of relevance at finite temperature:(i) the measurement of the chiral condensate and its first 4 moments near the finite-temperature de-confinement transition. The moments are necessary to construct the Binder cumulant (or kurtosis),whose value at the phase transition is known by universality;?e-mail: forcrand@phys.ethz.charXiv:1710.07305v1  [hep-lat]  19 Oct 2017(ii) the measurement of the Taylor coefficients of the expansion of the pressure in µ/T .In the latter situation, it is not unusual to take 1000 noise vectors or more per gauge configuration.Can one be more efficient?2 Related recent progressApplied mathematicians have proposed a novel approach to determine the spectrum of a Hermitianmatrix A [1]. Instead of obtaining the eigenvalues directly by a Krylov-based method, the idea is todetermine the “eigenvalue count”, namely the number µ(a, b) of eigenvalues in some interval [a, b].Upon varying a and b, the positions of the jumps in the eigenvalue count give the desired eigenvalues.Thus, one needs a sliding bandpass filter f (x, a, b) which is = 1 if x ∈ [a, b], 0 otherwise. In practice:- The function f is approximated by Chebyshev polynomials and applied to the matrix A. Then,the “eigenvalue count” µ(a, b) is Tr f (A, a, b).- The trace is approximated by using a small number of noise vectors following eq.(1).One great advantage is that, by varying the coefficients of the polynomial, a and b can be varied atwill without any additional matrix-vector multiplication, and thereby the complete spectrum of A canbe determined at low cost.Two lattice groups have tested this approach so far:• In [2], a degree-8000 Chebyshev approximation was constructed and applied to D/ †D/ , using 20noise vectors for the trace. This required 80000 Dirac-matrix-vector multiplications. These are verymoderate resources, considering that the lattice size was 643 × 96, and that the whole spectrum couldbe extracted. After averaging over 10 gauge configurations, a precise determination of 〈ψ̄ψ〉 could beobtained.• In [3], an overlap Dirac operator was used, on the same lattice size, with a Chebyshev approximationof the same degree 8000, and 1 noise vector per configuration only. The gauge ensemble consisted of50 configurations. Again, a precise determination of 〈ψ̄ψ〉 could be obtained.These promising first results were obtained at zero temperature. We are interested in applying asimilar approach at T ∼ Tc, the deconfinement temperature, where the “gauge noise” is expected to bemuch greater. And we want, in addition to ψ̄ψ =∑i λ−1i , higher negative moments of the eigenvalues,i.e.∑i λ−mi ,m = 1, 2, 3, 4, ...3 Application I: Binder cumulant for the finite-temperature transitionA suitable order parameter for the finite-T transition is the quark condensate ψ̄ψ. To probe the orderof the phase transition, it is customary to measure its Binder cumulant (or equivalently its kurtosis)B4(ψ̄ψ) ≡〈(δψ̄ψ)4〉〈(δψ̄ψ)2〉2, δψ̄ψ = ψ̄ψ − 〈ψ̄ψ〉 (2)Thus, the first 4 moments (ψ̄ψ)k, k = 1, .., 4 must be measured, or estimated in an unbiased way. Thelatter approach is normally chosen, using noise vectors. One must then decide how many noise vectorsper configuration to use: using too few will increase the “stochastic noise” on B4 and the associatedstatistical error; using too many is wasteful of computing resources.Here, we chose the most economical approach: 4 noise vectors η1, .., η4, forming unbiased estimators 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0  0.2  0.4  0.6  0.8  1  1.2  1.4Estimator pbpExact pbp 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0  0.2  0.4  0.6  0.8  1  1.2  1.4Estimator pbp2Exact pbp2 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0  0.2  0.4  0.6  0.8  1  1.2  1.4Estimator pbp3Exact pbp3 0 0.2 0.4 0.6 0.8 1 1.2 1.4 0  0.2  0.4  0.6  0.8  1  1.2  1.4Estimator pbp4Exact pbp4 0 100 200 300 400 500 600 0.4  0.6  0.8  1  1.2Distribution of psibarpsiExactNoiseFigure 1. Comparison of the first 4 powers of ψ̄ψ with their stochastic estimators eq.(3), on an ensemble of 63×4SU(3),N f = 4 configurations with staggered fermions, near the finite-temperature deconfinement transition. Thebottom figure compares the resulting distribution of the exactly measured ψ̄ψ with that of its estimator: thetwo-bump distribution signals the vicinity of the finite-T transition; the "stochastic noise" barely increases thevariance of the distribution.asψ̄ψ ←−14(η†1D/−1η1 + η†2D/−1η2 + . . .)(ψ̄ψ)2 ←−16(η†1D/−1η1 · η†2D/−1η2 + . . .)(ψ̄ψ)3 ←−14(η†1D/−1η1 · η†2D/−1η2 · η†3D/−1η3 + . . .)(ψ̄ψ)4 ←−(η†1D/−1η1 · η†2D/−1η2 + η†3D/−1η3 · η†4D/−1η4)(3)In Fig. 1, we compare the estimators above with the exact values obtained by diagonalization of D/ ,for an ensemble of 63 × 4 SU(3) configurations, where the parameters (N f = 4, am = 0.07) have beenchosen to be in the vicinity of a sharp temperature-driven crossover. Indeed, the exact distributionof ψ̄ψ in the last panel shows a clear two-peak distribution. What is striking is that the distributionof the stochastic estimator of ψ̄ψ, in the same figure, is hardly broader, even though the number ofnoise vectors has been kept to its minimum value 4. Thus, in this example, the “stochastic noise” isnegligible compared to the “gauge noise”. Whether or not this persists on larger lattices needs to bechecked.0.2 0.4 0.6 0.8 1.0Chebychevapprox.-1.0 -0.5 0.5 1.0200400600800μ(-1,x)Figure 2. Chebyshev bandpass filter (left) and eigenvalue count µ(−1, x) as a function of x, for the rescaled,staggered operator (iD/ ) on a 44 lattice.The novel method outlined in Sec. 2 can give us nearly exact results at a potentially reduced cost,compared with exact diagonalization. A Chebyshev bandpass filter is applied to (iD/ ), so that thenumber µ of eigenvalues in [a, b] can be approximated asµ(a, b) =1nn∑j=1pmax∑p=0gp(a, b) 〈〈η†jTp(iD/ ) η j〉〉 (4)where Tp is the Chebyshev polynomial of degree p, gp are the coefficients of the filter approximation,and 〈〈· · · 〉〉 represents the averaging over the gauge field. Fig. 2 illustrates the bandpass filter (left) andthe eigenvalue count µ(−1, x) versus x (right) 1, for one 44 configuration (D/ has been rescaled so thatits spectrum is contained in [−1,+1]). µ(−1, x) is a piecewise constant function, increasing by one asx passes each eigenvalue. We approximate the eigenvalues as solutions of µ(−1, x) being half-integer.Since all eigenvalues are obtained in this way, we can form arbitrary moments Tr(D/ −k), k = 1, 2, . . ..This is illustrated in Table 1 below. Note the remarkable accuracy, even for high moments.k exact approx. approx./exact2 −453.767 −451.148 0.9942294 3871.66 3844.25 0.9929210 −8.136 · 107 −8.091 · 107 0.994547Table 1. Tr(D/ −k)using 4 noise vectors and a degree-16k polynomial on a 44 lattice.4 Application II: Taylor expansion of free energyQCD at non-zero quark chemical potential µ represents a long-standing challenge for lattice simu-lations: the Dirac determinant becomes complex when µ , 0, which prevents its interpretation as asampling probability. Several strategies have been proposed to circumvent this “sign” (or phase) prob-lem. One of the first is that of Taylor-expanding the QCD free energy, or the pressure P, in powers of(µ/T ) around µ = 0:P(µ,T ) = P(µ = 0,T ) +∑k=1ck(T )(µT)k. (5)Charge conjugation symmetry implies that odd-order coefficients ck vanish. However, products of oddderivatives do contribute in the presence of multiple quark flavors, each with its chemical potential.Now, ck is simply the kth derivative of the free energy evaluated at µ = 0, so it can be measured in con-ventional, sign-problem free simulations. The µ-dependence appears only in the Dirac determinant,so that the generic expression for ck is the trace of a polynomial expression in D/ −1 and in derivatives∂mD/ /∂µm. The number of terms grows exponentially with the order k, with positive or negative signs.This makes the calculation of high-order coefficients ck extremely challenging: no calculation hasever yielded information at 10th order. Improvements in efficiency are needed.One such improvement has been proposed recently by Gavai and Sharma [5, 6]. The suggestionis to change the way the chemical potential is introduced on the lattice, and adopt a naive, linearprescription rather than the exponential prescription of Hasenfratz and Karsch [7]. This is rathercounter-intuitive, since it re-introduces a UV divergence in the pressure. Gavai and Sharma argue thatthis divergence can be removed by a numerical subtraction. The gain is that derivatives ∂mD/ /∂µmnow vanish for m > 1. This greatly simplifies the structure of the ck’s, which become the traces ofhomogeneous polynomials of degree k in A−1, with A = D/(∂D/∂µ)−1(note that A is neither hermitian,nor anti-hermitian – see Fig. 3).At this stage, our approach may become advantageous. Suppose we can determine the [complex]spectrum of A. Then, ck will be obtained as a sum of products of factors Tr (A−m), which can all bedetermined in one go!1Note the similarity with Lüscher’s “mode number” [4].-4 -2 2 4-4-2244x4³-4 -2 2 4-4-2244x6³-4 -2 2 4-4-2244x8³Figure 3. Complex spectrum of(D/(∂D/∂µ)−1), for lattices of increasing spatial size near the finite-T transition.Here, we take first steps in that promising direction. First, we generalize the notion of eigenvaluecount to the complex plane. The number of eigenvalues µ(Γ) inside a complex contour Γ can beobtained by a Cauchy integral [8]:µ(Γ) =12πi∮Γdz Tr[(z1 − A)−1]=12πi∮Γdz∑j1z − λ j(6)For Γ, we choose a circle of center c and radius r, and discretize the integral using the trapezoidal rule:µ(r, c) ∼1NN−1∑k=0r e2πiN k Tr[(zk1 − A)−1]with zk = c + r e2πiN k (7)where the trace is estimated using noise vectors. The so-obtained value of µ(r, c) will be nearly aninteger equal to the number of eigenvalues inside the circle. One can now discard the circles whichcontain zero eigenvalue, and divide by 2 the radius of non-empty circles. This procedure is repeatedrecursively, until each circle contains one and only one eigenvalue, and has a radius smaller than someprescribed tolerance ε. In this fashion, the complete spectrum of A is determined. Arbitrary momentsTrA−m can be evaluated.The crucial point is that the estimation of the trace in eq.(7) proceeds by solving linear systems(zk1 − A) x = η, using a multi-shift linear solver. Then the cost of solving many such systems isnot much greater than that of solving just one, because these systems differ only by the shift zk, andthe Krylov space built by the linear solver is shift-invariant. In addition, each iteration of the solverinvolves matrix-vector multiplication by A = D/(∂D/∂µ)−1. At first sight, the inversion of ∂D/∂µlooksdaunting. But this operator is completely local in space, so that the corresponding matrix is made ofsmall diagonal blocks, each of which is easy to invert.We illustrate our approach in Figs. 3 and 4. Fig. 3 shows the complex spectrum of A, for 3configurations obtained with the same action, near the finite-T transition, on lattices of increasingspatial sizes.Fig. 4 shows 3 successive stages of our recursive procedure: at each stage, circles empty of eigen-values are discarded, and non-empty circles are replaced by 4 circles with half the radius. The finaloutput is an estimate of all eigenvalues of A, from which successive moments Tr(D/ ∂D/∂µ−1)−mcan beestimated as well. Illustrative results are provided in Table 2 below.-4 -2 2 4-4-224-4 -2 2 4-4-224-4 -2 2 4-4-224Figure 4. Circular contours are recursively refined until each contour encloses one eigenvalue (or a few).m exact approx. ratio1 4.10306i 0.0352058 + 4.10474i 1.000412 142.497 142.607 − 0.00748385i 1.000773 10.8641i 0.04117 + 10.8765i 1.001144 −14.4443 −14.3171 − 0.0751784i 0.99119510 −133.716 −135.113 − 3.93177i 1.0104520 −79320.5 −81511.3 + 979.464i 1.02762Table 2. Tr (A−m) for 100 quadrature points and 100 noise vectors for ε = 10−3, on a 44 gauge configuration.5 OutlookWe have described a novel approach to obtain estimates of moments Tr (A−m) of a sparse matrix A,in the Hermitian and in the non-Hermitian case. Two applications were tested: measuring the Bindercumulant of the quark condensate (A = iD/ ), and the coefficients of the Taylor expansion of thepressure in µ/T (A = D/(∂D/∂µ)−1). In both cases, the regime of interest is near the finite-temperaturedeconfinement transition. In this regime, fluctuations of the gauge field are large, so that the “gaugenoise” may well exceed the “stochastic noise” caused by the limited number of noise vectors used toestimate Tr (A−m). Indeed, this is what happens in the first application: it is then sufficient to use 4noise vectors per configuration, at least for the small lattice sizes considered in this test.In the second application, we have generalized the eigenvalue count to the complex plane. Theviability of our approach rests on a massive usage of non-Hermitian multi-shift linear solvers. Greatsimplifications in the contour integrals are possible and are under study. However, while we obtainnearly exact results for the first 20 moments Tr (A−m), we have not measured the “gauge noise”,which presumably increases quickly with m. This may represent the ultimate, irreducible origin of thecomputing cost of the Taylor expansion of the pressure, as foreseen in [9].6 AcknowledgementsBJ was supported by the Schweizerischer Nationalfonds (SNF) under grant 200020-162515. PdFthanks the Kavli Institute for Theoretical Physics for hospitality. This research was supported in partby the National Science Foundation under Grant No. NSF PHY-1125915.⇒ ⇒References[1] E. Di Napoli, E. Polizzi, Y. Saad, Numerical Linear Algebra with Applications 23, 674 (2016)[2] Z. Fodor, K. Holland, J. Kuti, S. Mondal, D. Nogradi, C.H. Wong, PoS LATTICE2015, 310(2016), 1605.08091[3] G. Cossu, H. Fukaya, S. Hashimoto, T. Kaneko, J.I. Noaki, PTEP 2016, 093B06 (2016),1607.01099[4] L. Giusti, M. Lüscher, JHEP 03, 013 (2009), 0812.3638[5] R.V. Gavai, S. Sharma, Phys. Lett. B749, 8 (2015), 1406.0474[6] R.V. Gavai, S. Sharma, PoS CPOD2014, 015 (2015), 1503.03575[7] P. Hasenfratz, F. Karsch, Physics Letters B 125, 308 (1983)[8] Y. Futamura, H. Tadano, T. Sakurai, JSIAM Letters 2, 127 (2010)[9] P. de Forcrand, PoS LAT2009, 010 (2009), 1005.0539

https://findresearcher.sdu.dk/ws/files/129361665/1710.07305v1.pdf

Alternatives to the stochastic "noise vector" approach

Abstract

Similar works

Full text

Available Versions

EDP Sciences OAI-PMH repository (1.2.0)

Syddansk Universitets Forskerportal

Syddansk Universitets Forskerportal

CERN Document Server

Repository for Publications and Research Data

Crossref

Syddansk Universitets Forskerportal