How to compute Green's Functions for entire Mass Trajectories within Krylov Solvers

TL;DR

This paper introduces a novel numerical method integrated within Krylov solvers to efficiently compute multiple Green's functions and their derivatives for matrices of the form A=D-m, enabling rapid evaluation across a range of mass parameters in lattice QCD.

Contribution

The paper presents a new general approach to compute many Green's functions simultaneously within one iteration, applicable to matrices of a specific structure, and demonstrates its effectiveness in lattice QCD applications.

Findings

Efficient computation of Green's functions for multiple masses using a single inversion.

Ability to derive derivatives and Taylor expansions of solutions with respect to parameters.

Application to lattice QCD shows significant computational savings.

Abstract

The availability of efficient Krylov subspace solvers play a vital role for the solution of a variety of numerical problems in computational science. Here we consider lattice field theory. We present a new general numerical method to compute many Green's functions for complex non-singular matrices within one iteration process. Our procedure applies to matrices of structure $A=D-m$, with $m$ proportional to the unit matrix, and can be integrated within any Krylov subspace solver. We can compute the derivatives $x^{(n)}$ of the solution vector $x$ with respect to the parameter $m$ and construct the Taylor expansion of $x$ around $m$. We demonstrate the advantages of our method using a minimal residual solver. Here the procedure requires $1$ intermediate vector for each Green's function to compute. As real life example, we determine a mass trajectory of the Wilson fermion matrix for lattice QCD. Here we find that we can obtain Green's functions at all masses $\geq m$ at the price of one inversion at mass $m$.

Equations

$$Ax=(D-m)x=r_0$$

$$x_{i+1}$$

$$r_{i+1}$$

$${\alpha }_i=\frac{(A{r}_i{)}^{†}{r}_i}{(A{r}_i{)}^{†}A{r}_i}.$$

$$A{x}^{(n)}-n{x}^{(n-1)}=0.$$

$$A{x}_i^{(1)}=x_i$$

$$A(x_{i+1}^{(1)}-x_i^{(1)})=(x_{i+1}-x_i)$$

$$x_{i+1}^{(1)}-x_i^{(1)}={\alpha }_i(x-x_i).$$

$$x_{i+1}^{(1)}-x_i^{(1)}={\alpha }_i\left[\sum _{j=0}^N{\alpha }_j{r}_j-\sum _{j=0}^{i-1}{\alpha }_j{r}_j\right]$$

$$x_{N+1}^{(1)}-x_0^{(1)}=\sum _{i=0}^N{\alpha }_i\sum _{j=i}^N{\alpha }_j{r}_j.$$

$$x_{N+1}^{(1)}=\sum _{i=0}^N{\alpha }_i\left[\sum _{j=0}^i{\alpha }_j\right]r_i.$$

$${\alpha }_i^{(1)}={\alpha }_i\sum _{j=0}^i{\alpha }_j$$

$$x_{i+1}^{(1)}=x_i^{(1)}+{\alpha }_i^{(1)}{r}_i.$$

$${\alpha }_i^{(n+1)}={\alpha }_i\sum _{j=0}^i{\alpha }_j^{(n)}.$$

$$x_{i+1}^{(n)}=x_i^{(n)}+n!{\alpha }_i^{(n)}{r}_i.$$

$$x^E=\sum _{n=0}^N\frac{(\Delta m{)}^n}{n!}{x}^{(n)}.$$

$$x_{i+1}^E=x_i^E+\sum _{n=0}^N(\Delta m{)}^n{\alpha }_i^{(n)}{r}_i:=x_i^E+f_i{\alpha }_i{r}_i,$$

$$f_i=1+\frac{1}{{\alpha }_i}\sum _{n=1}^N\Delta m^n{\alpha }_i^{(n)}.$$

$${\alpha }_i^{(n)}=\frac{{\alpha }_i}{{\alpha }_{i-1}}{\alpha }_{i-1}^{(n)}+{\alpha }_i{\alpha }_i^{(n-1)}.$$

$$f_i$$

$$f_0$$

$$M=m-D$$

$$D_{xy}=\sum _{\mu =1}^4(1-{\gamma }_{\mu })U_{\mu }(x){\delta }_{x,y-\mu }+(1+{\gamma }_{\mu })U_{\mu }^{†}(x-\mu ){\delta }_{x,y+\mu }$$

$$Mx=\Phi ,$$

$$\kappa =\frac{1}{m}=0.1575,0.1570,0.1565,0.1555,0.1530,0.1400,0.1000.$$

