# On the Neuberger overlap operator

## Abstract

We compute Neuberger’s overlap operator by the Lanczos algorithm applied to the Wilson-Dirac operator. Locality of the operator for quenched QCD data and its eigenvalue spectrum in an instanton background are studied.

PACS No.: 11.15Ha, 11.30.Rd, 11.30.Fs

Key Words: Lattice QCD, Chiral fermions, Algorithms

1. Although brute force calculations of the quenched lattice QCD with Wilson fermions have been able to approach the chiral limit , there are increased efforts to make the chiral symmetry exact on the lattice . <sup>1</sup><sup>1</sup>1For a recent review on the topic see .

There are different starting points to formulate lattice actions with exact lattice chiral symmetry, but all of them seem to obey the Ginsparg-Wilson condition :

$$\gamma _5D^1+D^1\gamma _5=a\gamma _5\alpha ^1,$$
(1)   
where $`a`$ is the lattice spacing, $`D`$ is the lattice Dirac operator and $`\alpha ^1`$ is a local operator and trivial in the Dirac space.

A candidate is the overlap operator of Neuberger :

$$D=1\gamma _5\text{sign}(H),H=\gamma _5(\sigma aD_W)$$
(2)   
where $`\sigma `$ is a shift parameter in the range $`(0,2)`$, which we have fixed at one and $`D_W`$ is the Wilson-Dirac operator,

$$D_W=\frac{1}{2}\underset{\mu }{}[\gamma _\mu (_\mu ^{}+_\mu )a_\mu ^{}_\mu ]$$
(3)   
and $`_\mu `$ and $`_\mu ^{}`$ are the nearest-neighbor forward and backward difference operators.

The locality has been shown for smooth background fields and no violation has been observed in quenched samples simulated at moderate couplings .

2. So far, all the methods devised to compute the overlap operator by usual iterative solvers have been based on (rational) polynomial approximations of the inverse square root or the $`sign`$ function . (Mathematical foundations of these methods are reviewed in .) But they may exceed the storage limits in some machines. This is not the case with Legendre  and Chebyshev  polynomials, which on the other hand are not optimal . <sup>2</sup><sup>2</sup>2After submission of this paper, the rational polynomial approximation method  was improved with respect to memory limits  by running twice the Conjugate Gradient (CG) iteration, as it is the case here (see below) for the Lanczos algorithm.

In the present work we propose a new method, which uses the outcome of the Lanczos algorithm on $`H`$. The Lanczos iteration is known to approximate the spectrum of the underlying matrix in an optimal way and, in particular, it requires a constant memory .

Let $`Q_n=[q_1,\mathrm{},q_n]`$ be the set of orthonormal vectors, such that

$$HQ_n=Q_nT_n,q_1=\rho _1b,\rho _1=1/b_2$$
(4)   
where $`T_n`$ is a tridiagonal and symmetric matrix. Here $`b`$ stands for an arbitrary vector.

By writing down the above decomposition in terms of the vectors $`q_i,i=1,\mathrm{},n`$ and the matrix elements of $`T_n`$, we arrive at a three term recurrence that allows to compute these vectors in increasing order, starting from the vector $`q_1`$. This is called the Lanczos algorithm, which constructs a basis for the so called Krylov subspace: $`span(b,Hb,\mathrm{},H^{n1}b)`$ .

In the last equation, it has been assumed that after $`n`$ steps of the Lanczos algorithm, the Krylov subspace remains invariant. The task is the computation of $`x=(H^2)^{1/2}b`$. Our method is based on the following observations: Let $`(H^2)^{1/2}=f(H^2)`$ be a matrix-valued function, for example Robert’s integral formula :

$$f(H^2)=\frac{2}{\pi }_0^{\mathrm{}}𝑑t(t^2+H^2)^1$$
(5)   
Then, clearly:

$$f(H^2)Q_n=Q_nf(T_n^2)$$
(6)   
Since, on the other hand,

$$b=Q_ne_1^{(n)}/\rho _1,$$
(7)   
where $`e_1^{(n)}`$ denotes the unit vector with $`n`$ elements in the direction $`1`$, we get:

$$x=f(H^2)b=Q_nf(T_n^2)e_1^{(n)}/\rho _1$$
(8)   
There are some remarks to be made here:

