Comparison of Inversion Algorithms for Wilson Fermions on the CM5

TL;DR

This paper compares iterative algorithms like MR, CG, BiCGγ₅, and BiCGStab for inverting the Wilson fermion Dirac matrix on the CM-5, finding preconditioned over-relaxed MR and BiCGStab to be most effective and robust.

Contribution

It provides a comparative analysis of various iterative inversion algorithms for Wilson fermions on the CM-5, highlighting the efficiency and robustness of preconditioned over-relaxed MR and BiCGStab.

Findings

Preconditioned over-relaxed MR is simple, memory-efficient, and comparable to BiCGStab in performance.

MR and BiCGStab are equally robust, insensitive to initial solutions and round-off errors.

The study was conducted on 32^3×64 lattices with Wilson fermion action, with and without clover term.

Abstract

This talk presents results of a comparitive study of iterative algorithms like minimal residue ($MR$) and conjugate gradient ($CG$, $BiCG\gamma_5$, and \bicgstab) used for inverting the Dirac matrix $M$. The tests were done on the Connection Machine CM-5 using $32^3 \times 64$ lattices. The fermion action used is of Wilson type, both with and without the clover term. The overall conclusion is that preconditioned over-relaxed $MR$ is the simplest, uses the least memory, and is comparable in performance to \bicgstab. We find these two algorithms to be equally robust, $i.e.$ insensitive to the starting solution and to round-off errors.

Tables

Table 1. Table 1: Operations per iteration in the various algorithms and number of stored vectors.

Algorithm	Mults	saxpy$+$	dot$+$	Storage
	by $M$	caxpy	cdot
$M{R}_2$	$2$	$0+2$	$1+1$	3
$M{R}_4$	$4$	$0+2$	$1+1$	3
$C{G}_2$	$4$	$1+2$	$2+0$	4
$BiCG{\gamma }_5$	$4$	$3+0$	$0+2$	3
${BiCGstab}_1$	$2$	$0+6$	$1+3$	7
${BiCGstab}_2$	$4$	$0+6$	$1+3$	7

Equations

$$M=1+\kappa D(r)$$

$$m_{q}={1\over 2a}\bigg{(}{1\over\kappa}-{1\over\kappa_{c}}\bigg{)}$$

$$(1+\kappa D(r))\chi =\varphi$$

$$(1-{\kappa }^2{D}^2)\chi$$

$$(1-{\kappa }^4{D}^4)\chi$$

$$doi=1,$$

$$\alpha$$

$${\chi }_i$$

$$r_i$$

$$enddo$$

$$doi=1,$$

$$\alpha$$

$${\chi }_i$$

$$r_i$$

$$\beta$$

$$g_i$$

$$enddo$$

$$doi=1,$$

$$\alpha$$

$${\chi }_i$$

$$r_i$$

$$\beta$$

$$g_i$$

$$enddo$$

$$doi=1,$$

$$\beta$$

$$p_i$$

$$v_i$$

$$\alpha$$

$$s_i$$

$$t_i$$

$${\omega }_i$$

$${\chi }_i$$

$$r_i$$

$$enddo$$

Full text

1 INTRODUCTION

The solution of the Dirac equation for an arbitrary source constitutes the dominant fraction of computer time used in the numerical simulation of both quenched and unquenched QCD. In quenched simulations, the Dirac propagator is used in the construction of correlation functions of observables involving quark fields. Depending on how many observables are measured using a given set of quark propagators, the matrix inversion is $40-70\%$ of the total CPU time. In full QCD simulations this number jumps to $>90\%$ as matrix inversion is also needed in the update. Large scale simulations (quenched and unquenched) currently consume the equivalent of many giga-flop years on a variety of supercomputers. Thus, fast solvers are the key to progress, and we consider algorithm improvements by $>20\%$ significant. In this talk we catalogue and compare the efficiency of commonly used algorithms (Minimal Residue($MR$), conjugate gradient ($CG$, $BiCG\gamma_{5}$, and $BiCGstab$)) for Wilson and Clover class of actions. Further details of the lattice generation, the calculation of quark propagators, and the analysis of the meson and baryon spectrum are given in Ref. [1].