$$\frac{||r_i||}{||x_i||}<10^{-5}$$

$$\left(\begin{array}[]{cc}m&-D_{eo}\\ -D_{oe}&m\end{array}\right)\left(\begin{array}[]{c}x_{e}\\ x_{o}\end{array}\right)=\left(\begin{array}[]{c}\Phi_{e}\\ \Phi_{o}\end{array}\right).$$

$$M^{\prime}=\left(\begin{array}[]{cc}m&D_{eo}\\ D_{oe}&m\end{array}\right)$$

$$\left(\begin{array}[]{cc}m^{2}-D_{eo}D_{oe}&0\\ 0&m^{2}-D_{oe}D_{eo}\end{array}\right)\left(\begin{array}[]{c}y_{e}\\ y_{o}\end{array}\right)=\left(\begin{array}[]{c}\Phi_{e}\\ \Phi_{o}\end{array}\right).$$

$$x=M^{\prime }y.$$

Full text

We present a new general numerical method to compute many Green’s functions for complex non-singular matrices within one iteration process. Our procedure applies to matrices of structure $A=D-m$, with $m$ proportional to the unit matrix, and can be integrated within any Krylov subspace solver. We can compute the derivatives $x^{(n)}$ of the solution vector $x$ with respect to the parameter $m$ and construct the Taylor expansion of $x$ around $m$.

We demonstrate the advantages of our method using a minimal residual solver. Here the procedure requires $1$ intermediate vector for each Green’s function to compute. As real life example, we determine a mass trajectory of the Wilson fermion matrix for lattice QCD. Here we find that we can obtain Green’s functions at all masses $\geq m$ at the price of one inversion at mass $m$.

1 Introduction

Lattice Quantum Chromodynamics demands for enormous compute power and efficient numerical algorithms in order to solve huge sparse systems of linear equations. These linear systems originate from a discretized differential operator, the fermionic matrix, with varying gluonic background fields acting as coefficient functions. The latter can be generated by Monte Carlo methods. From the solutions, i.e. the Green’s functions or propagators, one can determine the properties of hadrons, as e.g. spectrum, weak decay constants or weak matrix elements [1, 2, 3].

Recently, in Ref. [4], it has been demonstrated that it can be very advantageous to exploit the structure of the Wilson fermion matrix $M=m-D$ [5], which is widely used in lattice calculations. In particular, the computation of so-called mass trajectories, i.e., multiple inversions of $m-D$ for a set of ‘bare masses’111The ’physical’ bare masses differ from $m$ by a constant shift, which is not relevant in this context $m_{0},m_{1},\ldots,m_{N}$ with $m_{i}<m_{i+1}$, could be carried out within one iteration process, using the quasi minimum residual algorithm (QMR) [6, 7, 8]. This useful property is a direct consequence of the non-symmetric Lanczos process [9], where the Lanczos part had to be gone through only once as it depends solely on $D$. Therefore, the method presented in Ref. [4] is restricted to QMR. Unfortunately the application of the QMR-MULT algorithm [4] is rather memory intensive: QMR itself needs $6$ intermediate vectors of the size of one propagator component and requires to store $3$ additional vectors for each of the members of the mass-trajectory. Computations of propagators as of today, however, tend to be memory bounded even on the largest distributed memory machines available.

In this paper, we propose a new numerical approach to compute the entire mass trajectory simultaneously. This method is not restricted to a special solver such as QMR and requires only a modest additional amount of workspace.

For definiteness and simplicity we present our idea as applied within the minimal residual algorithm (MR), one of the favorite standard solvers in lattice QCD that belongs to the class of the general conjugate residual methods[10, 11]. MR is a GCR algorithm with recursion depth [math].

Our key observation is that the derivatives $x^{(n)}=d^{n}x/dm^{n}$ of the solution $x$ can be related recursively to the coefficients that occur within the iteration. Thus the Taylor expansion of the propagator $x(m+\Delta m)$ around $m$ can be computed almost for free during the iteration process.

We find numerically that the procedure renders $x(m_{i}>m_{0})$ automatically with higher accuracy than $x(m_{0})$. The accuracy of the latter is determined as usual by the stopping residual’s norm. This enables us to calculate the entire mass trajectory at once, provided that we compute at $m_{0}$. On the other side, the procedure can provide good estimates for propagators with $x(m_{i}<m_{0})$ to be used as educated guesses for further inversions.

