The Lanczos algorithm for matrix tridiagonalisation suffers from strong
numerical instability in finite precision arithmetic when applied to evaluate
matrix eigenvalues. The mechanism by which this instability arises is well
documented in the literature. A recent application of the Lanczos algorithm
proposed by Bai, Fahey and Golub allows quadrature evaluation of inner products
of the form $\psi^\dagger g(A) \psi$. We show that this quadrature evaluation
is numerically stable and explain how the numerical errors which are such a
fundamental element of the finite precision Lanczos tridiagonalisation
procedure are automatically and exactly compensated in the Bai, Fahey and Golub
algorithm. In the process, we shed new light on the mechanism by which roundoff
error corrupts the Lanczos procedureComment: 3 pages, Lattice 99 contributio

Alan Irving

Cahill

Christopher Johnston

Davis

Duncan

Eamonn Cahill

James Sexton

Parlett

English

arXiv

Crossref

Numerical stability of Lanczos methods

The Lanczos algorithm for matrix tridiagonalisation suffers from strong numerical instability in finite precision arithmetic when applied to evaluate matrix eigenvalues. The mechanism by which this instability arises is well documented in the literature. A recent application of the Lanczos algorithm proposed by Bai, Fahey and Golub allows quadrature evaluation of inner products of the form $\psi^\dagger g(A) \psi$. We show that this quadrature evaluation is numerically stable and explain how the numerical errors which are such a fundamental element of the finite precision Lanczos tridiagonalisation procedure are automatically and exactly compensated in the Bai, Fahey and Golub algorithm. In the process, we shed new light on the mechanism by which roundoff error corrupts the Lanczos procedure

Cahill, Eamonn

Irving, Alan

Johnston, Christopher

Sexton, James

The UKQCD Collaboration

University of Liverpool Repository

