f11-3

11.3EigenvaluesandEigenvectorsofaTridiagonalMatrix 469Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).CITED REFERENCES AND FURTHER READING: Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press), §5.1. [1] Smith, B.T., et al. 1976, Matrix Eigensystem Routines — EISPACK Guide , 2nd ed., vol. 6 of Lecture Notes in Computer Science (New York: Springer-Verlag). Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra , vol. II of Handbook for Automatic Com- putation(New York: Springer-Verlag). [2] 11.3 Eigenvalues and Eigenvectors of a Tridiagonal Matrix Evaluation ofthe Characteristic Polynomial Onceouroriginal,real,symmetricmatrixhasbeenreducedtotridiagonalform, onepossiblewaytodetermineitseigenvaluesistofindtherootsofthecharacteristicpolynomial p n(λ)directly. Thecharacteristicpolynomialofatridiagonalmatrixcan be evaluated for any trial value of λby an efficient recursion relation (see [1], for example). The polynomials of lower degree producedduring the recurrenceform a Sturmian sequence that can be used to localize the eigenvalues to intervals on the real axis. A root-finding method such as bisection or Newton’s method can thenbe employed to refine the intervals. The corresponding eigenvectors can then be found by inverse iteration (see §11.7). Procedures based on these ideas can be found in [2,3]. If, however, more than a small fraction of all the eigenvalues and eigenvectors are required, then the factorization method next considered is much more efficient. The QRand QLAlgorithms The basic idea behind the QRalgorithm is that any real matrix can be decomposed in the form A=Q·R (11.3.1 ) whereQis orthogonal and Ris upper triangular. For a general matrix, the decompositionisconstructedbyapplyingHouseholdertransformationstoannihilatesuccessive columns of Abelow the diagonal (see §2.10). Now consider the matrix formed by writing the factors in (11.3.1) in the opposite order: A /prime=R·Q (11.3.2 ) SinceQis orthogonal,equation(11.3.1)gives R=QT·A. Thus equation (11.3.2) becomes A/prime=QT·A·Q (11.3.3 ) 470 Chapter11. EigensystemsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).We see that A/primeis an orthogonal transformation of A. You can verify that a QRtransformation preserves the following properties of amatrix: symmetry,tridiagonalform,andHessenbergform(tobedefinedin §11.5). There is nothing special about choosing one of the factors of Ato be upper triangular; one could equally well make it lower triangular. This is called the QL algorithm, since A=Q·L (11.3.4 ) whereLis lower triangular. (The standard, but confusing, nomenclature RandL stands for whether the rightorleftof the matrix is nonzero.) RecallthatintheHouseholderreductiontotridiagonalformin §11.2,westarted in the nth (last) column of the original matrix. To minimize roundoff, we then exhorted you to put the biggest elements of the matrix in the lower right-hand corner, if you can. If we now wish to diagonalize the resulting tridiagonal matrix,theQLalgorithm will have smaller roundoff than the QRalgorithm, so we shall useQLhenceforth. TheQLalgorithm consists of a sequence of orthogonaltransformations: A s=Qs·Ls As+1=Ls·Qs (=QT s·As·Qs)(11.3.5 ) The following (nonobvious!) theorem is the basis of the algorithm for a general matrixA: (i)IfAhaseigenvaluesofdifferentabsolutevalue |λi|,thenAs→[lower triangular form] as s→∞. The eigenvalues appear on the diagonal in increasing order of absolute magnitude. (ii) If Ahas an eigenvalue |λi|of multiplicity p, As→[lower triangular form] as s→∞, except for a diagonal block matrix of order p, whose eigenvalues →λi. The proof of this theorem is fairly lengthy; see, for example, [4]. The workload in the QLalgorithm is O(n3)per iteration for a general matrix, which is prohibitive. However, the workload is only O(n)per iteration for a tridiagonal matrix and O(n2)for a Hessenberg matrix, which makes it highly efficient on these forms. Inthissectionweareconcernedonlywiththecasewhere Aisareal,symmetric, tridiagonal matrix. All the eigenvalues λiare thus real. According to the theorem, if any λihas a multiplicity p, then there must be at least p−1zeros on the sub- and superdiagonal. Thus the matrix can be split into submatrices that can be diagonalized separately, and the complication of diagonal blocks that can arise in the general case is irrelevant. In the proof of the theorem quoted above, one finds that in general a super- diagonal element converges to zero like a(s) ij∼/parenleftbiggλi λj/parenrightbiggs (11.3.6 ) Although λi<λ j, convergence can be slow if λiis close to λj. Convergence can be accelerated by the technique of shifting:I fkis any constant, then A−k1has 11.3EigenvaluesandEigenvectorsofaTridiagonalMatrix 471Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).eigenvalues λi−k. If we decompose As−ks1=Qs·Ls (11.3.7 ) so that As+1=Ls·Qs+ks1 =QT s·As·Qs(11.3.8 ) then the convergence is determined by the ratio λi−ks λj−ks(11.3.9 ) The idea is to choose the shift ksat each stage to maximize the rate of convergence. A good choice for the shift initially would be ksclose to λ1, the smallest eigenvalue. Then the first row of off-diagonalelements would tend rapidly to zero. However, λ1is not usually known a priori. A very effective strategy in practice (althoughthere is no proof that it is optimal) is to compute the eigenvaluesof the leading 2×2diagonal submatrix of A. Then set k sequal to the eigenvalue closer to a11. Moregenerally,supposeyouhavealreadyfound r−1eigenvaluesof A. Then youcandeflatethematrixbycrossingoutthe first r−1rows andcolumns,leaving A= 0 ··· ··· 0 ··· 0... d r er... ... erdr+1 ··· 0 dn−1en−1 0 ··· 0en−1dn (11.3.10 ) Choose k sequaltotheeigenvalueoftheleading 2×2submatrixthatiscloserto dr. One can show that the convergence of the algorithm with this strategy is generally cubic (and at worst quadratic for degenerate eigenvalues). This rapid convergence is what makes the algorithm so attractive. Note that with shifting, the eigenvalues no longer necessarily appear on the diagonal in order of increasing absolute magnitude. The routine eigsrt(§11.1) can be used if required. As we mentionedearlier,the QLdecompositionofa generalmatrixis effected byasequenceofHouseholdertransformations. Foratridiagonalmatrix,however,itis moreefficienttouseplanerotations Ppq. Oneusesthesequence P12,P23,...,Pn−1,n to annihilate the elements a12,a23,...,a n−1,n. By symmetry, the subdiagonal elements a21,a32,...,a n,n−1will be annihilated too. Thus each Qsis a product of plane rotations: QT s=P(s) 1·P(s) 2···P(s) n−1 (11.3.11 ) 472 Chapter11. EigensystemsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).wherePiannihilates ai,i+1. Notethatit is QTinequation(11.3.11),not Q,because we defined L=QT·A. QLAlgorithmwithImplicitShifts The algorithm as described so far can be very successful. However, when the elements of Adiffer widely in order of magnitude, subtracting a large ks from the diagonal elements can lead to loss of accuracy for the small eigenvalues. This difficulty is avoided by the QLalgorithm with implicit shifts . The implicit QLalgorithm is mathematically equivalent to the original QLalgorithm, but the computation does not require ks1to be actually subtracted from A. The algorithm isbased on the following lemma: If Aisa symmetric nonsingular matrix andB=QT·A·Q, whereQis orthogonal and Bis tridiagonal with positive off-diagonal elements, then QandBare fully determined when the last row of QTis specified. Proof: LetqT idenote the ith row vector of the matrix QT. Thenqiis the ith column vector of the matrixQ. The relation B·QT=QT·Acan be written  β 1γ1 α2β2γ2 ... αn−1βn−1γn−1 αn βn · q T 1 qT2 ... qT n−1 qT n = q T 1 qT2 ... qT n−1 qT n ·A (11.3.12 ) The nth row of this matrix equation is α nqT n−1+βnqT n=qT n·A (11.3.13 ) SinceQis orthogonal, qT n·qm=δnm (11.3.14 ) Thus if we postmultiply equation (11.3.13) by qn, we find βn=qT n·A·qn (11.3.15 ) which is known since qnis known. Then equation (11.3.13) gives αnqT n−1=zT n−1 (11.3.16 ) where zT n−1≡qT n·A−βnqT n (11.3.17 ) is known. Therefore α2 n=zT n−1zn−1, (11.3.18 ) or αn=|zn−1| (11.3.19 ) and qT n−1=zT n−1/α n (11.3.20 ) (where αnis nonzero by hypothesis). Similarly, one can show by induction that if we know qn,qn−1,...,qn−jand the α’s,β’s, and γ’s up to level n−j, one can determine the quantities at level n−(j+1 ). To apply the lemma in practice, suppose one can somehow find a tridiagonal matrix As+1such that As+1=QT s·As·Qs (11.3.21 ) whereQT sis orthogonal and has the same last row as QT sin the original QLalgorithm. ThenQs=QsandAs+1 =As+1. 11.3EigenvaluesandEigenvectorsofaTridiagonalMatrix 473Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Now, in the original algorithm, from equation (11.3.11) we see that the last row of QT s is the same as the last row of P(s) n−1. But recall that P(s) n−1is a plane rotation designed to annihilate the (n−1,n)element of As−ks1. A simple calculation using the expression (11.1.1) shows that it has parameters c=dn−ks/radicalbig e2n+(dn−ks)2,s =−en−1/radicalbig e2n+(dn−ks)2(11.3.22 ) The matrix P(s) n−1·As·P(s)T n−1is tridiagonal with 2 extra elements:  ··· ××× ××× x ××× x×× (11.3.23 ) We must now reduce this to tridiagonal form with an orthogonal matrix whose last row is [0,0,..., 0,1]so that the last row of QT swill stay equal to P(s) n−1. This can be done by a sequence of Householder or Givens transformations. For the special form of the matrix(11.3.23), Givens is better. We rotate in the plane (n−2,n−1)to annihilate the (n−2,n) element. [By symmetry, the (n, n−2)element will also be zeroed.] This leaves us with tridiagonal form except for extra elements (n−3,n−1)and (n−1,n−3). We annihilate these with a rotation in the (n−3,n−2)plane, and so on. Thus a sequence of n−2 Givens rotations is required. The result is that QT s=QT s=P(s) 1·P(s) 2···P(s) n−2·P(s) n−1 (11.3.24 ) where the P’s are the Givens rotations and Pn−1is the same plane rotation as in the original algorithm. Then equation (11.3.21) gives the next iterate of A. Note that the shift ksenters implicitly through the parameters (11.3.22). Thefollowingroutine tqli(“Tridiagonal QLImplicit”),basedalgorithmically on the implementations in [2,3], works extremely well in practice. The number of iterations for the first few eigenvalues might be 4 or 5, say, but meanwhile the off-diagonal elements in the lower right-hand corner have been reduced too. The latereigenvaluesareliberatedwithverylittlework. Theaveragenumberofiterations per eigenvalue is typically 1.3−1.6. The operation count per iteration is O(n), with a fairly large effective coefficient, say, ∼20n. The total operation count for the diagonalization is then ∼20n×(1.3−1.6)n∼30n2. If the eigenvectors are required, the statements indicated by comments are included and there is an additional, much larger, workload of about 3n3operations. SUBROUTINE tqli(d,e,n,np,z) INTEGER n,np REAL d(np),e(np),z(np,np) C USES pythag QL algorithm with implicit shifts, to determine the eigenvalues and eigenvectors of a real, symmetric, tridiagonal matrix, or of a real, symmetric matrix previously reduced by tred2 §11.2. dis a vector of length np. On input, its first nelements are the diagonal elements of the tridiagonal matrix. On output, it returns the eigenvalues. The vector einputs the sub- diagonal elements of the tridiagonal matrix, with e(1) arbitrary. On output eis destroyed. When finding only the eigenvalues, several lines may be omitted, as noted in the comments.If the eigenvectors of a tridiagonal matrix are desired, the matrix z(nbynmatrix stored innpbynparray) is input as the identity matrix. If the eigenvectors of a matrix that has been reduced by tred2 are required, then zis input as the matrix output by tred2 .I n either case, the kth column of zreturns the normalized eigenvector corresponding to d(k) . INTEGER i,iter,k,l,m REAL b,c,dd,f,g,p,r,s,pythag 474 Chapter11. EigensystemsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).do11i=2,n Convenient to renumber the elements of e. e(i-1)=e(i) enddo 11 e(n)=0. do15l=1,n iter=0 1d o 12m=l,n-1 Look for a single small subdiagonal element to split the matrix. dd=abs(d(m))+abs(d(m+1)) if (abs(e(m))+dd.eq.dd) goto 2 enddo 12 m=n 2 if(m.ne.l)then if(iter.eq.30)pause ’too many iterations in tqli’ iter=iter+1 g=(d(l+1)-d(l))/(2.*e(l)) Form shift. r=pythag(g,1.) g=d(m)-d(l)+e(l)/(g+sign(r,g)) This is dm−ks. s=1.c=1.p=0. do 14i=m-1,l,-1 A plane rotation as in the original QL,f o l - lowed by Givens rotations to restore tridi-agonal form.f=s*e(i) b=c*e(i) r=pythag(f,g) e(i+1)=rif(r.eq.0.)then Recover from underflow. d(i+1)=d(i+1)-p e(m)=0. goto 1 endif s=f/r c=g/rg=d(i+1)-pr=(d(i)-g)*s+2.*c*b p=s*r d(i+1)=g+pg=c*r-b C Omit lines from here ... do 13k=1,n Form eigenvectors. f=z(k,i+1)z(k,i+1)=s*z(k,i)+c*f z(k,i)=c*z(k,i)-s*f enddo 13 C ... to here when finding only eigenvalues. enddo 14 d(l)=d(l)-pe(l)=ge(m)=0. goto 1 endif enddo 15 return END CITED REFERENCES AND FURTHER READING: Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), pp. 331–335. [1] Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra , vol. II of Handbook for Automatic Com- putation(New York: Springer-Verlag). [2] 11.4HermitianMatrices 475Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-1992 by Numerical Recipes Software. Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).Smith, B.T., et al. 1976, Matrix Eigensystem Routines — EISPACK Guide , 2nd ed., vol. 6 of Lecture Notes in Computer Science (New York: Springer-Verlag). [3] Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §6.6.6. [4] 11.4 Hermitian Matrices The complex analog of a real, symmetric matrix is a Hermitian matrix, satisfying equation (11.0.4). Jacobi transformationscan be used to find eigenvaluesandeigenvectors,asalsocanHouseholderreductiontotridiagonalformfollowedby QLiteration. Complexversions of the previousroutines jacobi,tred2, and tqli are quite analogousto their real counterparts. Forworkingroutines,consult [1,2]. An alternative, using the routines in this book, is to convert the Hermitian problem to a real, symmetric one: If C=A+iBis a Hermitian matrix, then the n×ncomplex eigenvalue problem (A+iB)·(u+iv)=λ(u+iv)( 11.4.1 ) is equivalent to the 2n×2nreal problem /bracketleftbigg A−B BA/bracketrightbigg ·/bracketleftbigg u v/bracketrightbigg =λ/bracketleftbigg u v/bracketrightbigg (11.4.2 ) Note that the 2n×2nmatrix in (11.4.2) is symmetric: AT=AandBT=−B ifCis Hermitian. Corresponding to a given eigenvalue λ, the vector /bracketleftbigg −v u/bracketrightbigg (11.4.3 ) is also an eigenvector, as you can verify by writing out the two matrix equa- tions implied by (11.4.2). Thus if λ1,λ2,...,λ nare the eigenvalues of C, then the 2neigenvalues of the augmented problem (11.4.2) are λ1,λ1,λ2,λ2,..., λn,λ n; each, in other words, is repeated twice. The eigenvectors are pairs of the formu+ivandi(u+iv);thatis,theyarethesameuptoaninessentialphase. Thus wesolvetheaugmentedproblem(11.4.2),andchooseoneeigenvalueandeigenvectorfromeachpair. Thesegivetheeigenvaluesandeigenvectorsoftheoriginalmatrix C. Working with the augmented matrix requires a factor of 2 more storage than the original complex matrix. In principle, a complex algorithm is also a factor of 2 more efficient in computer time than is the solution of the augmented problem. In practice, most complex implementations do not achieve this factor unless they arewritten entirely in real arithmetic. (Good library routines always do this.) CITED REFERENCES AND FURTHER READING: Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra , vol. II of Handbook for Automatic Com- putation(New York: Springer-Verlag). [1] Smith, B.T., et al. 1976, Matrix Eigensystem Routines — EISPACK Guide , 2nd ed., vol. 6 of Lecture Notes in Computer Science (New York: Springer-Verlag). [2]