Distance preconditioning for lattice Dirac operators

We propose a preconditioning of the Dirac operator based on the factorisation of a predefined function related to the decay of the propagator with the distance. We show that it can improve the accuracy of correlators involving heavy quarks at large distances and accelerate the computation of light quark propagators.


I. INTRODUCTION
A key ingredient of lattice QCD simulations is the inversion of the Dirac operator which enters the generation of unquenched gauge field configurations and the computation of hadronic observables.One needs to solve numerically a linear system of the form where D[U ] is the chosen discretization of the massless interacting Dirac operator, M is the quark mass in lattice units, η(x) is a source vector that is different from zero on a single time-slice (that without any loss we shall assume to be at x 0 = 0).The solution S(y) is obtained by iterative numerical algorithms, solvers, devised to invert so-called sparse matrices, like the matrices that result from the discretization of differential equations by finite differences methods.In this letter we shall not discuss the details of any particular solver (see ref. [1] for a complete review and for an updated list of references).For any solver one checks if the condition is satisfyed.Here S n (y) is the tentative solution at iteration number n, the norm is any good norm in field space and r, the residue, is the global numerical accuracy requested for the solution.Typically r is a small number of the order of the arithmetic precision allowed by the computer architecture.Depending on the values of the quark mass the solution of eq. ( 1) poses different numerical problems.For light quarks the matrix (D[U ] + M ) x,y is badly conditioned and its numerical inversion requires a big number of iterations.At the other extreme, the number of iterations required for heavy quark masses is small but there may be problems with the numerical accuracy resulting for the time-slices far away from the source (y 0 0).Indeed eq. ( 2) is a global condition while for heavy quark propagators the time-slices far away from the source are exponentially suppresed by a factor of the order of exp(−M y 0 ) and give a negligible contribution to the norm on the left side of eq. ( 2).When this problem arises one cannot trust numerical results at large times and it becomes impossible to extract physical informations by fitting the leading exponentials contributing to correlation functions.
In order to alleviate both difficulties, we propose a preconditioning of the Dirac operator that factorises from the propagator a function aiming to modify its leading decay with the distance 1 .The simplest choice is to factorize a function α(y 0 ), to solve numerically the preconditioned equation, and to restore the original propagator by multiplying each time slice for 1/α(y 0 ).α(y 0 ) is defined such that all the different time-slices give comparable contributions to the calculation of the residue in the preconditioned case.Our preconditioning is inspired to what is usually done in deriving the Eichten and Hill [2] lattice HQET action but of course does not introduce any approximation.Indeed, the choice above is suited for heavy quark propagators, while for light quark masses we will introduce a generalisation of the factorised function.

II. PRECONDITIONING HEAVY QUARK PROPAGATORS
We work with the O(a)-improved Wilson lattice Dirac operator but the numerical problems that we address arise also with alternative discretisations of the continuum action and the proposed solution can as well be easily implemented in those cases.We have tested our preconditioning scheme both for heavy and for light quark masses and in the free and in the interacting case.We start with the results for the heavy quarks.We first want to pick up a case where the problem arises.As an example, we have calculated the correlation function by solving eq. ( 1) for a heavy quark propagator of mass M 0.5 in the free theory for different choices of the residue.More precisely the red points in FIG. 1 have been obtained with a residue r = 10 −11 while the black points with a residue r = 10 −6 and the two black lines correspond to the squares of these two values of r.As is clearly visible from FIG. 1, and from FIG. 2 where we show the effective masses of the correlations shown in FIG. 1, the black points start to deviate from the red ones for y 0 18, i.e. when the correlator, which in this case is just the square module of the propagator, becomes smaller than the square of the "loose" residue r = 10 −6 .
If the time extent of the lattice is not too large the problem can be solved by brute force by lowering the residue and the results obtained in the preset case with r = 10 −11 can be considered as exact.If instead the time extent of the lattice is rather large the brute force approach cannot be considered because the required residues would be smaller than what is allowed on doubleprecision architectures, even in the case of moderately heavy quarks.In the case under consideration, by choosing a loose precision, i.e. a residue r = 10 −6 , we make the numerical problem evident and we show that also such an "extreme" situation can be recovered by using our proposal.Moreover, we notice that a residue r = 10 −6 is the smallest allowed on single-precision architectures that presently are considerably much faster than doubleprecision ones.
We now come to the proposed solution.We redefine the quark fields and the propagators as follows Once the previous expressions are inserted in eq. ( 1) we get the preconditioned system that we solve numerically in place of eq. ( 1).In order to write the preconditioned Dirac operator it is sufficient to modify the forward and backward lattice covariant derivatives in the time direction accoring to ∇ 0 S(y) = U 0 (y)S(y + 0) − S(y) → α(y 0 + 1) α(y 0 ) U 0 (y)S (y + 0) − S (y) Particular care has to be used at the boundaries of the lattice in order to respect the boundary conditions originally satisfied by the quark fields.If as in the case of FIG. 1 S(y) satisfies anti-periodic boundary conditions along the time direction, it follows from eq. ( 4) that The blue points in FIG. 1 correspond to the correlation function obtained by solving eq. ( 5) with the loose residue r = 10 −6 but after having factorized the function by setting α 0 = 0.4.As expected, the preconditioned correlator stays above the line of the loose precision residue and the "exact" result can be back recovered as follows In FIG. 2 the blue points correspond to the effective mass of the preconditioned correlator after the "restoration" of eq. ( 10) and fall exactly on top of the red ones in spite of the fact that they have been obtained with the same loose precision that affected the non preconditioned black points.
In FIG. 3 we show the same plot as in FIG. 2 but in the interacting theory.The gauge ensamble used correspond to the entry E5 in TABLE I.The size of the lattice is L 0 L 1 L 2 L 3 = 64 × 32 3 and the hopping parameter of the sea quarks is k sea = 0.13625 corresponding to a PCAC quark mass of about am P CAC sea 0.07.The data shown in FIG. 3 correspond to a pseudoscalar-pseudoscalar correlator, as in the free theory case, of two degenerate heavy quarks with hopping parameters k h = 0.125 corresponding to a PCAC quark mass of about am P CAC h 0.35.The unpreconditioned correlators decay approximately as fast as in the free theory case and from the difference of the black (unpreconditioned, r = 10 −6 ) and red (unpreconditioned, r = 10 −11 ) sets of data we see the same distortion of FIG. 2. The blue points have been obtained by solving eq. ( 5) after having factorized α(y 0 ) with α 0 = 0.4 and by restoring the results according to eq. ( 10).Also in the interacting theory the blue points are identical to red points though they have been obtained with the same loose residue r = 10 −6 used to obtain the black points.
We close this section by observing that our preconditioning technique may be particularly useful when working with the Schrödinger Functional [5,6] formulation of the theory.In this case, countrary to the case of periodic boundary conditions along the time direction, the correlators decay exponantially over the whole time extent of the lattice and one has to choose very small residues also in computing relatively light quark propagators.We have performed several succesful experiments with our preconditioning technique also in the Schrödinger Functional case by using α(x 0 ) = exp(−α 0 x 0 ).