2 TECHNICAL DETAILS

The Dirac matrix discretised $a\ la$ Wilson on a 4-dimensional hypercubic grid can be written as

[TABLE]

where the matrix $D(r)$ connects a given site to its 8 nearest neighbours, $r$ is the Wilson parameter introduced to remove the doubling problem, and $\kappa$ is a parameter in terms of which the quark mass is defined as

[TABLE]

where $a$ is the lattice spacing and $\kappa_{c}$ denotes the chiral limit at which the matrix becomes singular. The quark propagator $\chi$ is the solution to the set of linear equations

[TABLE]

where $\phi$ is an arbitrary source vector. For a delta function source it represents the propagation of a quark from that point to all other points on the lattice.

2.1 Properties of the Dirac matrix

The Dirac matrix satisfies the identity $\gamma_{5}M\gamma_{5}=M^{\dagger}$ which is retained by both the Wilson and Staggered discretized versions. As a result the eigenvalues of $M$ come in complex-conjugate pairs. For the staggered theory, they lie along a line in the imaginary direction at $m_{q}$. For the Wilson theory they are distributed as shown in Fig. 1. The eigenvalues lying along the real axis in the confined phase acquire an imaginary part at the chiral/deconfinement transition. I believe that these real eigenvalues play a central role in the rate of convergence, and a study of their distribution and change with preconditioning may provide clues to developing better algorithms.

2.2 Convergence Criteria

We define the remainder at the $n^{th}$ iteration to be $r_{n}=M\chi_{n}-\phi$ and in the figures plot the convergence condition ${\cal C}\equiv|r|^{2}/|\chi|^{2}$. On $32^{3}\times 64$ lattices, 32 bit IEEE precision restricts the precision to which ${\cal C}$ can be measured to $\approx 10^{-14}$ due to round-off errors. In the figures we plot the iterated residue, so we compare algorithms using the same preconditioned matrix. Also, since our final stopping criteria is with respect to the true residue, a comparison of the total number of multiplies by the original Dirac matrix $M$ is also meaningful, irrespective of the preconditioning.

2.3 Polynomial Preconditioning

The following pre-conditioned systems are used to solve for $\chi$ [2]

[TABLE]

where we define $M^{\prime}=1-\kappa D(r)$. We label these second and fourth order pre-conditioned systems by adding subscripts to the algorithm name. The second-order form is equivalent to red/black pre-conditioning and one can halve the number of operations by solving for only the red (or black) sites and reconstructing the other from these. On the CM5 the necessary reorganization of data arrays to make use of this simplification has not been done. The tests have been performed solving for the full matrix as they do not affect the number of matrix multiplies needed (the red and black systems are solved independently and simultaneously).

2.4 Staggered Fermions

The staggered version of the Dirac matrix has the form $m_{q}+iA$, where $A$ is hermitian. Thus, the hermitian product $M^{\dagger}M=(m_{q}^{2}+A^{2})$ is also second order preconditioned. For such systems $CG$ has proven to be the optimal algorithm [4].

3 ALGORITHMS

In this section we give the precise implementation of the various algorithms since the speed of convergence can depend on the details. We also give the performance characteristics of these algorithms, saving a comparitive study for the next section.

3.1 Minimal Residue

The minimal residue algorithm consists of the following steps:

[TABLE]

where $\alpha$ is a complex number of $O(1)$, and the order of the hermitian dot product is important. (The algorithm fails to converge if we use $\alpha^{*}$ instead of $\alpha$, and is about $50\%$ slower if we use real$(\alpha)$). For initial guess one can use $\chi_{0}=0$ or any other trial vector. The only restriction seems to be that $\langle Mr_{i-1},r_{i-1}\rangle$ should not be small. This situation arises for $\chi_{n}$ produced by the $CG$ class of algorithms. Thus, $MR$ should NOT be used following even a few $CG$ steps.

