We describe in some detail our numerical treatment of Neuberger's lattice Dirac operator as implemented in a practical application. We discuss the improvements we have found to accelerate the numerical computations and give an estimate of the expense when using this operator in practice

B Bunk

F Niedermayer

KG Wilson

L Fox

M Löscher

P Hasenfratz

P Hernández

PH Ginsparg

RG Edwards

S Chandrasekharan

T Blum

WH Press

Y Saad

arXiv

Pilar Hernández

Karl Jansen

Laurent Lellouche

Crossref

A Numerical Treatment of Neuberger’s Lattice Dirac Operator

We describe in some detail our numerical treatment of Neuberger's lattice Dirac operator as implemented in a practical application. We discuss the improvements we have found to accelerate the numerical computations and give an estimate of the expense when using this operator in practice.We describe in some detail our numerical treatment of Neuberger's lattice Dirac operator as implemented in a practical application. We discuss the improvements we have found to accelerate the numerical computations and give an estimate of the expense when using this operator in practice

Hernandez, Pilar

Jansen, Karl

Lellouch, Laurent

CERN Document Server

CERN/TH-99-376CPT-99/PE.3916LAPTH-Conf-771/99A numerical treatment of Neuberger’s latticeDirac operatorPilar Hernandezy, Karl JansenzCERN, 1211 Geneva 23, Switzerlandand Laurent LellouchxLAPTH, Chemin de Bellevue, B.P. 110,F74941 Annecy-le-Vieux Cedex, FranceJanuary 13, 2000AbstractWe describe in some detail our numerical treatment of Neuberger’s lat-tice Dirac operator as implemented in a practical application. We discussthe improvements we have found to accelerate the numerical computationsand give an estimate of the expense when using this operator in practice.∗Talk given by K.J. at the “Interdisciplinary Workshop on Numerical Challenges to LatticeQCD”, Wuppertal, August 22-24,1999†On leave from Departamento de F´ısica Teo´rica, Universidad de Valencia.‡Heisenberg Foundation Fellow§On leave from Centre de Physique The´orique, CNRS Luminy, F-13288 Marseille Cedex 9,France.11 Lattice formulation of QCDToday, we believe that the world of quarks and gluons is described theoreti-cally by quantum chromodynamics (QCD). This model shows a number of non-perturbative aspects that cannot be adequately addressed by approximation schemessuch as perturbation theory. The only way to evaluate QCD, addressing both itsperturbative and non-perturbative aspects at the same time, is lattice QCD. Inthis approach the theory is put on a 4-dimensional Euclidean space-time latticeof nite physical length L, with a non-vanishing value of the lattice spacing a.Having only a nite number of grid points, physical quantities can be computednumerically by solving a high-dimensional integral by Monte Carlo methods, mak-ing use of importance sampling.The introduction of a lattice spacing regularizes the theory and is an intermediatestep in the computation of physical observables. Eventually, the regularizationhas to be removed and the value of the lattice spacing has to be sent to zero toreach the target theory, i.e. continuum QCD. In fact, in conventional formulationsof lattice QCD [1], the introduction of the lattice spacing renders the theoryon the lattice somewhat dierent from the continuum analogue and a numberof properties of the continuum theory are only very dicult and cumbersometo establish in the lattice regularized theory. One of the main reasons for thisdiculty is that in conventional lattice QCD the regularization breaks a particularsymmetry of the continuum theory, which plays a most important role there,namely chiral symmetry.However, the last few years have seen a major breakthrough in that we now haveformulations of lattice QCD that have an exact lattice chiral symmetry [5]. Inthis approach, many properties of continuum QCD are preserved even at a non-vanishing value of the lattice spacing [2, 3, 4, 5, 6]. This development followedthe rediscovery [7] of the so-called Ginsparg{Wilson (GW) relation [8] which isfullled by any operator with the exact lattice chiral symmetry of [5]. It is not theaim of this contribution to discuss the physics consequences of the GW relation.We have to refer the interested reader to reviews [9, 10] about these topics. Herewe would like to discuss the numerical treatment of a particular lattice operatorthat satises the GW relation, namely Neuberger’s solution [4]. This solution hasa complicated structure and is challenging to implement numerically. Thus, the2large theoretical advantage of an operator satisfying the GW relation must beweighed against the very demanding computational eort required to implementit.This contribution is organized as follows. After discussing Neuberger’s latticeDirac operator we want to show how we evaluated the operator in our practi-cal application [11] and what kind of improvements we found to accelerate thenumerical computations. For alternative ideas for improvements, see the contri-butions of H. Neuberger [12] and A. Borici [13] to this workshop. We nally givesome estimates of the computational expense of using Neuberger’s operator.2 Neuberger’s lattice Dirac operatorThe operator we have used acts on elds (complex vectors) (x) where x =(x0; x1; x2; x3) and the xµ;  = 0; 1; 2; 3, are integer numbers denoting a 4-dimensionalgrid point in a lattice of size N4 with N = L=a. The elds (x) carry in additiona \colour" index  = 1; 2; 3 as well as a \Dirac" index i = 1; 2; 3; 4. Hence,  isa N4  3  4 complex vector.In order to reach the expression for Neuberger’s operator we rst introduce thematrix AA = 1 + s− a2γµ(rµ +rµ− arµrµ} ; (1)where rµ and rµ are the nearest-neighbour forward and backward derivatives,the precise denition of which can be found in the Appendix. The parameters is to be taken in the range of jsj < 1 and serves to optimize the localizationproperties [14] of Neuberger’s operator, which is then given byD =1an1− A (AyA−1/2o : (2)Through the appearance of the square root in eq. (2), all points on the lattice areconnected with one another, giving rise to a very complicated, multi-neighbouraction. However, the application of D to a vector  will only contain applicationsof A or AyA on this vector. Since these matrices are sparse, as only nearest-neighbour interactions are involved, we will never have to store the whole matrix.In the computation of physical quantities, the inverse of D, applied to a givenvector, is generically needed. Hence one faces the problem of having to compute3a vector X = D−1, with  a prescribed vector (the \source") as required by theparticular problem under investigation. Fortunately, a number of eciently work-ing algorithms for computing X = D−1 are known, such as conjugate gradient,BiCGstab, or variants thereof [16]. In conventional approaches to lattice QCD anoperator ~D is used that is very similar to the matrix A in eq. (1). Computing thevector ~X = ~D−1 requires a number niter of iterations of some particular method,say BiCGstab. Employing Neuberger’s operator D in computing X = D−1, itturns out that the number of iterations needed is of the same order of magnitudeas when using ~D. At the same time, in each of these iterations, the square roothas to be evaluated. When this is done by some polynomial approximation, it isfound that the required degree of this polynomial is roughly of the same order asthe number of iterations needed for computing the vector X. Hence, with respectto the conventional case, the numerical eort is squared and the price to pay forusing the operator D is high.On the other hand, any solution of the Ginsparg{Wilson relation gives us a toolby which particular problems in lattice QCD can be studied, which would beextremely hard to address with conventional approaches. It is for these casesthat the large numerical eort is justied, but clearly, we would like to haveclever ideas coming from areas such as Applied Mathematics, to decrease thenumerical expense or even overcome this bottleneck.3 Approximation of(AyA−1=2For computing the square root that appears in eq. (2), we have chosen a Cheby-shev approximation [15] by constructing a polynomial Pn,(x) of degree n, whichhas an exponential convergence rate in the interval x 2 [; 1]. Outside this inter-val, convergence is still found but it will not be exponential. The advantages ofusing this Chebyshev approximation are the well-controlled exponential t accu-racy as well as the possibility of having numerically very stable recursion relations[17] to construct the polynomial, allowing for large degrees. In order to have anoptimal approximation, it is desirable to know the lowest and the highest eigen-value of AyA. A typical example of the eigenvalues of AyA is shown in g. 1,where we show the 11 lowest eigenvalues as obtained on a number of congu-4Figure 1: Monte Carlo time evolution of the eleven lowest eigenvalues of AyA at = 5:85. The lowest eigenvalue for each conguration is the open square.rations using the Ritz functional method [19]. There is a wide spread and verylow-lying eigenvalues appear. Choosing  to be the value of the lowest of theseeigenvalues would result in a huge degree n of the polynomial Pn,. We thereforecomputed O(10) lowest-lying eigenvalues of AyA as well as their eigenfunctionsand projected them out of the matrix AyA. The approximation is then only per-formed for the matrix with a reduced condition number, resulting in a substantialdecrease of the degree of the polynomial. In addition, we computed the highesteigenvalue of AyA and normalized the matrix A such that kAyAk . 1.Since our work [11], aiming at the physical question of spontaneous chiral sym-5metry breaking in lattice QCD, has been one of the rst of its kind, we wantedto exclude possible systematic errors and demanded a very high precision for theapproximation to the square root:kX − Pn,(AyA)AyAPn,(AyA)Xk2=k2Xk2 < 10−16 (3)where X is a gaussian random vector. In our practical applications we xed thisprecision beforehand and set  to be the 11th lowest eigenvalue of AyA. This thendetermines the degree of the polynomial n and hence our approximation Dn tothe exact Neuberger operator D. We checked that the precision we required forthe approximation of the square root is directly related to the precision by whichthe GW relation itself is fullled. Choosing n such that the accuracy in eq. (4)is reached results ink [γ5Dn + Dnγ5 −Dnγ5Dn] Xk2=kXk2  10−16 : (4)In addition, we nd that the deviations from the exact GW relation decreaseexponentially fast with increasing n.4 The inverse of Neuberger’s operatorAs mentioned above, in physics applications a vector D−1 has to be computed,with  a prescribed source vector. Not only is the computation of this vector verycostly, there also appears to be a conceptual problem: in inspecting the lowesteigenvalue of DynDn, very small eigenvalues are often found as shown in g. 2.These very small eigenvalues belong to a given chiral sector of the theory, i.e.their corresponding eigenfunctions  are eigenfunctions of γ5 with γ5 = . Infact, these modes play an important physical role as they are associated withtopological sectors of the theory [3, 2, 5, 4].As far as the practical applications are concerned, it is clear that in the presenceof such a small eigenvalue, the inversion of Dn will be very costly, as the conditionnumber of the problem is then very high. In order to address this problem, wefollowed two strategies:(i) We compute the lowest eigenvalue of DynDn and its eigenfunction (usingagain the Ritz functional method [19]) and if it is a zero mode {in which6Figure 2: Monte Carlo time evolution of the lowest eigenvalue of DynDn. Theeigenvalues belong to given chiral sectors of the theory denoted as  for chiralityplus (full squares) and minus (open squares). Data are obtained at  = 5:85choosing s = 0:6. Whenever there is a zero mode of DynDn, the value of thelowest eigenvalue is set to 10−8.7case it is also a zero mode of Dn{ we project this mode out of Dn andinvert only the reduced matrix; this is then well conditioned, as the verysmall eigenvalues appear to be isolated. In this strategy, the knowledge ofthe eigenfunction must be very precise and an accuracy of approximatingthe square root as indicated in eq. (3) is mandatory.(ii) Again we determine the lowest eigenvalue of DynDn and the chirality of thecorresponding zero mode, if there is any. We then make use of the fact thatDynDn commutes with γ5. This allows us to perform the inversion in thechiral sector without zero modes. In this strategy, the accuracy demanded ineq. (3) could be relaxed and this strategy, which essentially follows ref. [18],is in general much less expensive than following strategy (i).However, even adopting strategy (ii), solving the system DnX =  is still costly.We therefore tried two ways of improving on this. We rst note that instead ofsolving h1− A=pAyAiX =  (5)we can equally well solve hAy −pAyAiX = Ay : (6)In practice, however, we found no real advantage in using the formulation ofeq. (6). We have further considered two acceleration schemes.Scheme (a)We choose two dierent polynomials (now approximatingpAyA and not theinverse) Pn, and Pm,, m < n, such thatPn, = Pm, +  (7)with  a \small" correction. Then we haveAy − Pn,−1=Ay − Pm, −−1h1 +(Ay − Pm,−1i (Ay − Pm,−1: (8)This leads us to the following procedure of solving DnX = :8(1) rst solve(Ay − Pm,Y =  ; (9)(2) then solve(Ay − Pm,X0 =  + Y ; (10)(3) use X0 as a starting vector to nally solve(Ay − Pn,X =  : (11)The generation of the starting vector X0 in steps (1) and (2) is only a smalloverhead. In g. 3 we plot the relative residuum 2stop = kDnX − k2=kXk2 as afunction of the number of applications of Dn. In this case n = 100 and m = 30.We show the number of applications of the matrix Dn for the case of a randomstarting vector (dotted line) and the case where X0 was generated according tothe above procedure (solid line). The gain is of approximately a factor of two.Scheme (b)In the second approach, we use a sequence of polynomials to solve DnX = . Tothis end we rst solve Ay − Pm1,X1 = Ay (12)by choosing a polynomial Pm1, and a stopping criterion for the solver (1)stop suchthatm1 < n; (1)stop > stop : (13)The value of (1)stop is chosen such that it is roughly of the same order of magnitudeas the error that the polynomial of degree m1 itself induces. The solution X19Figure 3: The residuum as a function of the number of applications of the matrixDn. The dotted line corresponds to a random starting vector. The solid line toa starting vector generated following scheme (a).10is then used as a starting vector for the next equation, employing a polynomialPm2, and stopping criterion (2)stop withm1 < m2 < n; (1)stop > (2)stop > stop : (14)This procedure is then repeated until we reach the desired polynomial Pn, andstopping criterion stop to solve the real equationAy − Pn,X = Ay : (15)As for scheme (a), we gain a factor of about two in the numerical eort. Wenally remark that some rst tests using the scheme proposed in [13] resulted ina similar performance gain as the two schemes presented above.In table 1 we give a typical example of the expense of a simulation followingstrategy (ii). We list both the cost of computing the lowest eigenvalue of DynDnin terms of the number of iterations to minimize the Ritz functional [19] and thenumber of iterations to solve DnX = . In both applications, a polynomial ofdegree n is used to approximate the square root. The numbers in table 1 indicatethat a quenched calculation, employing Neuberger’s operator, leads to a compu-tational cost that is comparable with a dynamical simulation using conventionaloperators.N n nev ninvert8 190 170 8010 250 325 20012 325 700 300Table 1: N is the number of lattice sites along a side of the hypercube; n, thedegree of polynomial; nev, the number of iterations required to obtain the lowesteigenvalue of DynDn; and ninvert, the number of iterations necessary to computeX = D−1n .115 ConclusionsThe theoretical advance that an exact chiral symmetry brings to lattice gaugetheory is accompanied by the substantial increase in numerical eort that isrequired to implement operators satisfying the GW relation. Thus, while theNielsen{Ninomiya theorem has been circumvented, the \no free lunch theorem"has not. Whether alternative formulations, such as domain wall fermions, canhelp in this respect remains to be seen.AppendixWe give here the explicit denitions needed in eq. (1). The forward and backwardderivatives rµ; rµ act on a vector (x) asrµ(x) = 1a[U(x; )(x + a^)− (x)]rµ(x) =1a(x) − U(x− a^; )−1(x− a^) ;where ^ denotes the unit vector in direction . The (gauge) eld U(x; ) 2 SU(3)lives on the links connecting lattice points x and x + a^ and acts on the colourindex  = 1; 2; 3 of the eld . Finally, the Dirac matrices γµ;  = 0; 1; 2; 3 arehermitean 44 matrices acting on the Dirac index i of the eld . Their explicitform is given byγµ = 0 eµeyµ 0!(16)withe0 = 1; ek = −ik (17)and1 = 0 11 0!2 = 0 −ii 0!3 = 1 00 −1!: (18)12With the choice of the γ matrices given above, the matrix γ5 = γ0γ1γ2γ3 isdiagonal and given byγ5 = 1 00 −1!: (19)We nally note that whenever repeated indices appear, they are summed over.References[1] K.G. Wilson, Phys. Rev. D10 (1974) 2445.[2] P. Hasenfratz, V. Laliena and F. Niedermayer, Phys.Lett.B427 (1998) 125.[3] P. Hasenfratz, Nucl. Phys. B525 (1998) 401.[4] H. Neuberger, Phys. Lett. B417 (1998) 141 and B427 (1998) 353.[5] M. Lu¨scher, Phys. Lett. B428 (1998) 342.[6] S. Chandrasekharan, Phys. Rev. D60 (1999) 074503.[7] P. Hasenfratz, Nucl. Phys. B (Proc.Suppl.) 63A-C (1998) 53.[8] P.H. Ginsparg and K.G. Wilson, Phys. Rev. D25 (1982) 2649.[9] F. Niedermayer, Nucl. Phys. B (Proc.Suppl.) 73 (1999) 105.[10] T. Blum, Nucl.Phys.B (Proc.Suppl.) 73 (1999) 167.[11] P. Hernandez, K. Jansen and L. Lellouch, CERN preprints CERN-TH/99-197, hep-lat/9907022 and CERN-TH/99-273, hep-lat/9909026.[12] H. Neuberger, hep-lat/9910040.[13] A. Borici, hep-lat/9910045.[14] P. Hernandez, K. Jansen and M. Lu¨scher, Nucl. Phys. B552 (1999) 363.[15] W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery, NumericalRecipes, Second Edition, Cambridge University Press, Cambridge, 1992.13[16] Y. Saad, Iterative methods for sparse linear systems, PWS Publishing Com-pany, Boston, 1996.[17] L. Fox and I. B. Parker, Chebyshev polynomials in numerical analysis, OxfordUniversity Press, London, 1968.[18] R.G. Edwards, U.M. Heller and R. Narayanan, Phys. Rev. D59 (1999)094510.[19] B. Bunk, K. Jansen, M. Lu¨scher and H. Simma, DESY report (September1994); T. Kalkreuter and H. Simma, Comp.Phys.Comm. 93 (1996) 33.14

A numerical treatment of Neuberger's lattice Dirac operator

A numerical treatment of Neuberger's lattice Dirac operator

Abstract

Similar works

Full text

Available Versions

Crossref

CERN Document Server