III. PRECONDITIONING LIGHT QUARK PROPAGATORS
In this section we shall briefly discuss how the ideas developed and discussed in the previous section can be used to accelerate the numerical calculation of light quark propagators.We start our discussion by generalizing eq. ( 4) as follows S(y) = β(y 0 , y 1 , y 2 , y 3 ) S (y) In the following we shall consider the particular choice and the preconditioned lattice Dirac operator can be obtained as easily as before by changing all the covariant The values corresponding to α0 = 0.0 have been obtained without using our preconditioning technique.
derivatives according to and by changing accordingly the boundary conditions in all directions as done in eqs.(7) for the time direction.
An important difference of the present case with respect to the one discussed in the previous section is that the restoration of the true propagator must be performed before making the contractions needed to build correlation functions by using the first of eqs.(11).
Here the preconditioning is not to gain precision, but to accelerate the convergence of the inversion.Therefore, by applying eqs.(11), ( 12) and (13) to the calculation of a light quark propagator one aims to make the propagator to decay faster than the original unpreconditioned operator.By judiciously chosing the parameter α 0 it is possible to change the condition number of the preconditioned system without compromising the numerical accuracy of the solution, an operation that should be performed on double-precision computer architectures.
In TABLE I we quantify the gain in computational time that can be achieved by showing the number of iterations of the SAP+GCR solver required to solve the lattice Dirac equation for light quarks with and without our preconditioning.The SAP+GCR solver has been introduced and explained in details by the author in ref. [7].The gauge ensambles used to perform this test have been generated within the CLS agreement [8] with the parameters given in the table.In the case of the E-lattices the SAP+GCR solver has been ran on 128 processors of a cluster of PC's by dividing the global lattices into blocks of 4 4 points.In the case of the D-lattice the SAP+GCR solver has been ran on 32 processors of a cluster of PC's by dividing the global lattices into blocks of 6 4 points.The table shows that by increasing the value of the parameter α 0 the number of iterations goes down with a time gain that can easily reach the 30%.In the case under discussion, we checked that higher values of α 0 would induce a "heavy quark" like behavior and produce distorted results for the reasons discussed at lenghty in the previous section.
The method discussed in this letter can be generalised by adding some Dirac structure in the factorised function, an option presently under investigation.

FIG. 1 :FIG. 2 :
FIG.1:The red points correspond to the correlation fuction −CP P (y0) obtained by inverting the lattice Dirac equation in the free theory with M 0.5 and with a residue r = 10 −11 .The black points correspond to the same quantity but have been obtained with a residue r = 10 −6 .The blue points correspond to the correlation function −C P P (y0) obtained by solving the preconditioned lattice Dirac equation with M 0.5 and α0 = 0.4.The two black lines correspond to r 2 for the two values of the residue used in the calculations.We use logarithmic scale on the y-axis.

FIG. 3 :
FIG. 3: The red points correspond to the effective mass of the correlation fuction −CP P (y0) obtained by inverting the lattice Dirac equation in the interacting theory with am pcac h 0.35 and with a residue r = 10 −11 .The black points correspond to the same quantity but have been obtained with a residue r = 10 −6 .The blue points correspond to the effective mass of the restored correlation function −C P P (y0) obtained by solving the preconditioned lattice Dirac equation with am pcac h 0.35 and α0 = 0.4.