a) By applying the Lanczos iteration on $`H`$, the problem of computing $`(H^2)^{1/2}`$ reduces to the problem of computing $`(T_n^2)^{1/2}`$ which is typically a much smaller problem than the original one. It can be solved for example by using the full decomposition of $`T_n`$ in its eigenvalues and eigenvectors; in fact this is the method we have employed too, for its compactness and the small overhead for moderate $`n`$.

b) In the floating point arithmetic, there is a danger that once the Lanczos polynomial (algorithm) has approximated well some part of the spectrum, the iteration reproduces vectors which are rich in that direction . As a consequence, the orthogonality of the Lanczos vectors is spoiled with an immediate impact on the history of the iteration.

c) In general, there is no guarantee that the algorithm will converge at smaller $`n`$, unless $`n=\text{rank}(H)`$ in exact arithmetic . Therefore, for a given $`n`$ the equations (4) and (8) hold approximately.

Therefore, in practical implementations one should be satisfied with a stopping criterium such as:

$$\text{error}(n)=|(\rho _1Hx_n_2)^21|^{1/2}$$
(9)   
is made small enough.

It is worth writing down the error in terms of the Lanczos matrix; straightforward algebra gives:

$$\text{error}(n)=|\beta _nz_n^{(n)}|,z_n=(T_n^2)^{1/2}e_1^{(n)}$$
(10)   
where $`\beta _n`$ is the element $`(n+1,n)`$ of the matrix $`T_{n+1}`$ and $`z_n^{(n)}`$ is the last component of the vector $`z_n`$.

Since $`H`$ and $`(H^2)^{1/2}`$ are equally conditioned in 2-norm, we expect, that once the system $`Hx=b`$ is solved, the system $`(H^2)^{1/2}x=b`$ is also solved. In this context, it is desirable to compare the error (10) with the residual error of the original system, $`r_n`$. As before, in terms of the Lanczos matrix, it is given by:

$$\rho _1r_n_2=|\beta _ny_n^{(n)}|,y_n=T_n^1e_1^{(n)}$$
(11)   
As long as the orthogonality between Lanczos vectors is sufficiently maintained, equations (10-11) should hold to a good accuracy.

To implement the result (8), we first construct the Lanczos matrix and then compute $`z_n`$. By repeating the iteration, we compute Lanczos vectors and obtain the result. We saved the scalar products, though it was not necessary. If we call $`1/\rho _i,i=1,\mathrm{}`$ the norm of the residual error of the system $`Hx=b`$, it is easy to show that

$$\rho _{i+1}\beta _i+\rho _i\alpha _i+\rho _{i1}\beta _{i1}=0$$
(12)   
Therefore, we have the following algorithm for solving the system $`(H^2)^{1/2}x=b`$:

$$\begin{array}{c}\beta _0=0,\rho _1=1/b_2,q_0=o,q_1=\rho _1b\hfill \\ fori=1,\mathrm{}\hfill \\ v=Hq_i\hfill \\ \alpha _i=q_i^{}v\hfill \\ v:=vq_i\alpha _iq_{i1}\beta _{i1}\hfill \\ \beta _i=v_2\hfill \\ q_{i+1}=v/\beta _i\hfill \\ \rho _{i+1}=\frac{\rho _i\alpha _i+\rho _{i1}\beta _{i1}}{\beta _i}\hfill \\ if\frac{1}{|\rho _{i+1}|}<tol,n=i,endfor\hfill \\ \\ Set(T_n)_{i,i}=\alpha _i,(T_n)_{i+1,i}=(T_n)_{i,i+1}=\beta _i,otherwise(T_n)_{i,j}=0\hfill \\ z_n=(T_n^2)^{1/2}=U_n(\mathrm{\Lambda }_n^2)^{1/2}U_n^Te_1^{(n)}\hfill \\ \\ q_0=o,q_1=\rho _1b,x_0=o\hfill \\ fori=1,\mathrm{},n\hfill \\ x_i=x_{i1}+q_iz_n^{(i)}/\rho _1\hfill \\ v=Hq_i\hfill \\ v:=vq_i\alpha _iq_{i1}\beta _{i1}\hfill \\ q_{i+1}=v/\beta _i\hfill \end{array}$$
(13)   
where by $`o`$ we denote a vector with zero entries and $`U_n,\mathrm{\Lambda }_n`$ the matrices of the egienvectors and eigenvalues of $`T_n`$. Note that there are only four large vectors necessary to store: $`q_{i1},q_i,v,x_i`$.