The parameter $\omega$ in Eq. 5 is an over-relaxation parameter. We find that, to within few percent, the convergence rate is insensitive to $\omega$ in the range $1.1-1.35$, and gives the envelope of best convergence. For $1.0\leq\omega\leq 1.08$ the convergence fluctuates between the $\omega=1$ (poorest) and $\omega=1.1$ (best) cases as shown in Fig. 2.

We find that higher orders of polynomial pre-conditioning extend the range of applicability (towards smaller quark masses) of the $MR$ algorithm. There is a $\sim 10\%$ degradation in CPU time with each order, $MR_{4}$ takes roughly $10\%$ more computer time but can be used for lighter quark masses [2].

3.2 Conjugate Gradient

The $CG$ method for a non-hermitian matrix is

[TABLE]

where both $\alpha,\beta$ are complex numbers and we set $g_{0}=M^{\dagger}r_{0}$. For $M$ we use the second order preconditioned (non-hermitian) matrix defined in 4. The $CG$ algorithm shows two distinct rates of convergence, a fast initial drop, and a slower asymptotic rate that is relevant for light quarks and makes $CG$ less competitive than $MR$ or $BiCGstab$. The only cases where we have seen some gain due to the fast initial convergence behavior of $CG$ is on switching from $MR$ or $BiCGstab$ at late stages of convergence, and even then only for light quarks. This is exemplified in Fig. 3.

3.3 $BiCG\gamma_{5}$

The $BiCG\gamma_{5}$ algorithm exploits the property $\gamma_{5}M\gamma_{5}=M^{\dagger}$ of Dirac fermions to yield a stabilized version of $CG$ algorithm for the matrix $M$ itself [4]

[TABLE]

where $\alpha$ and $\beta$ are real and $g_{0}=r_{0}$. This algorithm, in our opinion, has two drawbacks as shown in Fig. 4; it is sensitive to the initial $\chi_{0}$, for example the choice $\chi_{0}=0$ is bad as the algorithm actually blows up. A second pecularity of this algorithm is that the convergence shows very large fluctuations with a period of $\sim 10$ iterations. Overall, we do not see any advantage to using it since the asymptotic convergence rate is no better than $MR$ or $BiCGstab$.

Forcrand and Borici have raised the possibility that the problems could be due to round-off errors in the accumulations. For this purpose we ran our tests with dot products done in 64 bit precision. This change did not improve the convergence.

3.4 $BiCGstab$

The steps in the $BiCGstab$ algorithm are [5]

[TABLE]

where $\alpha,\beta,\omega$ are complex numbers initialized to $(1,0)$, $v_{0}=p_{0}=0$, and $\hat{r}_{0}=r_{0}=\phi-M\chi_{0}$. We have tested this algorithm for two cases of $M$, the original Dirac matrix and the second-order preconditioned version and label these as \hbox{BiCGstab}_{1} and \hbox{BiCGstab}_{2} respectively.

The sharp spikes at count 600 and 1200 in Fig. 5 result from refreshing the remainder. The data show that the convergence rate is not affected by these spikes. The sharp drop at count 1900 is caused by a switch to $CG$ and is apparent only for the smallest quark mass.

4 Clover action

In Fig. 6 we show the convergence behavior of $MR_{2}$ and $BiCGstab_{2}$ algorithms when inverting the tadpole improved clover action at $\beta=6.0$. The dependence on gauge configurations is exhibited by using 6 distinct lattices and the lightest quark mass ($m_{\pi}=0.185$). Our conclusion is that the performance of both algorithms is comparable. If one is concerned about the stability of $MR_{2}$, then one can use $MR_{4}$. The overhead of the clover term is $\approx 15\%$.

Acknowledgements

We thank A. Borici, C. Davies, and P. de Forcrand for discussions on the distribution of eigenvalues and for providing data. We thank S. Sharpe for comments on this work. R.G. is very grateful to Klaus Schilling for organizing a very informative meeting and for inviting him. We acknowledge the support of the Advanced Computing Lab at Los Alamos for these simulations.