We note that our method can be applied to the case of even-odd-preconditioning [12]. For general sources, however, two separate inversions are necesssary. The additional expense of a factor of two easily pays off against large numbers of Green’s functions to be computed.

As a by-product we are in principle in the position to compute all derivatives $x^{(n)}$ up to order $n=N$, the number where the iteration stops.

The paper is organized as follows: In section 2.1, we present the MR algorithm followed by the computation of the derivatives in section 2.2. The construction of the Taylor series is explained in section 2.3. Furthermore we will present the results of numerical tests that deal with the inversion of the Wilson fermion matrix on a real full QCD configuration at $\beta=5.6$ and $\kappa=0.1575$ on a $16^{3}\times 32$ lattice.

2 M3R: Multiple Masses Minimal Residual

2.1 Standard Minimal Residual Algorithm

Consider the linear system:

[TABLE]

with $A$ being a complex nonsingular matrix. The search direction in the Minimal Residual Algorithm is taken to be the direction of the residual. At every iteration step the following operation is performed:

[TABLE]

The coefficient $\alpha_{i}$ is determined by a minimization of the 2-norm of the new residual $r_{i+1}$:

[TABLE]

2.2 Computation of the Derivatives

It is possible to extract in addition to the solution an arbitrary number of derivatives $x^{(n)}=d^{n}x(m)/dm^{n}$ at almost no extra computational cost, provided the matrix is of the special form of equation 1 and the source $r_{0}$ is independent of $m$.

$n$-fold differentiation of eq. 1 with respect to $m$ leads to a linear equation for the $n$-th derivative

[TABLE]

Let $x_{i}$ denote the approximation of the solution $x$ after $i$ iterations of the Minimal Residual algorithm. We are interested in a similar approximation for the 1st derivative $x_{i}^{(1)}$. Therefore we demand eq. 4 for $n=1$ to hold at each iteration step:

[TABLE]

It is obvious that the approximants $x_{i}^{(1)}$, defined in this way, converge towards $x^{(1)}$ as $x_{i}$ approaches $x$. The next step is to calculate the recursion relation for $x_{i}^{(1)}$. To this end we make use of the recursion relation for $x_{i}$ itself which is particularly simple in the case of the Minimal Residual algorithm:

[TABLE]

We then find the formal recursion relation

[TABLE]

At this stage we replace the true solution $x$ by it’s best approximation $x_{N}$.

[TABLE]

A summation over $i$ then leads to the following expression for $x_{N+1}^{(1)}$:

[TABLE]

It is important to start the inversion with an initial guess $x_{0}=0$ since the nonsingularity of $A$ then implies that $x_{0}^{(1)}$ vanishes as well222A nontrivial starting guess for $x_{0}$ would require a starting guess for $x_{0}^{(1)}$ computed from $Ax_{0}^{(1)}=x_{0}$ at the price of an extra inversion. After a resummation of the r.h.s we find

[TABLE]

Note that we can recursively build up the r.h.s. of this equation in terms of the coefficients

[TABLE]

at the price of only one extra vector-update per iteration step

[TABLE]

One can extend this procedure to higher derivatives by use of the following recursion relation for the generalized coefficients:

[TABLE]

In this way one ends up with the general updating-step for the derivatives

[TABLE]

We remark that for every additional derivative only one extra vector of storage is necessary.

2.3 Recombination of the Taylor-series

Usually one is interested in the solution of eq. 1 for a given set of different values of $m$. In this section we will show that it is possible to perform a summation over all available approximate derivatives of $x(m)$ and hence obtain a Taylor extrapolation to a different value of the parameter $m$.

We would like to perform a summation of the N-truncated Taylor series to an extrapolated solution $x^{E}\equiv x(m+\Delta m)$:

[TABLE]

This Taylor series can be constructed recursively by use of equation 13:

[TABLE]

where we have introduced scale-factors

[TABLE]

In order to compute these scale-factors $f_{i}$ we recast eq. 12 in recursive form:

[TABLE]

Inserted into equation 16 this yields:

[TABLE]

We close with the comment that the entire procedure is selfconsistent up to terms of $O(m^{N+1})$ and $O(\Delta m^{N+1})$.

2.4 Template for the $M^{3}R$-Algorithm

Next we write down a template for our Multiple Masses Minimal Residual $(M^{3}R)$ algorithm. For shortness we apply one extrapolated solution $x^{E}(\Delta m)$, more solutions can be included easily.

[TABLE]