arXiv:hep-lat/9909131v1  17 Sep 1999Numerical Stability of Lanczos MethodsLTH457TCDMATH 99-10hep-lat/9909131Eamonn Cahill, a Alan Irving, b Christopher Johnson, b James Sexton, a and The UKQCDCollaborationa School of Mathematics, Trinity College, Dublin 2, Ireland and Hitachi Dublin Laboratory, Dublin 2,Irelandb Theoretical Physics Division, Department of Mathematical Sciences, University of Liverpool, PO Box147, Liverpool L69 3BX, UKThe Lanczos algorithm for matrix tridiagonalisation suffers from strong numerical instability in finite precisionarithmetic when applied to evaluate matrix eigenvalues. The mechanism by which this instability arises is welldocumented in the literature. A recent application of the Lanczos algorithm proposed by Bai, Fahey and Goluballows quadrature evaluation of inner products of the form ψ† · g(A) ·ψ. We show that this quadrature evaluationis numerically stable (accurate to machine precision) and explain how the numerical errors which are such afundamental element of the finite precision Lanczos tridiagonalisation procedure are automatically and exactlycompensated in the Bai, Fahey and Golub algorithm. In the process, we shed new light on the mechanism bywhich roundoff error corrupts the Lanczos procedure.1. THE BAI, FAHEY AND GOLUBMETHODThe determinant of the fermion matrix playsa central role in dynamical Lattice QCD sim-ulations. The direct approach [1] to evaluat-ing fermion determinants typically makes use ofthe Lanczos procedure for eigenvalue estimationof Hermitian matrices. Given an n × n Hermi-tian matrix A, and an orthonormal seed vectorf1, the Lanczos procedure is an iterative processwhich (in exact arithmetic) generates a sequenceof orthogonal vectors fi, i = 1, . . . ,m. If we de-fine an n × m matrix Fm = (f1, . . . , fm) whosecolumns are given by the vectors fi, then we have:F†mFm = Im, Tm = F†mAFm is tridiagonal, and theeigenvalues of Tm converge to the eigenvalues ofA as m→ n.An alternative stochastic approach to evaluat-ing determinants was proposed in 1996 by Bai,Fahey and Golub (BFG) [2] who observe that avector-matrix-vector of the form ψ†·g(A)·ψ can beexpressed as an integral of the function g(·) over amodified spectral measure. An m-point Gaussianquadrature approximation to this integral is thengiven by ψ† ·g(A)·ψ ≈∑mi=1 wig(θi). The remark-able result [3] applied in the BFG method is thatthe abscissa, θi, of this quadrature rule are givenby the eigenvalues of the tridiagonal matrix Tmwhich is generated by a Lanczos procedure ap-plied with seed vector ψ, and that the weights,wi, are given as the squares of the first compo-nents of the corresponding eigenvectors of Tm.Figure 1 shows an example of the results whichthe BFGmethod generates on a 123×24 quenchedgauge configuration. The configuration is aquenched configuration with β = 5.7 producedwithin the current UKQCD research program. Arandom Gaussian fermion vector, ψ, was gener-ated, and the BFG method applied to evaluateψ† · log(M†κMκ) · ψ for Wilson hopping parame-ter κ = .1650 as a function of the order m of theLanczos tridiagonal matrix and of the resultingGaussian quadrature approximation. To generatean estimate of the “true” value for this matrix-vector inner-product for comparison, we arbitrar-ily averaged the values for order 1991 to order2000. The plot then shows the log of the differ-ence of the m-th order BFG estimate and this“true” value plotted against m. The plot showsa rapid convergence to the “true” value (conver-gence has occurred at about order 200), followedby what is clearly fluctuations at the level of algo-21e-201e-1011e+10400 800 1200 1600 2000ErrorFigure 1. Convergence history for the evaluationof ψ† · log(M†κMκ) ·ψ on a 123×24 quenched con-figuration as a function of the number of Lanczositerations.rithm precision. Note in particular that the longconverged region is perfectly flat, and shows nodrift. Checks for a toy model, and the mathe-matical analysis which we briefly describe here,indicate that the convergence is in fact, conver-gence to the exact answer, and that finite preci-sion arithmetic introduces no systematic error.2. THE EFFECT OF FINITE PRECI-SION ARITHMETICAt first sight, Figure 1 might be interpreted asgood phenomenological evidence that the BFGmethod is reliable and robust. A second sighthowever should disturb a person familiar withthe Lanczos procedure at finite precision. In fi-nite precision arithmetic, the Lanczos proceduregenerates clusters of almost degenerate eigenval-ues of the tridiagonal matrix Tm where only asingle eigenvalue would have occurred in an in-finite precision calculation. These clusters couldpotentially destroy the convergence of the BFGquadrature sum to the correct infinite precisionresult since the corresponding abscissa are in-cluded multiple times rather than singly in thatsum. We have found however, that when clustersarise, the corresponding weights in the quadra-ture rule adjust so that the net contribution ofthe cluster of multiple eigenvalues and weights isequal to that of the corresponding single eigen-value and weight which would arise in infiniteprecision.This effect is illustrated in Table 1 which showsthe weights corresponding to a given cluster fora Lanczos procedure applied to a model problemwhich has been chosen so all quantities can be cal-culated both numerically and exactly and directcomparisons can be made. At iteration 40, thecluster in question contains only a single eigen-value, and the relative error between the value ofthe weight as generated by the Lanczos procedureand an exact calculation is seen to be O(10−14).By iteration 160, the cluster expands to includefour almost degenerate eigenvalues (which typi-cally differ only at about O(10−8)). The individ-ual weights are seen now to differ quite signifi-cantly from the correct value, but the sum of thethree weights does match the correct value, andthe contribution of the cluster to the quadratureis seen to be independent of the size of the cluster.3. AN EXPLANATIONThe eigenvalue clustering effect which is wellknown [4], and the invariance of weights whichwe have just described imposes strict constraintson the mechanism by which roundoff error cor-rupts the Lanczos procedure. At finite precision,the Lanczos procedure still generates a sequenceof vectors fi, i = 1 . . .m. However these vectorsno longer remain orthogonal, and no longer tridi-agonalise A. We have instead that1. F†mFm 6= Im,2. F†mAFm = Hm, symmetric but not tridiago-nal,3. The tridiagonal matrix, Tm which is, inpractice, used to estimate eigenvalues of A isgiven by the truncation of Hm to tridiagonalform (i.e. the all off-tridiagonal elements ofHm are set to zero).Ignoring the very small differences be-tween eigenvalues in a degenerate cluster,the eigenvalue spectrum of Tm is given as3Table 1This table illustrates the invariance of the sum of the weights for an eigenvalue cluster generated from afinite precision Lanczos procedure.m Weights in Cluster Sum of Weights Relative Error40 0.121274864637753313E-02 0.121274864637753313E-02 0.250E-14120 0.239992265988143331E-04 0.121274864637753270E-02 0.286E-140.298384753074422057E-030.890364666704296305E-03160 0.128211917795651532E-05 0.121274864637753248E-02 0.303E-140.118900026992380976E-030.470001178219875870E-030.622565321987319131E-03{θ1(n1), . . . , θj(nj), θj+1(1), . . . θp(1)}, where niis the multiplicity of the cluster with eigenvalueθi, j counts the non-unit-multiplicity clusters,and p counts the different distinct eigenvalues.We denote the ni eigenvectors with degenerateeigenvalue θi by s(k)i , k = 1 . . . ni. The n-rowby m column matrix Fm relates these eigenvec-tors to “Ritz” vectors yi which are n-componentvectors in the space spanned by the columns ofFm as yi = Fms(k)i . A well established result[4] for the Lanczos procedure is that, in a non-unit-multiplicity cluster, the eigenvectors s(k)i allgenerate the same Ritz vector yi (approximately).Further this Ritz vector is a true eigenvector ofthe matrix A, and the corresponding eigenvalueis a true eigenvalue of A (approximately). Thematrix Fm is therefore seen to be rank-deficientsince there are independent linear combinationsof columns which equal the same vector.We have found [5] that the best way to ex-pose this rank deficiency is to execute a singu-lar value decomposition of Fm. Using this singu-lar value decomposition, it is possible to find aspecial set of orthonormal basis vectors for thespace spanned by the columns of Fm. In thisspecial basis, the matrix Hm takes a block diag-onal form with blocks corresponding to the non-unit-multiplicity eigenvalues of Tm having entriesat all positions equal to the corresponding eigen-value. When Hm is truncated to diagonal form inthe normal way, these blocks in the special basisrepresentation become diagonal (the effect of thetruncation is to set all the non-diagonal entries forthe non-unit-multiplicity blocks equal to zero).We thus obtain the usual degenerate eigenvalueclusters of the truncated matrix Tm. The reasonthat the Lanczos procedure at finite precision pro-duces good eigenvalue estimates is thus seen to bea combination of the specially structured natureof the space spanned by the columns of Fm, andof the very special way in which the truncation todiagonal form acts. The invariance of the sum ofweights for a cluster is then explained by the ob-servation that this sum gives the projection of thestarting vector ψ seeding the Lanczos procedureonto the corresponding Ritz vector.REFERENCES1. A. Duncan, E. Eichten, R. Roskies, H.Thacker, Phys. Rev. D 60 (1999), hep-lat/9902015.2. Z. Bai, M. Fahey and G. Golub, “Some largescale matrix computation problems”, Techni-cal Report SCCM-95-09, School of Engineer-ing, Stanford University, 1995.3. P. Davis and P. Raninowitz, “Methods of Nu-merical Integration”, Academic Press, NewYork, 1984.4. B. Parlett, The Symmetric Eigenvalue Prob-lem, 2nd edition, Prentice-Hall, EnglewoodCliffs, NJ, 1998.5. E. Cahill, “Krylov Subspace Methods and theCalculation of Matrix Determinants”, Ph.D.Thesis, School of Mathematics, Trinity Col-lege, Dublin, 1999.

Numerical Stability of Lanczos Methods

Numerical Stability of Lanczos Methods

Abstract

Similar works

Full text

Available Versions

Crossref

University of Liverpool Repository