11.3 Eigenvalues and Eigenvectors of a Tridiagonal Matrix 475 f+=e[j]*a[i][j]; hh=f/(h+h); Form K,equation (11.2.11). for(j=1;j<1;j+)[ Form q and store in e overwriting p. f=a[i][j]; e[j]=gme[j]-hh*f; for (k=1;k<=j;k++) Reduce a,equation (11.2.13) a[j][k]-(f*e[k]+g*a[i][k]); else e[i]=a[i][1]; d[i]=h; /Next statement can be omitted if eigenvectors not wanted*/ 83g d[1]=0.0: granted for 19881992 e[1]=0.0: /Contents of this loop can be omitted if eigenvectors not 1.200 wanted except for statement d[i]=a[i][i];/ for (i=1;i<=n;i++){ Begin accumulation of transformation ma- 1=1-1; trices. if (a[i]){ This block skipped when i=1 for(j=1;j<=1;j+){ g=0.0; for(k=1;k<=1;k++) Use u and u/H stored in a to form P.Q. (North g+=a[i]k]*a[k][]: for(k=1;k<=1;k++) America computer, THE ak][j]-=g*a[k][1]; ART d[i]=a[i][i]; This statement remains. 9 Programs a[i][i]=1.0; Reset row and column of a to identity for(j=1;j<=1;j++)a[j][1]=a[i][j]=0.0; matrix for next iteration. IENTIFIC to dir CITED REFERENCES AND FURTHER READING: Golub,G.H.,and Van Loan,C.F.1989,Matrix Computations,2nd ed.(Baltimore:Johns Hopkins University Press),85.1.[1] COMPUTING(ISBN 198918920 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). 10-621 Wilkinson,J.H.,and Reinsch,C.1971,Linear Algebra,vol.Il of Handbook for Automatic Com- putation (New York:Springer-Verlag).[2] Numerical Recipes 43108 11.3 Eigenvalues and Eigenvectors of a (outside Tridiagonal Matrix North Software. Amer ying of visit website Evaluation of the Characteristic Polynomial machine Once our original,real,symmetric matrix has been reduced to tridiagonal form, one possible way to determine its eigenvalues is to find the roots of the characteristic polynomial pn(A)directly.The characteristic polynomial of a tridiagonal matrix can be evaluated for any trial value of A by an efficient recursion relation(see [1],for example).The polynomials of lower degree produced during the recurrence form a
11.3 Eigenvalues and Eigenvectors of a Tridiagonal Matrix 475 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). f += e[j]*a[i][j]; } hh=f/(h+h); Form K, equation (11.2.11). for (j=1;j<=l;j++) { Form q and store in e overwriting p. f=a[i][j]; e[j]=g=e[j]-hh*f; for (k=1;k<=j;k++) Reduce a, equation (11.2.13). a[j][k] -= (f*e[k]+g*a[i][k]); } } } else e[i]=a[i][l]; d[i]=h; } /* Next statement can be omitted if eigenvectors not wanted */ d[1]=0.0; e[1]=0.0; /* Contents of this loop can be omitted if eigenvectors not wanted except for statement d[i]=a[i][i]; */ for (i=1;i<=n;i++) { Begin accumulation of transformation mal=i-1; trices. if (d[i]) { This block skipped when i=1. for (j=1;j<=l;j++) { g=0.0; for (k=1;k<=l;k++) Use u and u/H stored in a to form P·Q. g += a[i][k]*a[k][j]; for (k=1;k<=l;k++) a[k][j] -= g*a[k][i]; } } d[i]=a[i][i]; This statement remains. a[i][i]=1.0; Reset row and column of a to identity for (j=1;j<=l;j++) a[j][i]=a[i][j]=0.0; matrix for next iteration. } } CITED REFERENCES AND FURTHER READING: Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns Hopkins 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 Computation (New York: Springer-Verlag). [2] 11.3 Eigenvalues and Eigenvectors of a Tridiagonal Matrix Evaluation of the Characteristic Polynomial Once our original, real, symmetric matrix has been reduced to tridiagonal form, one possible way to determine its eigenvalues is to find the roots of the characteristic polynomial pn(λ) directly. The characteristic polynomial of a tridiagonal matrix can be evaluated for any trial value of λ by an efficient recursion relation (see [1], for example). The polynomials of lower degree produced during the recurrence form a
476 Chapter 11.Eigensystems 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 then be 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.31.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 QR and QL Algorithms The basic idea behind the OR algorithm is that any real matrix can be 且品 decomposed in the form A=Q·R (11.3.1) where Q is orthogonal and R is upper triangular.For a general matrix,the E 2 decomposition is constructed by applying Householder transformations to annihilate successive columns of A below the diagonal (see 82.10). Now consider the matrix formed by writing the factors in (11.3.1)in the opposite order: A'=R.Q (11.3.2) 豆。2 Since Q is orthogonal,equation(11.3.1)gives R=Q.A.Thus equation(11.3.2) becomes A'=QF·A·Q (11.3.3) We see that A'is an orthogonal transformation of A. You can verify that a QR transformation preserves the following properties of o a matrix:symmetry,tridiagonal form,and Hessenberg form(to be defined in $11.5). Numerical There is nothing special about choosing one of the factors of A to be upper 43106 triangular;one could equally well make it lower triangular.This is called the QL algorithm,since A=Q· (11.3.4) where L is lower triangular.(The standard,but confusing,nomenclature R and L stands for whether the right or lefi of the matrix is nonzero.) Recall that in the Householder reduction to tridiagonal form in 811.2,we started 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. the QL algorithm will have smaller roundoff than the QR algorithm,so we shall use OL henceforth
476 Chapter 11. Eigensystems Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 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 then be 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 QR and QL Algorithms The basic idea behind the QR algorithm is that any real matrix can be decomposed in the form A = Q · R (11.3.1) where Q is orthogonal and R is upper triangular. For a general matrix, the decomposition is constructed by applying Householder transformations to annihilate successive columns of A below the diagonal (see §2.10). Now consider the matrix formed by writing the factors in (11.3.1) in the opposite order: A = R · Q (11.3.2) Since Q is orthogonal, equation (11.3.1) gives R = QT · A. Thus equation (11.3.2) becomes A = QT · A · Q (11.3.3) We see that A is an orthogonal transformation of A. You can verify that a QR transformation preserves the following properties of a matrix: symmetry, tridiagonal form, and Hessenberg form (to be defined in §11.5). There is nothing special about choosing one of the factors of A to be upper triangular; one could equally well make it lower triangular. This is called the QL algorithm, since A = Q · L (11.3.4) where L is lower triangular. (The standard, but confusing, nomenclature R and L stands for whether the right or left of the matrix is nonzero.) Recall that in the Householder reduction to tridiagonal form in §11.2, we started 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, the QL algorithm will have smaller roundoff than the QR algorithm, so we shall use QL henceforth
11.3 Eigenvalues and Eigenvectors of a Tridiagonal Matrix 477 The QL algorithm consists of a sequence of orthogonal transformations: As=Q。·Lg (11.3.5) As+1=LgQ。(=Q·Ag·Q) The following (nonobvious!)theorem is the basis of the algorithm for a general matrix A:(i)If A has eigenvalues of different absolute value il,then As[lower triangular form]as s-oo.The eigenvalues appear on the diagonal in increasing order of absolute magnitude.(ii)If A has an eigenvalue of multiplicity p, 81 As-[lower triangular form]as s-oo,except for a diagonal block matrix of order p,whose eigenvaluesAi.The proof of this theorem is fairly lengthy; see,for example,[4]. The workload in the QL algorithm 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. RECIPES I In this section we are concerned only with the case where A is a real,symmetric, tridiagonal matrix.All the eigenvalues A;are thus real.According to the theorem, if any Ai has a multiplicity p,then there must be at least p-1 zeros on the sub-and superdiagonal.Thus the matrix can be split into submatrices that can be d n 2 Press. 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 OF SCIENTIFIC( (11.3.6) 入 Although Ai<入，convergence can be slow if入:is close toλj.Convergence can be accelerated by the technique of shifring:If k is any constant,then A-k1 has eigenvalues入，-k.If we decompose 10-621 Ag-kg1=Qg·Lg (11.3.7) Numerical so that 431 Ag+1=Ls·Q。+ks1 Recipes (11.3.8) =Q·AaQ North then the convergence is determined by the ratio 入-ks 显 (11.3.9) 入)-ks The idea is to choose the shift ks at each stage to maximize the rate of convergence.A good choice for the shift initially would be k close to A1,the smallest eigenvalue.Then the first row of off-diagonal elements would tend rapidly to zero.However,A1 is not usually known a priori.A very effective strategy in practice(although there is no proof that it is optimal)is to compute the eigenvalues
11.3 Eigenvalues and Eigenvectors of a Tridiagonal Matrix 477 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). The QL algorithm consists of a sequence of orthogonal transformations: As = 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 matrix A: (i) If A has eigenvalues of different absolute value |λi|, then As → [lower triangular form] as s → ∞. The eigenvalues appear on the diagonal in increasing order of absolute magnitude. (ii) If A has 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 QL algorithm 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. In this section we are concerned only with the case where A is a real, symmetric, tridiagonal matrix. All the eigenvalues λi are thus real. According to the theorem, if any λi has a multiplicity p, then there must be at least p − 1 zeros 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 superdiagonal element converges to zero like a(s) ij ∼ λi λj s (11.3.6) Although λi < λj , convergence can be slow if λi is close to λj . Convergence can be accelerated by the technique of shifting: If k is any constant, then A − k1 has 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 ks at each stage to maximize the rate of convergence. A good choice for the shift initially would be k s close to λ1, the smallest eigenvalue. Then the first row of off-diagonal elements would tend rapidly to zero. However, λ1 is not usually known a priori. A very effective strategy in practice (although there is no proof that it is optimal) is to compute the eigenvalues
478 Chapter 11.Eigensystems of the leading 2 x 2 diagonal submatrix of A.Then set ks equal to the eigenvalue closer to a11. More generally,suppose you have already foundr-1 eigenvalues of A.Then you can deflate the matrix by crossing out the first r-1 rows and columns,leaving 「0 0 0 dr er A= (11.3.10) er dr+1 0 dn-1 En-1 10 0 en-1 dn 欲 Cam Choose ks equal to the eigenvalue of the leading 2 x 2 submatrix that is closer to 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(811.1) 人 Press. can be used if required. As we mentioned earlier,the QL decomposition of a general matrix is effected by a sequence of Householder transformations.For a tridiagonal matrix,however,it is more efficient to use plane rotations POne uses the sequence P12,P23,...,Pn-1.n SCIENTIFIC( to annihilate the elements a12,a23,...,an-1.n By symmetry,the subdiagonal elements a21,a32,...,an.will be annihilated too.Thus each Q is a product of plane rotations: Q=Pa.P9.…P@1 (11.3.11) where Pi annihilates ai.1.Note that it is Q in equation(11.3.11),not Q,because we defined L=Q.A. 子 Numerica 10.621 4310 QL Algorithm with Implicit Shifts (outside Recipes The algorithm as described so far can be very successful.However,when Software. 首 the elements of A differ widely in order of magnitude,subtracting a large k from the diagonal elements can lead to loss of accuracy for the small eigenvalues. This difficulty is avoided by the QL algorithm with implicit shifis.The implicit QL algorithm is mathematically equivalent to the original QL algorithm,but the computation does not require ks 1 to be actually subtracted from A. The algorithm is based on the following lemma:If A is a symmetric nonsingular matrix and B=QA.Q,where Q is orthogonal and B is tridiagonal with positive off-diagonal elements,then Q and B are fully determined when the last row of Q is specified.Proof: Let q denote the ith row vector of the matrix Q.Then q;is the ith column vector of the
478 Chapter 11. Eigensystems Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). of the leading 2 × 2 diagonal submatrix of A. Then set k s equal to the eigenvalue closer to a11. More generally, suppose you have already found r − 1 eigenvalues of A. Then you can deflate the matrix by crossing out the first r − 1 rows and columns, leaving A = 0 ··· ··· 0 ··· 0 . . . dr er . . . . . . er dr+1 ··· 0 dn−1 en−1 0 ··· 0 en−1 dn (11.3.10) Choose ks equal to the eigenvalue of the leading 2 × 2 submatrix that is closer to d r. 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 mentioned earlier, the QL decomposition of a general matrix is effected by a sequence of Householder transformations. For a tridiagonal matrix,however, it is more efficient to use plane rotations Ppq. One uses the sequenceP12, P23,..., Pn−1,n to annihilate the elements a12, a23,...,an−1,n. By symmetry, the subdiagonal elements a21, a32,...,an,n−1 will be annihilated too. Thus each Qs is a product of plane rotations: QT s = P(s) 1 · P(s) 2 ··· P(s) n−1 (11.3.11) where Pi annihilates ai,i+1. Note that it is QT in equation (11.3.11), not Q, because we defined L = QT · A. QL Algorithm with Implicit Shifts The algorithm as described so far can be very successful. However, when the elements of A differ widely in order of magnitude, subtracting a large k s from the diagonal elements can lead to loss of accuracy for the small eigenvalues. This difficulty is avoided by the QL algorithm with implicit shifts. The implicit QL algorithm is mathematically equivalent to the original QL algorithm, but the computation does not require ks1 to be actually subtracted from A. The algorithm is based on the following lemma: If A is a symmetric nonsingular matrix and B = QT · A · Q, where Q is orthogonal and B is tridiagonal with positive off-diagonal elements, then Q and B are fully determined when the last row of QT is specified. Proof: Let qT i denote the ith row vector of the matrix QT . Then qi is the ith column vector of the
11.3 Eigenvalues and Eigenvectors of a Tridiagonal Matrix 479 matrix Q.The relation B.Q=QT.A can be written 1Y1 9 a2 22 q2 A (11.3.12) an-1 Bn-1 n-1 The nth row of this matrix equation is anq7-1+Bnq7=q7·A (11.3.13) 三 81 Since Q is orthogonal, q·qm=dnm (11.3.14) Thus if we postmultiply equation(11.3.13)by q,we find 4州 Bn=qA·qn (11.3.15) which is known since qr is known.Then equation (11.3.13)gives anq7-1=z不-1 (11.3.16) where z-1=q·A-3nqd (11.3.17) is known.Therefore a品=z7-1Zn-1, (11.3.18) or 9 an Zn-1 (11.3.19) and gn-1 =t-1/an (11.3.20) CIENTIFIC (wheren is nonzero by hypothesis).Similarly,one can show by induction that if we know 4n4n-1...,qn-and the a's,B's,and y's up to level n-j,one can determine the 6 quantities at level n-(j+1). To apply the lemma in practice,suppose one can somehow find a tridiagonal matrix A+1 such that 五+1=可.五0 (11.3.21) where O is orthogonal and has the same last row as Q?in the original L algorithm Then QQ.and As+1 =A.+1. Numerical 10621 Now,in the original algorithm,from equation(11.3.11)we see that the last row ofQ? 43106 is the same as the last row of P But recall that P is a plane rotation designed to Recipes annihilate the (n-1,n)element of A.-1.A simple calculation using the expression (11.1.1)shows that it has parameters (outside Software. dn-ks -en-1 North e=√层+d-k7 8= (11.3.22) Ve品+(dn-ks The matrixA.is tridiagonal with 2extra elements: + (11.3.23) We must now reduce this to tridiagonal form with an orthogonal matrix whose last row is 00.]so that the last row ofwill stay equal topThis can be done by
11.3 Eigenvalues and Eigenvectors of a Tridiagonal Matrix 479 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). matrix Q. The relation B · QT = QT · A can be written β1 γ1 α2 β2 γ2 . . . αn−1 βn−1 γn−1 αn βn · qT 1 qT 2 . . . qT n−1 qT n = qT 1 qT 2 . . . 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) Since Q is 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 qn is known. Then equation (11.3.13) gives αnqT n−1 = z T n−1 (11.3.16) where z T n−1 ≡ qT n · A − βnqT n (11.3.17) is known. Therefore α2 n = z T n−1zn−1, (11.3.18) or αn = |zn−1| (11.3.19) and qT n−1 = z T n−1/αn (11.3.20) (where αn is nonzero by hypothesis). Similarly, one can show by induction that if we know qn, qn−1,..., qn−j and 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+1 such that As+1 = QT s · As · Qs (11.3.21) where QT s is orthogonal and has the same last row as QT s in the original QL algorithm. Then Qs = Qs and As+1 = As+1. 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−1 is 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 e2 n + (dn − ks)2 , s = −en−1 e2 n + (dn − ks)2 (11.3.22) The matrix P(s) n−1 · As · P(s)T n−1 is 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 s will stay equal to P(s) n−1. This can be done by
480 Chapter 11.Eigensystems 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 Q==.…2p21 (11.3.24) where the P's are the Givens rotations and Pn-is the same plane rotation as in the original algorithm.Then equation (11.3.21)gives the next iterate of A.Note that the shift ks enters implicitly through the parameters (11.3.22). 5常 81 The following routine tqli("Tridiagonal QL Implicit"),based algorithmically on the implementations in [2,31,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 later eigenvalues are liberated with very little work.The average number of iterations per eigenvalue is typically 1.3-1.6.The operation count per iteration is O(n), RECIPES I with a fairly large effective coefficient,say,~20n.The total operation count for the diagonalization is then 20n x(1.3-1.6)n ~30n2.If the eigenvectors are required,the statements indicated by comments are included and there is an Press. additional,much larger,workload of about 3n3 operations. #include #include "nrutil.h" SCIENTIFIC void tqli(f1oatd0,float e☐，intn,f1oat**z) QL algorithm with implicit shifts,to determine the eigenvalues and eigenvectors of a real,sym- metric,tridiagonal matrix,or of a real,symmetric matrix previously reduced by tred2 $11.2.On input,d[1..n]contains the diagonal elements of the tridiagonal matrix.On output,it returns the eigenvalues.The vector e[1..n]inputs the subdiagonal elements of the tridiagonal matrix, with e [1]arbitrary.On output e is destroyed.When finding only the eigenvalues,several lines may be omitted,as noted in the comments.If the eigenvectors of a tridiagonal matrix are de- COMPUTING (ISBN 188810920 sired,the matrix z[1..n][1..n]is input as the identity matrix.If the eigenvectors of a matrix that has been reduced by tred2 are required,then z is input as the matrix output by tred2. In either case,the kth column of z returns the normalized eigenvector corresponding to d[k] Numerica 10.621 float pythag(float a,float b); uctio Recipes 43108 int m,l,iter,i,k; float s,r,p,g,f,dd,c,b; for(i=2;1<n;1++)e[i-1]=e[i]; Convenient to renumber the el- (outside e[n]=0.0; ements of e. Software. for(1=1;1<=n;1++){ North iter=0; do for(m1;m<=n-1;m++) Look for a single small subdi- dd=fabs(d[m])+fabs(d[m+1]); agonal element to split if ((float)(fabs(e[m])+dd)==dd)break; the matrix. 2 1f(m1=1)[ if (iter++==30)nrerror("Too many iterations in tqli"); g=(d[1+1]-d[1])/(2.0*e[1]); Form shift. r=pythag(g,1.0); g-d[m]-d[l]+e[1]/(g+SIGN(r,g)); This is dm-ks s=C=1.0; p=0.0;
480 Chapter 11. Eigensystems Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). 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−1 is the same plane rotation as in the original algorithm. Then equation (11.3.21) gives the next iterate of A. Note that the shift ks enters implicitly through the parameters (11.3.22). The following routine tqli (“Tridiagonal QL Implicit”), based algorithmically 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 later eigenvalues are liberated with very little work. The average number of iterations 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 3n3 operations. #include #include "nrutil.h" void tqli(float d[], float e[], int n, float **z) 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. On input, d[1..n] contains the diagonal elements of the tridiagonal matrix. On output, it returns the eigenvalues. The vector e[1..n] inputs the subdiagonal elements of the tridiagonal matrix, with e[1] arbitrary. On output e is 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[1..n][1..n] is input as the identity matrix. If the eigenvectors of a matrix that has been reduced by tred2 are required, then z is input as the matrix output by tred2. In either case, the kth column of z returns the normalized eigenvector corresponding to d[k]. { float pythag(float a, float b); int m,l,iter,i,k; float s,r,p,g,f,dd,c,b; for (i=2;i<=n;i++) e[i-1]=e[i]; Convenient to renumber the ele[n]=0.0; ements of e. for (l=1;l<=n;l++) { iter=0; do { for (m=l;m<=n-1;m++) { Look for a single small subdiagonal element to split the matrix. dd=fabs(d[m])+fabs(d[m+1]); if ((float)(fabs(e[m])+dd) == dd) break; } if (m != l) { if (iter++ == 30) nrerror("Too many iterations in tqli"); g=(d[l+1]-d[l])/(2.0*e[l]); Form shift. r=pythag(g,1.0); g=d[m]-d[l]+e[l]/(g+SIGN(r,g)); This is dm − ks. s=c=1.0; p=0.0;
11.4 Hermitian Matrices 481 for(i=m-1;1>=1;1--）[ A plane rotation as in the origi- f=s*e[i]; nal OL.followed by Givens b=c*e[i]; rotations to restore tridiag- e[i+1]=(r=pythag(f,g)); onal form. 1f(r==0.0){ Recover from underflow. d[i+1]-=p; e[m]=0.0; break; s-f/r; c=g/r; gd[1+1]-p; r=(d[i]-g)*s+2.0*c*b; d[i+1]=g+(p=s*r）; g=c*r-b; 83g /Next loop can be omitted if eigenvectors not wanted*/ 19881992 for(k=1;k=1)continue; d[1]-=p; e[1]=g; e[m]=0.0; (North America 州bMe se make one paper University Press. THE while (m !1); ART 22 ictly proh Programs CITED REFERENCES AND FURTHER READING: Acton,F.S.1970,Numerica/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.ll of Handbook for Automatic Com- OF SCIENTIFIC COMPUTING(ISBN putation (New York:Springer-Verlag).[2] Smith,B.T..et al.1976,Matrix Eigensystem Routines-EISPACK Guide,2nd ed.,vol.6 of 1988-1992 Lecture Notes in Computer Science (New York:Springer-Verlag).[3] Stoer,J.,and Bulirsch,R.1980,Introduction to Numerical Analysis(New York:Springer-Verlag). Fuurggoglrion Numerica 10621 56.6.6.[4 43106 11.4 Hermitian Matrices (outside Recipes North Software. The complex analog of a real,symmetric matrix is a Hermitian matrix, satisfying equation(11.0.4).Jacobi transformations can be used to find eigenvalues and eigenvectors,as also can Householder reduction to tridiagonal form followed by visit website machine OL iteration.Complex versions of the previous routines jacobi,tred2,and tqli are quite analogous to their real counterparts.For working routines,consult [1.21. An alternative,using the routines in this book,is to convert the Hermitian problem to a real,symmetric one:If C=A+iB is a Hermitian matrix,then the n x n complex eigenvalue problem (A+iB)·(u+iv)=λ(u+iv) (11.4.1)
11.4 Hermitian Matrices 481 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). for (i=m-1;i>=l;i--) { A plane rotation as in the original QL, followed by Givens rotations to restore tridiagonal form. f=s*e[i]; b=c*e[i]; e[i+1]=(r=pythag(f,g)); if (r == 0.0) { Recover from underflow. d[i+1] -= p; e[m]=0.0; break; } s=f/r; c=g/r; g=d[i+1]-p; r=(d[i]-g)*s+2.0*c*b; d[i+1]=g+(p=s*r); g=c*r-b; /* Next loop can be omitted if eigenvectors not wanted*/ for (k=1;k= l) continue; d[l] -= p; e[l]=g; e[m]=0.0; } } while (m != l); } } CITED REFERENCES AND FURTHER READING: Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), pp. 331–335. [1] Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol. II of Handbook for Automatic Computation (New York: Springer-Verlag). [2] 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., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: 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 transformations can be used to find eigenvalues and eigenvectors, as also can Householder reduction to tridiagonal form followed by QL iteration. Complex versions of the previous routines jacobi, tred2, and tqli are quite analogous to their real counterparts. For working routines, 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 + iB is a Hermitian matrix, then the n × n complex eigenvalue problem (A + iB) · (u + iv) = λ(u + iv) (11.4.1)