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.


Introduction
Measuring observables on an ensemble of gauge field configurations can require computing resources similar to those required to generate the configurations. However, comparatively little effort has been devoted to optimizing the computing strategy for such measurements. An important category of fermionic observables can be expressed as the trace of some function of the Dirac operator, or some other related operator. For example, the quark condensateψψ is obtained from the trace of the inverse of the Dirac operator D / . Generically, the trace for a given gauge field configuration is estimated by using "noise vectors" η k , k = 1, .., n, via with η i † k η j k = δ i j . Here, i and j label components of a noise vector η, and .. means averaging over noise vectors (for a given gauge field). Clearly, there are two sources of noise here: the "stochastic noise" 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 from configuration to configuration. The balance between these two types of noise depends on the observable, and on the parameters of the Monte Carlo simulation: in particular, fluctuations coming from the 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 deconfinement 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.ch (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?

Related recent progress
Applied mathematicians have proposed a novel approach to determine the spectrum of a Hermitian matrix A [1]. Instead of obtaining the eigenvalues directly by a Krylov-based method, the idea is to determine 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 at will without any additional matrix-vector multiplication, and thereby the complete spectrum of A can be 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 20 noise vectors for the trace. This required 80000 Dirac-matrix-vector multiplications. These are very moderate resources, considering that the lattice size was 64 3 × 96, and that the whole spectrum could be extracted. After averaging over 10 gauge configurations, a precise determination of ψ ψ could be obtained.
• In [3], an overlap Dirac operator was used, on the same lattice size, with a Chebyshev approximation of the same degree 8000, and 1 noise vector per configuration only. The gauge ensemble consisted of 50 configurations. Again, a precise determination of ψ ψ could be obtained.
These promising first results were obtained at zero temperature. We are interested in applying a similar approach at T ∼ T c , the deconfinement temperature, where the "gauge noise" is expected to be much greater. And we want, in addition toψψ = i λ −1 i , higher negative moments of the eigenvalues, i.e. i λ −m i , m = 1, 2, 3, 4, ...

Application I: Binder cumulant for the finite-temperature transition
A suitable order parameter for the finite-T transition is the quark condensateψψ. To probe the order of the phase transition, it is customary to measure its Binder cumulant (or equivalently its kurtosis) Thus, the first 4 moments (ψψ) k , k = 1, .., 4 must be measured, or estimated in an unbiased way. The latter approach is normally chosen, using noise vectors. One must then decide how many noise vectors per configuration to use: using too few will increase the "stochastic noise" on B 4 and the associated statistical error; using too many is wasteful of computing resources.
Here, we chose the most economical approach: 4 noise vectors η 1 , .., η 4 , forming unbiased estimators asψ In Fig. 1, we compare the estimators above with the exact values obtained by diagonalization of D / , for an ensemble of 6 3 × 4 SU(3) configurations, where the parameters (N f = 4, am = 0.07) have been chosen to be in the vicinity of a sharp temperature-driven crossover. Indeed, the exact distribution ofψψ in the last panel shows a clear two-peak distribution. What is striking is that the distribution of the stochastic estimator ofψψ, in the same figure, is hardly broader, even though the number of noise vectors has been kept to its minimum value 4. Thus, in this example, the "stochastic noise" is negligible compared to the "gauge noise". Whether or not this persists on larger lattices needs to be checked.  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 the number µ of eigenvalues in [a, b] can be approximated as where T p is the Chebyshev polynomial of degree p, g p are the coefficients of the filter approximation, and · · · represents the averaging over the gauge field. Fig. 2 illustrates the bandpass filter (left) and the eigenvalue count µ(−1, x) versus x (right) 1 , for one 4 4 configuration (D / has been rescaled so that its spectrum is contained in [−1, +1]). µ(−1, x) is a piecewise constant function, increasing by one as x 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

Application II: Taylor expansion of free energy
QCD at non-zero quark chemical potential µ represents a long-standing challenge for lattice simulations: the Dirac determinant becomes complex when µ 0, which prevents its interpretation as a sampling probability. Several strategies have been proposed to circumvent this "sign" (or phase) problem. One of the first is that of Taylor-expanding the QCD free energy, or the pressure P, in powers of (µ/T ) around µ = 0: Charge conjugation symmetry implies that odd-order coefficients c k vanish. However, products of odd derivatives do contribute in the presence of multiple quark flavors, each with its chemical potential. Now, c k is simply the k th derivative of the free energy evaluated at µ = 0, so it can be measured in conventional, sign-problem free simulations. The µ-dependence appears only in the Dirac determinant, so that the generic expression for c k is the trace of a polynomial expression in D / −1 and in derivatives ∂ m D / /∂µ m . The number of terms grows exponentially with the order k, with positive or negative signs. This makes the calculation of high-order coefficients c k extremely challenging: no calculation has ever yielded information at 10 th order. Improvements in efficiency are needed. One such improvement has been proposed recently by Gavai and Sharma [5,6]. The suggestion is to change the way the chemical potential is introduced on the lattice, and adopt a naive, linear prescription rather than the exponential prescription of Hasenfratz and Karsch [7]. This is rather counter-intuitive, since it re-introduces a UV divergence in the pressure. Gavai and Sharma argue that this divergence can be removed by a numerical subtraction. The gain is that derivatives ∂ m D / /∂µ m now vanish for m > 1. This greatly simplifies the structure of the c k 's, which become the traces of homogeneous 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, c k will be obtained as a sum of products of factors Tr (A −m ), which can all be determined in one go! Here, we take first steps in that promising direction. First, we generalize the notion of eigenvalue count to the complex plane. The number of eigenvalues µ(Γ) inside a complex contour Γ can be obtained by a Cauchy integral [8]: For Γ, we choose a circle of center c and radius r, and discretize the integral using the trapezoidal rule: where the trace is estimated using noise vectors. The so-obtained value of µ(r, c) will be nearly an integer equal to the number of eigenvalues inside the circle. One can now discard the circles which contain zero eigenvalue, and divide by 2 the radius of non-empty circles. This procedure is repeated recursively, until each circle contains one and only one eigenvalue, and has a radius smaller than some prescribed tolerance . In this fashion, the complete spectrum of A is determined. Arbitrary moments TrA −m can be evaluated. The crucial point is that the estimation of the trace in eq.(7) proceeds by solving linear systems (z k 1 − A) x = η, using a multi-shift linear solver. Then the cost of solving many such systems is not much greater than that of solving just one, because these systems differ only by the shift z k , and the Krylov space built by the linear solver is shift-invariant. In addition, each iteration of the solver involves matrix-vector multiplication by A = D / ∂D / ∂µ −1 . At first sight, the inversion of ∂D / ∂µ looks daunting. But this operator is completely local in space, so that the corresponding matrix is made of small 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 3 configurations obtained with the same action, near the finite-T transition, on lattices of increasing spatial sizes.  Table 2 below.   Table 2. Tr (A −m ) for 100 quadrature points and 100 noise vectors for = 10 −3 , on a 4 4 gauge configuration.

Outlook
We ). In both cases, the regime of interest is near the finite-temperature deconfinement transition. In this regime, fluctuations of the gauge field are large, so that the "gauge noise" may well exceed the "stochastic noise" caused by the limited number of noise vectors used to estimate Tr (A −m ). Indeed, this is what happens in the first application: it is then sufficient to use 4 noise 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. The viability of our approach rests on a massive usage of non-Hermitian multi-shift linear solvers. Great simplifications in the contour integrals are possible and are under study. However, while we obtain nearly 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 the computing cost of the Taylor expansion of the pressure, as foreseen in [9].