f11-6

480 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).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] Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), §6.5.4. [3] 11.6 The QR Algorithm for Real Hessenberg Matrices Recall the following relations for the QRalgorithm with shifts: Qs·(As−ks1)=Rs (11.6.1 ) whereQis orthogonal and Ris upper triangular, and As+1=Rs·QT s+ks1 =Qs·As·QT s(11.6.2 ) TheQRtransformationpreservesthe upper Hessenbergform of the original matrix A≡A1, and the workload on such a matrix is O(n2)per iteration as opposed toO(n3)on a general matrix. As s→∞,Asconverges to a form where the eigenvaluesareeitherisolatedonthediagonalorareeigenvaluesofa 2×2submatrix on the diagonal. As we pointed out in §11.3, shifting is essential for rapid convergence. A key differencehereisthatanonsymmetricrealmatrixcanhavecomplexeigenvalues. This means that good choices for the shifts ksmay be complex,apparently necessitating complex arithmetic. Complex arithmetic can be avoided, however, by a clever trick. The trick dependsonaresultanalogoustothelemmaweusedforimplicitshiftsin §11.3. The lemma we need here states that if Bis a nonsingular matrix such that B·Q=Q·H (11.6.3 ) whereQisorthogonaland HisupperHessenberg,then QandHarefullydetermined bythefirstcolumnof Q. (Thedeterminationisuniqueif Hhaspositivesubdiagonal elements.) The lemma can be proved by induction analogously to the proof given for tridiagonal matrices in §11.3. The lemma is used in practice by taking two steps of the QRalgorithm, either with two real shifts ksandks+1, or with complex conjugate values ksand ks+1 =ks*. This gives a real matrix As+2, where As+2=Qs+1·Qs·As·QT s·QT s+1· (11.6.4 ) 11.6TheQR AlgorithmforRealHessenbergMatrices 481Sample 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).TheQ’s are determined by As−ks1=QT s·Rs (11.6.5 ) As+1=Qs·As·QT s (11.6.6 ) As+1−ks+11=QT s+1·Rs+1 (11.6.7 ) Using (11.6.6), equation (11.6.7) can be rewritten As−ks+11=QT s·QT s+1·Rs+1·Qs (11.6.8 ) Hence, if we define M=(As−ks+11)·(As−ks1)( 11.6.9 ) equations (11.6.5) and (11.6.8) give R=Q·M (11.6.10 ) where Q=Qs+1·Qs (11.6.11 ) R=Rs+1·Rs (11.6.12 ) Equation (11.6.4) can be rewritten As·QT=QT·As+2 (11.6.13 ) Thussupposewe cansomehowfindanupperHessenbergmatrix Hsuch that As·QT=QT·H (11.6.14 ) whereQis orthogonal. If QThas the same first columnas QT(i.e.,Qhas the same first row as Q), thenQ=QandAs+2 =H. The first row of Qis found as follows. Equation (11.6.10) shows that Qis the orthogonal matrix that triangularizes the real matrix M. Any real matrix can be triangularized by premultiplying it by a sequence of Householder matrices P1 (acting on the first column), P2(acting on the second column), ...,Pn−1. Thus Q=Pn−1···P2·P1, and the first row of Qis the first row of P1sincePiis an (i−1)×(i−1)identity matrix in the top left-hand corner. We now must find Q satisfying (11.6.14) whose first row is that of P1. The Householder matrix P1is determined by the first column of M. SinceAs is upper Hessenberg, equation (11.6.9) shows that the first column of Mhas the form [p1,q1,r1,0, ..., 0]T, where p1=a2 11−a11(ks+ks+1)+ksks+1+a12a21 q1=a21(a11+a22−ks−ks+1) r1=a21a32(11.6.15 ) 482 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).Hence P1=1−2w1·wT 1 (11.6.16 ) wherew1has only its first 3 elements nonzero (cf. equation 11.2.5). The matrix P1·As·PT 1is therefore upper Hessenberg with 3 extra elements: P1·A1·PT 1= ××××××× ××××××× x×××××× xx ××××× ×××× ××× ×× (11.6.17 ) ThismatrixcanberestoredtoupperHessenbergformwithoutaffectingthefirst row byasequenceofHouseholdersimilaritytransformations. ThefirstsuchHouseholder matrix,P 2, acts on elements 2, 3, and 4 in the first column, annihilating elements 3 and 4. This produces a matrix of the same form as (11.6.17), with the 3 extra elements appearing one column over:  ××××××× ××××××× ×××××× x××××× xx ×××× ××× ×× (11.6.18 ) Proceeding in this way up to P n−1, we see that at each stage the Householder matrixPrhas a vector wrthat is nonzero only in elements r,r+1, and r+2. These elements are determined by the elements r,r+1, and r+2in the (r−1)st column of the current matrix. Note that the preliminary matrix P1has the same structure as P2,...,Pn−1. The result is that Pn−1···P2·P1·As·PT 1·PT 2···PT n−1=H (11.6.19 ) whereHis upper Hessenberg. Thus Q=Q=Pn−1···P2·P1 (11.6.20 ) and As+2=H (11.6.21 ) The shifts of origin at each stage are taken to be the eigenvalues of the 2×2 matrix in the bottom right-hand corner of the current As. This gives ks+ks+2=an−1,n−1+ann ksks+1=an−1,n−1ann−an−1,nan,n−1(11.6.22 ) 11.6TheQR AlgorithmforRealHessenbergMatrices 483Sample 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).Substituting (11.6.22) in (11.6.15), we get p1=a21{[(ann−a11)(an−1,n−1−a11)−an−1,nan,n−1]/a21+a12} q1=a21[a22−a11−(ann−a11)−(an−1,n−1−a11)] r1=a21a32 (11.6.23 ) We have judiciously grouped terms to reduce possible roundoff when there are small off-diagonal elements. Since only the ratios of elements are relevant for a Householder transformation, we can omit the factor a21from (11.6.23). In summary, to carry out a double QRstep we construct the Householder matricesPr,r =1,...,n −1.F o rP1weuse p1,q1,and r1givenby(11.6.23). For theremainingmatrices, pr,qr,and rraredeterminedbythe (r, r−1),(r+1,r−1), and (r+2,r−1)elements of the current matrix. The number of arithmetic operations can be reduced by writing the nonzero elements of the 2w·wTpart of the Householder matrix in the form 2w·wT= (p±s)/(±s) q/(±s) r/(±s) ·[1 q/(p±s)r/(p±s)] ( 11.6.24 ) where s2=p2+q2+r2(11.6.25 ) (We have simply divided each element by a piece of the normalizing factor; cf. the equations in §11.2.) If we proceed in this way, convergence is usually very fast. There are two possiblewaysofterminatingtheiterationforaneigenvalue. First,if an,n−1becomes “negligible,”then annis an eigenvalue. We can thendelete the nth rowandcolumn ofthematrixandlookforthenexteigenvalue. Alternatively, an−1,n−2maybecome negligible. In this case the eigenvalues of the 2×2matrix in the lower right-hand corner may be taken to be eigenvalues. We delete the nth and (n−1)st rows and columns of the matrix and continue. Thetest forconvergencetoaneigenvalueiscombinedwithatest fornegligible subdiagonal elements that allows splitting of the matrix into submatrices. We findthe largest isuch that a i,i−1is negligible. If i=n, we have found a single eigenvalue. If i=n−1, we have found two eigenvalues. Otherwise we continue the iteration on the submatrix in rows iton(ibeing set to unity if there is no small subdiagonal element). Afterdetermining i,thesubmatrixinrows itonisexaminedtoseeifthe product ofanytwoconsecutivesubdiagonalelementsis smallenoughthatwecanworkwith an even smaller submatrix, starting say in row m. We start with m=n−2 and decrement it down to i+1, computing p,q, and raccording to equations (11.6.23)with 1 replaced by mand 2 by m+1. If these were indeed the elements of the special “first” Householder matrix in a double QRstep, then applying the Householder matrix would lead to nonzero elements in positions (m+1,m−1), (m+2,m−1),and (m+2,m). We requirethatthefirst twooftheseelementsbe 484 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).small compared with the local diagonal elements am−1,m−1,ammandam+1 ,m+1. A satisfactory approximate criterion is |am,m−1|(|q|+|r|)/lessmuch|p|(|am+1 ,m+1|+|amm|+|am−1,m−1|)(11.6.26 ) Very rarely, the procedure described so far will fail to converge. On such matrices, experience shows that if one double step is performed with any shifts that are of order the norm of the matrix, convergence is subsequently very rapid. Accordingly, if ten iterations occur without determining an eigenvalue, the usual shifts are replaced for the next iteration by shifts defined by ks+ks+1=1.5×(|an,n−1|+|an−1,n−2|) ksks+1=(|an,n−1|+|an−1,n−2|)2(11.6.27 ) The factor 1.5 was arbitrarily chosen to lessen the likelihood of an “unfortunate” choice of shifts. This strategy is repeated after 20 unsuccessful iterations. After 30unsuccessful iterations, the routine reports failure. Theoperationcountforthe QRalgorithmdescribedhereis ∼5k 2periteration, where kisthecurrentsizeofthematrix. Thetypicalaveragenumberofiterationsper eigenvalueis ∼1.8,sothetotaloperationcountforalltheeigenvaluesis ∼3n3. This estimateneglectsanypossibleefficiencyduetosplittingorsparsenessofthematrix. The following routine hqris based algorithmically on the above description, in turn following the implementations in [1,2]. SUBROUTINE hqr(a,n,np,wr,wi) INTEGER n,npREAL a(np,np),wi(np),wr(np) Finds all eigenvalues of an nbynupper Hessenberg matrix athat is stored in an npbynp array. On input acan be exactly as output from elmhes §11.5; on output it is destroyed. The real and imaginary parts of the eigenvalues are returned in wrandwi, respectively. INTEGER i,its,j,k,l,m,nn REAL anorm,p,q,r,s,t,u,v,w,x,y,zanorm=0. Compute matrix norm for possible use in locating single small subdiagonal element.do 12i=1,n do11j=max(i-1,1),n anorm=anorm+abs(a(i,j)) enddo 11 enddo 12 nn=nt=0. Gets changed only by an exceptional shift. 1 if(nn.ge.1)then Begin search for next eigenvalue. its=0 2d o 13l=nn,2,-1 Begin iteration: look for single small sub- diagonal element. s=abs(a(l-1,l-1))+abs(a(l,l)) if(s.eq.0.)s=anorm if(abs(a(l,l-1))+s.eq.s)then abs(a(l,l-1)=0.goto 3 endif enddo 13 l=1 3 x=a(nn,nn) if(l.eq.nn)then One root found. wr(nn)=x+twi(nn)=0. nn=nn-1 11.6TheQR AlgorithmforRealHessenbergMatrices 485Sample 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).else y=a(nn-1,nn-1) w=a(nn,nn-1)*a(nn-1,nn) if(l.eq.nn-1)then Two roots found... p=0.5*(y-x) q=p**2+w z=sqrt(abs(q))x=x+tif(q.ge.0.)then ...a real pair. z=p+sign(z,p) wr(nn)=x+zwr(nn-1)=wr(nn)if(z.ne.0.)wr(nn)=x-w/z wi(nn)=0. wi(nn-1)=0. else ...a complex pair. wr(nn)=x+p wr(nn-1)=wr(nn)wi(nn)=zwi(nn-1)=-z endif nn=nn-2 else No roots found. Continue iteration. if(its.eq.30)pause ’too many iterations in hqr’ if(its.eq.10.or.its.eq.20)then Form exceptional shift. t=t+xdo 14i=1,nn a(i,i)=a(i,i)-x enddo 14 s=abs(a(nn,nn-1))+abs(a(nn-1,nn-2)) x=0.75*s y=xw=-0.4375*s**2 endif its=its+1 do 15m=nn-2,l,-1 Form shift and then look for 2 consecu- tive small subdiagonal elements. z=a(m,m) r=x-z s=y-z p=(r*s-w)/a(m+1,m)+a(m,m+1) Equation (11.6.23). q=a(m+1,m+1)-z-r-s r=a(m+2,m+1) s=abs(p)+abs(q)+abs(r) Scale to prevent overflow or underflow. p=p/sq=q/s r=r/s if(m.eq.l)goto 4u=abs(a(m,m-1))*(abs(q)+abs(r)) v=abs(p)*(abs(a(m-1,m-1))+abs(z)+abs(a(m+1,m+1))) if(u+v.eq.v)goto 4 Equation (11.6.26). enddo 15 4d o 16i=m+2,nn a(i,i-2)=0. if (i.ne.m+2) a(i,i-3)=0. enddo 16 do19k=m,nn-1 D o u b l eQ Rs t e po nr o w s ltonnand columns mtonn. if(k.ne.m)then p=a(k,k-1) Begin setup of Householder vector. q=a(k+1,k-1) r=0. if(k.ne.nn-1)r=a(k+2,k-1)x=abs(p)+abs(q)+abs(r)if(x.ne.0.)then p=p/x Scale to prevent overflow or underflow. 486 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).q=q/x r=r/x endif endifs=sign(sqrt(p**2+q**2+r**2),p) if(s.ne.0.)then if(k.eq.m)then if(l.ne.m)a(k,k-1)=-a(k,k-1) else a(k,k-1)=-s*x endifp=p+s Equations (11.6.24). x=p/s y=q/s z=r/sq=q/p r=r/p do 17j=k,nn Row modification. p=a(k,j)+q*a(k+1,j)if(k.ne.nn-1)then p=p+r*a(k+2,j) a(k+2,j)=a(k+2,j)-p*z endif a(k+1,j)=a(k+1,j)-p*y a(k,j)=a(k,j)-p*x enddo 17 do18i=l,min(nn,k+3) Column modification. p=x*a(i,k)+y*a(i,k+1) if(k.ne.nn-1)then p=p+z*a(i,k+2) a(i,k+2)=a(i,k+2)-p*r endifa(i,k+1)=a(i,k+1)-p*qa(i,k)=a(i,k)-p enddo 18 endif enddo 19 goto 2 ...for next iteration on current eigenvalue. endif endif goto 1 ...for next eigenvalue. endif returnEND 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] Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press),§7.5. 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] 11.7EigenvaluesorEigenvectorsbyInverseIteration 487Sample 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).11.7 Improving Eigenvalues and/or Finding Eigenvectors by Inverse Iteration The basic idea behind inverse iteration is quite simple. Let ybe the solution of the linear system (A−τ1)·y=b (11.7.1 ) wherebis a random vector and τis close to some eigenvalue λofA. Then the solutionywill be close to the eigenvector corresponding to λ. The procedure can be iterated: Replace bbyyand solve for a new y, which will be even closer to the true eigenvector. We can see why this works by expanding both yandbas linear combinations of the eigenvectors xjofA: y=/summationdisplay jαjxjb=/summationdisplay jβjxj (11.7.2 ) Then (11.7.1) gives /summationdisplay jαj(λj−τ)xj=/summationdisplay jβjxj (11.7.3 ) so that αj=βj λj−τ(11.7.4 ) and y=/summationdisplay jβjxj λj−τ(11.7.5 ) Ifτis close to λn, say, then provided βnis not accidentally too small, ywill be approximately xn, upto a normalization. Moreover,eachiterationofthis procedure givesanotherpowerof λj−τinthedenominatorof(11.7.5). Thustheconvergence is rapid for well-separated eigenvalues. Suppose at the kth stage of iteration we are solving the equation (A−τk1)·y=bk (11.7.6 ) wherebkandτkare our current guesses for some eigenvector and eigenvalue of interest (let’s say, xnandλn). Normalize bkso thatbk·bk=1. The exact eigenvector and eigenvalue satisfy A·xn=λnxn (11.7.7 ) so (A−τk1)·xn=(λn−τk)xn (11.7.8 ) Sinceyof (11.7.6)is an improvedapproximationto xn, we normalizeit and set bk+1=y |y|(11.7.9 )