Obviously, the memory doesn’t grow with $`n`$. This is not the case for the shifted CG iterations () needed to compute the (rational) polynomial approximation of $`(H^2)^{1/2}`$. Since there is a one to one connection between CG and Lanczos, such approximations on the original matrix are transfered to the corresponding Lanczos matrix . This should be contrasted with the exact computation of $`(T_n^2)^{1/2}`$.

To test the above analysis, we have performed simulations of SU(3) gauge theory at $`\beta =6.0`$ on a $`8^316`$ lattice and picked up an equilibrated configuration.

In Fig.1 we show the residual error $`|(\rho _1Hx_n_2)^21|^{1/2}`$ computed directly from $`x_n`$ and compare with the same quantity given in terms of the Lanczos matrix, i.e. $`|\beta _nz_n^{(n)}|`$. It fluctuates between two branches: the upper one corresponds to odd $`n`$, the number of matrix-vector multiplications, and the lower one to even values of $`n`$. For large $`n`$, the computed and estimated errors deviate from each other, which may indicate accumulation of roundoff errors in the computed residual error. Note that we have employed 64 bit precision.

For comparison, we show in Fig.2 the residual error $`\rho _1r_n_2`$ of the system $`Hx=b`$ as computed directly and from $`|\beta _ny_n^{(n)}|`$. Again, we have two branches as explained above, but here there is no distinction between the computed and estimated errors. The appearance of two branches is not surprising since we are dealing with a non-definite matrix $`H`$.

In Fig.3 we compare the computed residual errors of the systems $`Hx=b`$ and $`(H^2)^{1/2}x=b`$ (upper branches). They are the same most of the time, unless $`n`$ becomes large and deviations become clearer. This behavior shows that both systems are solved at the same time, which should serve us as a guide, because the computation of $`x_i`$ at each step $`i`$ is very demanding.

We have compared the efficiency of the above method with that of the rational polynomial approximation . First, we solved the system $`Hx=b`$ by both Lanczos and CG with an $`10^5`$ accuracy for the residual error. We needed the same number of multiplications with $`H^2`$, 226. To compute the inverse square root with the same accuracy ($`10^5`$), we applied then the above method and the rational polynomial method with $`N=30`$, the number of the terms in the approximating sum (). (Smaller $`N`$ will give a lower accuracy: in Fig. 4 we display the norm of the residual error (9) as a function of $`N`$.) Therefore, if we stored the Lanczos vectors as in the rational polynomial method, both methods need about the same amount of work. To avoid memory restrictions, we increase the work by a factor of two. <sup>3</sup><sup>3</sup>3As it was stated above, after submission of this paper, the rational polynomial approximation method  was changed in the same fashion, by increasing the work by a facor of two .

3. An immediate application of the above method is to check the locality of Neuberger’s overlap operator $`D`$.

We used 100 equilibrated SU(3) configurations at $`\beta =5.7,6.0`$ on a $`8^316`$ lattice. For each configuration we computed the absolute value of $`D`$ elements in the first column with space, spin and color indices fixed at one, a selection that suffices to look for violations of the locality.

The average over 100 configurations is plotted in Fig. 5. At both values of $`\beta `$, there was no single configuration to show an exceptional behavior: the maximum deviation is slightly above the mean value. The values of $`D`$ decrease rapidly with the time slices. For large $`t/a`$ but away from the center, <sup>4</sup><sup>4</sup>4Throughout the paper we have used periodic boundary conditions in all directions. they fall off exponentially.

Recently, it has been shown that for $`\beta =6.0,6.2,6.4`$ the occurrence of configurations with exceptionally small eigenvalues of $`H^2`$ becomes vanishingly small for $`12^4`$ and $`16^4`$ large lattices, whereas for sufficiently smooth gauge fileds the locality is guaranteed .

