f11-2

462 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).SUBROUTINE eigsrt(d,v,n,np) INTEGER n,np REAL d(np),v(np,np) Giventheeigenvalues dandeigenvectors vasoutput from jacobi(§11.1)or tqli(§11.3), this routine sorts the eigenvalues into descending order, and rearranges the columns of v correspondingly. The method is straight insertion. INTEGER i,j,kREAL pdo 13i=1,n-1 k=i p=d(i)do 11j=i+1,n if(d(j).ge.p)then k=j p=d(j) endif enddo 11 if(k.ne.i)then d(k)=d(i)d(i)=p do 12j=1,n p=v(j,i)v(j,i)=v(j,k) v(j,k)=p enddo 12 endif enddo 13 returnEND CITED REFERENCES AND FURTHER READING: Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press),§8.4. 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). [1] Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra , vol. II of Handbook for Automatic Com- putation(New York: Springer-Verlag). [2] 11.2 Reduction of a Symmetric Matrix to Tridiagonal Form: Givens and Householder Reductions As already mentioned, the optimum strategy for finding eigenvalues and eigenvectors is, first, to reduce the matrix to a simple form, only then beginning an iterativeprocedure. Forsymmetricmatrices,thepreferredsimpleformistridiagonal. TheGivens reduction is a modification of the Jacobi method. Instead of trying to reduce the matrix all the way to diagonal form, we are content to stop when thematrix is tridiagonal. This allows the procedureto be carried out in a finite number of steps, unlike the Jacobi method, which requires iteration to convergence. 11.2ReductionofaSymmetric MatrixtoTridiagonalForm 463Sample 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).Givens Method For the Givens method, we choose the rotation angle in equation (11.1.1) so as to zero an element that is notat one of the four “corners,” i.e., not app,apq, oraqqin equation (11.1.3). Specifically, we first choose P23to annihilate a31 (and, by symmetry, a13). Then we choose P24to annihilate a41. In general, we choose the sequence P23,P24,...,P2n;P34,...,P3n;... ;Pn−1,n wherePjkannihilates ak,j−1. The methodworks becauseelements such as a/prime rpand a/prime rq, with r/negationslash=pr/negationslash=q, are linear combinationsof the old quantities arpandarq,b y equation (11.1.4). Thus, if arpandarqhave already been set to zero, they remain zero as the reduction proceeds. Evidently, of order n2/2rotations are required, and the number of multiplications in a straightforward implementation is of order 4n3/3, not counting those for keeping track of the product of the transformation matrices, required for the eigenvectors. The Householder method, to be discussed next, is just as stable as the Givens reductionandit is afactorof2moreefficient,sotheGivensmethodis notgenerallyused. Recentwork(see [1])hasshownthattheGivensreductioncanbereformulated to reduce the number of operations by a factor of 2, and also avoid the necessity of taking square roots. This appears to make the algorithm competitive with theHouseholder reduction. However, this “fast Givens” reduction has to be monitored to avoid overflows, and the variables have to be periodically rescaled. There does not seem to be any compelling reason to prefer the Givens reduction over the Householder method. Householder Method TheHouseholderalgorithmreducesan n×nsymmetricmatrix Atotridiagonal form by n−2orthogonal transformations. Each transformation annihilates the requiredpartofawholecolumnandwholecorrespondingrow. Thebasicingredient is a Householder matrix P, which has the form P=1−2w·wT(11.2.1 ) wherewis a real vectorwith |w|2=1. (In the present notation, the outeror matrix productoftwovectors, aandbiswritten a·bT,whilethe innerorscalarproductof the vectors is written as aT·b.) The matrix Pis orthogonal, because P2=(1−2w·wT)·(1−2w·wT) =1−4w·wT+4w·(wT·w)·wT =1(11.2.2 ) Therefore P=P−1. ButPT=P, and soPT=P−1, provingorthogonality. 464 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).RewritePas P=1−u·uT H(11.2.3 ) wherethescalar His H≡1 2|u|2(11.2.4 ) anducan now be any vector. Suppose xis the vector composed of the first column ofA. Choose u=x∓|x|e1 (11.2.5 ) wheree1is the unit vector [1,0,..., 0]T, and the choice of signs will be made later. Then P·x=x−u H·(x∓|x|e1)T·x =x−2u·(|x|2∓|x|x1) 2|x|2∓2|x|x1 =x−u =±|x|e1(11.2.6 ) This shows that the Householder matrix Pacts on a given vector xto zero all its elements except the first one. To reduce a symmetric matrix Ato tridiagonal form, we choose the vector x for the first Householder matrix to be the lower n−1elements of the first column. Then the lower n−2elements will be zeroed: P1·A= 1 00 ··· 0 0 0 ... (n−1)P1 0 · a 11 a12a13 ··· a1n a21 a31 ...irrelevant an1  = a 11 a12a13 ··· a1n k 0 ...irrelevant 0 (11.2.7 ) Here we have written the matrices in partitioned form, with (n−1)Pdenoting a Householder matrix with dimensions (n−1)×(n−1). The quantity kis simply plus or minus the magnitude of the vector [a21,...,a n1]T. 11.2ReductionofaSymmetric MatrixtoTridiagonalForm 465Sample 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).The complete orthogonal transformation is now A/prime=P·A·P= a 11 k 0 ··· 0 k 0 ... irrelevant 0 (11.2.8 ) We have used the fact that P T=P. Now choose the vector xfor the second Householder matrix to be the bottom n−2elements of the second column, and from it construct P2≡ 10 0 ··· 0 01 0 ··· 0 00 ......(n−2)P2 00 (11.2.9 ) Theidentityblockintheupperleftcornerinsuresthatthetridiagonalizationachieved in the first step will not be spoiled by this one, while the (n−2)-dimensional Householdermatrix (n−2)P2createsoneadditionalrowandcolumnofthetridiagonal output. Clearly, a sequence of n−2such transformations will reduce the matrix Ato tridiagonal form. Instead of actually carrying out the matrix multiplications in P·A·P,w e compute a vector p≡A·u H(11.2.10 ) Then A·P=A·(1−u·uT H)=A−p·uT A/prime=P·A·P=A−p·uT−u·pT+2Ku·uT where the scalar Kis defined by K=uT·p 2H(11.2.11 ) Ifwe write q≡p−Ku (11.2.12 ) thenwehave A/prime=A−q·uT−u·qT(11.2.13 ) This is the computationally useful formula. Following [2],theroutineforHouseholderreductiongivenbelowactuallystarts in the nth column of A, not the first as in the explanation above. In detail, the equationsareasfollows: Atstage m(m=1,2,...,n −2)thevector uhastheform uT=[ai1,a i2,...,a i,i−2,a i,i−1±√σ,0,..., 0] ( 11.2.14 ) 466 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).Here i≡n−m+1= n, n−1,..., 3( 11.2.15 ) and the quantity σ(|x|2in our earlier notation) is σ=(ai1)2+··· +(ai,i−1)2(11.2.16 ) We choose the sign of σin (11.2.14) to be the same as the sign of ai,i−1to lessen roundoff error. Variablesarethuscomputedinthefollowingorder: σ,u,H ,p,K ,q,A/prime.A ta n y stage m,Ais tridiagonal in its last m−1rows and columns. Iftheeigenvectorsofthefinaltridiagonalmatrixarefound(forexample,bythe routine in the next section), then the eigenvectorsof Acan be obtained by applying the accumulated transformation Q=P1·P2···Pn−2 (11.2.17 ) to those eigenvectors. We therefore form Qby recursion after all the P’s have been determined: Qn−2=Pn−2 Qj=Pj·Qj+1,j =n−3,..., 1 Q=Q1(11.2.18 ) The inputparametersforthe routinebelow are the n×nreal, symmetricmatrix a, stored in an np×nparray. On output, acontains the elements of the orthogonal matrix q. The vector dreturns the diagonal elements of the tridiagonal matrix A/prime, while the vector ereturns the off-diagonalelements in its components 2through n, with e(1)=0. Note that since ais overwritten,youshouldcopyit beforecallingthe routine, if it is required for subsequent computations. Noextrastoragearraysareneededfortheintermediateresults. Atstage m,the vectorspandqare nonzero only in elements 1,...,i(recall that i=n−m+1), whileuis nonzero only in elements 1,...,i −1. The elements of the vector eare being determined in the order n, n−1,..., so we can store pin the elements of e not already determined. The vector qcan overwrite poncepis no longer needed. We store uin the ith row of aandu/Hin the ith column of a. Once the reduction is complete, we compute the matrices Qjusing the quantities uandu/Hthat have been stored in a. SinceQjis an identity matrix in the last n−j+1rows and columns, we only need compute its elements up to row and column n−j. These canoverwritethe u’s andu/H’sinthecorrespondingrowsandcolumnsof a,which are no longer required for subsequent Q’s. Theroutine tred2,givenbelow,includesonefurtherrefinement. Ifthequantity σis zero or “small” at any stage, one can skip the corresponding transformation. A simple criterion, such as σ<smallest positivenumberrepresentableonmachine machineprecision 11.2ReductionofaSymmetric MatrixtoTridiagonalForm 467Sample 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).would be fine most of the time. A more careful criterion is actually used. Define the quantity /epsilon1=i−1/summationdisplay k=1|aik| (11.2.19 ) If/epsilon1=0to machineprecision,we skip the transformation. Otherwisewe redefine aikbecomes aik//epsilon1 (11.2.20 ) and use the scaled variables for the transformation. (A Householder transformation depends only on the ratios of the elements.) Notethatwhendealingwithamatrixwhoseelementsvaryovermanyordersof magnitude, it is important that the matrix be permuted, insofar as possible, so that the smaller elements are in the top left-hand corner. This is because the reductionis performedstarting fromthe bottomright-handcorner,and a mixtureof small and large elements there can lead to considerable rounding errors. Theroutine tred2isdesignedforusewiththeroutine tqliofthenextsection. tqlifinds the eigenvalues and eigenvectors of a symmetric, tridiagonal matrix. The combination of tred2andtqliis the most efficient known technique for finding all the eigenvalues and eigenvectors (or just all the eigenvalues) of a real, symmetric matrix. Inthelistingbelow,thestatementsindicatedbycommentsarerequiredonlyfor subsequentcomputationof eigenvectors. If onlyeigenvaluesare required,omission ofthecommentedstatementsspeedsuptheexecutiontimeof tred2byafactorof2 forlarge n. Inthelimitoflarge n,theoperationcountoftheHouseholderreduction is2n 3/3foreigenvaluesonly,and 4n3/3forboth eigenvaluesand eigenvectors. SUBROUTINE tred2(a,n,np,d,e) INTEGER n,np REAL a(np,np),d(np),e(np) Householder reduction of a real, symmetric, nbynmatrix a,s t o r e di na n npbynpphysical array. On output, ais replaced by the orthogonal matrix Qeffecting the transformation. d returns the diagonal elements of the tridiagonal matrix, and ethe off-diagonal elements, withe(1)=0. Severalstatements, asnoted incomments, canbeomittedifonlyeigenvalues are to be found, in which case acontains no useful information on output. Otherwise they are to be included. INTEGER i,j,k,l REAL f,g,h,hh,scaledo 18i=n,2,-1 l=i-1 h=0.scale=0.if(l.gt.1)then do 11k=1,l scale=scale+abs(a(i,k)) enddo 11 if(scale.eq.0.)then Skip transformation. e(i)=a(i,l) else do12k=1,l a(i,k)=a(i,k)/scale Use scaled a’s for transformation. h=h+a(i,k)**2 Form σinh. enddo 12 f=a(i,l) 468 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).g=-sign(sqrt(h),f) e(i)=scale*g h=h-f*g Now his equation (11.2.4). a(i,l)=f-g Storeuin the ith row of a. f=0. do15j=1,l C Omit following line if finding only eigenvalues a(j,i)=a(i,j)/h Storeu/Hinith column of a. g=0. Form an element of A·uing. do13k=1,j g=g+a(j,k)*a(i,k) enddo 13 do14k=j+1,l g=g+a(k,j)*a(i,k) enddo 14 e(j)=g/h Form element of pin temporarily unused f=f+e(j)*a(i,j) element of e. enddo 15 hh=f/(h+h) Form K, equation (11.2.11). do17j=1,l Formqand store in eoverwriting p. f=a(i,j) g=e(j)-hh*fe(j)=g do 16k=1,j Reduce a, equation (11.2.13). a(j,k)=a(j,k)-f*e(k)-g*a(i,k) enddo 16 enddo 17 endif else e(i)=a(i,l) endif d(i)=h enddo 18 C Omit following line if finding only eigenvalues. d(1)=0. e(1)=0.do 24i=1,n Beginaccumulation oftransformation matrices. C Delete lines from here ... l=i-1 if(d(i).ne.0.)then This block skipped when i=1. do22j=1,l g=0. do19k=1,l Useuandu/Hstored in ato formP·Q. g=g+a(i,k)*a(k,j) enddo 19 do21k=1,l a(k,j)=a(k,j)-g*a(k,i) enddo 21 enddo 22 endif C ... to here when finding only eigenvalues. d(i)=a(i,i) This statement remains. C Also delete lines from here ... a(i,i)=1. R e s e tr o wa n dc o l u m no f ato identity matrix for next iteration. do23j=1,l a(i,j)=0. a(j,i)=0. enddo 23 C ... to here when finding only eigenvalues. enddo 24 return END 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 )