Note that the standard overrelaxation parameter $\omega$ can be included as usual.

3 Applications in Lattice QCD

3.1 The Wilson Fermion Matrix

In 4 dimensional lattice Quantum Chromodynamics the Wilson fermion operator has the structure:

[TABLE]

with $m$ being the bare-mass and

[TABLE]

the off-diagonal hopping term. The coefficients ${U_{\mu}(x)}$ are SU(3) background gauge fields. Greens Functions in QCD are the solutions of the linear equation

[TABLE]

with $\Phi$ being a $\delta$-source or a ’smeared’ source [13, 4].

3.2 Results

The numerical tests were performed on the Quadrics Q4 parallel computer at Wuppertal University and the QH2 at IfH/DESY Zeuthen. We test the M3R algorithm on a full-QCD configuration of size $16^{3}\times 32$ at $\beta=5.6$ and the following mass trajectory:

[TABLE]

We calculated the propagator on the lightest mass, corresponding to $\kappa=0.1575$ and extrapolated with the M3R - technique to the six heavier masses. Our stopping criterion was

[TABLE]

with $||.||$ denoting the 2-norm and $r_{i}$ being the true residual $||M(m+\Delta m)x_{i}-r_{0}||$. In figure 1 the true residuals of the 7 inversions against the number of M3R iterations are shown.

It should be noted that the true residuals of the extrapolated heavier masses are always smaller than those corresponding to the lightest mass, for which the MR is actually inverting. We even observe a faster convergence with increasing $\Delta m$, despite the fact that we used a truncated Taylor-expansion to perform the extrapolation to larger masses. All residuals reach a plateau slightly below $10^{-5}$, which reflects the 32bit precision of the Quadrics computer.

The gain factor for our example on a 256-node QH2 is a factor of 3 in CPU-time compared to seperate inversions at all the masses of the trajectory. This gain factor, of course, depends on the number of Green’s functions that fit into memory and the mass range of the trajectory.

3.3 Even-Odd Preconditioning

In lattice quantum chromodynamics it has become standard to use the so called even-odd or red-black ordering scheme for the lattice-sites. This ordering scheme induces a block structure of the fermion matrix, which can be exploited to improve its condition number.

If one denotes the even lattice sites with $x_{e}$, the odd sites with $x_{o}$, equation 21 reads

[TABLE]

Right preconditioning with the matrix

[TABLE]

yields two seperate matrix equations, each with a smaller condition number resulting in a performance gain of a factor 2-3:

[TABLE]

The vector $y$ can be transformed back to $x$ using one additional matrix multiplication:

[TABLE]

The M3R-Template can be applied as shown above to the even-odd preconditioned system using $m^{2}$ as extrapolation parameter.

We test the preconditioned version of the algorithm on the same configuration, which has been reordered into even-odd form. We prefer to use a $\delta$-source on an even lattice-point, because it enables us to restrict the application to the upper half of the system saving a factor of two in CPU time. Using the trajectory of seven masses from section 3.2 we obtain the results shown in fig. 2.

The behaviour remains qualitatively the same. Due to the smaller condition number of the matrix the convergence speed increases by approximately a factor of 2.

The gain factor compared to a sequential MR inversion on every mass is equal to 3. For a general source, however, one has to carry out the inversion on both systems of equation 26 reducing the above gain factor for our example to 1.5 .

4 Summary and Outlook

We have presented a new numerical method to compute a trajectory of solutions $x$, of the linear system $Ax=r$ exploiting the shifted structure $A=m-D$ of the linear system. We have shown that in a realistic setup in the framework of lattice Quantum Chromodynamics the method appears to be numerically stable as well as more efficient than standard methods.

For our setup of masses the gain factor varies between 1.5 and 3 depending on the source. Furthermore it depends strongly on the set of masses to compute. Note that we have chosen a large range of masses to demonstrate the convergence properties of the algorithm. In a realistic QCD setup the gain factor would be higher. This holds in particular for situations where one would like to make use of the Feynman-Hellmann theorem with respect to the quark mass dependence as in the case of the Pion-Nucleon-Sigma term [14] in Quenched QCD.

The method can be generalized in many different directions: In principle it can be embedded into any known iterative solver like Conjugate Gradient (CG) or Generalized Conjugate Residual (GCR). Furthermore it seems possible to derive similar algorithms for matrices with more complex parameter-dependencies than the one discussed in this paper.

5 Acknowledgements

We appreciate the support of the staff at the computer center at IfH/DESY Zeuthen.