In fact at $`\beta =5.7`$ the complete spectrum of $`1D_W`$ is in the left half of the complex plane. <sup>5</sup><sup>5</sup>5I am grateful to Urs M. Heller for pointing out that normal spectroscopy with Wilson fermions can be done at $`\beta =5.7`$ and $`\kappa =0.1675`$ , the latter being greater than our corresponding $`\kappa =1/6`$. We computed partial spectra of $`\gamma _5H`$ by the implicitely restarted Arnoldi iteration  and found that indeed at this coupling there was no eigenvalue in the right half of the complex plane for all our 100 configurations. Therefore, there exist a $`\beta (5.7,6.0)`$ at which the operator $`1D_W`$ becomes singular. In this case, there is an unbounded number of near zero modes of $`H`$ and therefore $`D`$ is no longer local.

4. As another application we consider the lattice index theorem for an SU(2) instanton background on a small $`4^4`$ lattice. An instanton on the lattice can be prepared in various ways. We follow  and prepare an instanton with size $`\rho `$ in the center of the lattice in the singular gauge. The index of $`D`$ is given by:

$$\text{index}(D)=\frac{1}{2}\text{Tr}[\text{sign}(H)]$$
(14)   
Because of the O(a) lattice errors, we expect the instanton being observed for $`\rho a`$. We computed the smallest eigenvalues of $`D`$ by an (implicitly restarted) Arnoldi iteration with the non-converged Ritz values used as explicit shifts . We fixed the number of Arnoldi steps at 32 and have stopped the iteration when the next starting vector norm is smaller than $`10^6`$. We have checked the stability of the computed eigenvalues by increasing the cutoff beyond the number of the converged eigenvalues. The stability is observed unless the cutoff becomes too large, which means that a larger Arnoldi matrix should be employed.

We note that it is crucial for the eigenvalue computation to have a proper accuracy in the computation of $`D`$, which in our case has been set at a residual error norm (9) less than $`10^{10}`$.

We have computed eigenvalues for $`\rho /a[0.5,1.5]`$ with steps of $`0.1`$. For brevity, we show in Table 1 eigenvalues of a smaller set of $`\rho `$. For $`\rho a`$ there are two zero modes, whereas for $`\rho >a`$ there is a single zero mode present. As numerical accuracy is an issue here, we have perturbed the instanton background by applying a small fluctuating gauge field. The picture doesn’t change, but the zero modes for $`\rho a`$ become nearly zero modes with opposite chiralities.

We have also computed the eigenvalues exactly by standard QR algorithms. The instanton mode appears single in both approaches. While the values of the other smallest eigenvalues are reproduced exactly by the implicitly restarted Arnoldi algorithm, their multiplicity cannot be handled. Since $`D`$ is normal, the Arnoldi matrix is normal and tridiagonal and therefore irreducible, giving no information on the multiplicity. The latter is essential when $`D`$ has more exact zero modes and one must rely on block variants of the same algorithm.

Even on such a small lattice, the computation of the zero modes is not the fastest method to compute the topological charge. Tracing the crossings of the smallest eigenvalues of $`H`$ is more practical.

Nonetheless we note that an estimation of the topological charge can be made during the computation of $`D`$ as described in this work. Having computed the Lanczos matrix of $`H`$ one can estimate sign($`H`$), which for our SU(2) instanton gives excellent agreement for most of starting vectors $`b`$. In general, a separate study is needed to conclude on this method.

5. To conclude: we have computed with a new method the overlap operator based on the Lanczos algorithm applied on the Wilson-Dirac operator. Compared to the other methods , its main advantage is of being free from memory restrictions. Additionally, there is no approximation made in the computation of the inverse square root of the Lanczos matrix.

The locality of the overlap operator has been tested. We recommend to check it always before any other computation.

The computation of $`D`$ turns to be more difficult than $`D_W`$. However, the so-called classically perfect actions  may help substantially to work on moderate lattices. Further studies are needed for the dynamical implementation of $`D`$.

We are grateful to stimulating discussions with Ferenc Niedermayer on the topics covered by this work and to Roland Rosenfelder and Philippe de Forcrand for making critical remarks on this paper.

We are grateful to Herbert Neuberger for comments on the method presented in this work following the posting of its first version.

We thank PSI where this work was done and SCSC Manno for the allocation of computer time on the NEC SX4.