f11-5

476 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).11.5 Reduction of a General Matrix to Hessenberg Form The algorithms for symmetric matrices, given in the preceding sections, are highly satisfactory in practice. By contrast, it is impossible to design equallysatisfactory algorithms for the nonsymmetric case. There are two reasons for this. First,theeigenvaluesofanonsymmetricmatrixcanbeverysensitivetosmallchanges in the matrix elements. Second, the matrix itself can be defective, so that there isno complete set of eigenvectors. We emphasize that these difficulties are intrinsic propertiesofcertainnonsymmetricmatrices,andnonumericalprocedurecan“cure” them. Thebestwecanhopeforareproceduresthatdon’texacerbatesuchproblems. Thepresenceofroundingerrorcan onlymake thesituation worse. With finite- precision arithmetic, one cannot even design a foolproof algorithm to determinewhether a given matrix is defective or not. Thus current algorithms generally tryto findacomplete set ofeigenvectors,andrelyontheuser toinspecttheresults. Ifany eigenvectors are almost parallel, the matrix is probably defective. Apartfromreferringyoutotheliterature,andtothecollectedroutinesin [1,2],we are goingto sidestep the problemofeigenvectors,givingalgorithmsforeigenvalues only. Ifyourequirejusta feweigenvectors,youcanread §11.7andconsiderfinding them by inverse iteration. We consider the problem of finding alleigenvectors of a nonsymmetric matrix as lying beyond the scope of this book. Balancing The sensitivity of eigenvalues to rounding errors during the execution of some algorithms can be reduced by the procedure of balancing . The errors in the eigensystem found by a numerical procedure are generally proportional to theEuclidean norm of the matrix, that is, to the square root of the sum of the squares of the elements. The idea of balancing is to use similarity transformations to make corresponding rows and columns of the matrix have comparable norms, thus reducing the overall norm of the matrix while leaving the eigenvalues unchanged. A symmetric matrix is already balanced. Balancing is a procedure with of order N 2operations. Thus, the time taken by the procedure balanc, given below, should never be more than a few percent of the total time required to find the eigenvalues. It is therefore recommended thatyoualwaysbalance nonsymmetric matrices. It never hurts, and it can substantially improvethe accuracyof the eigenvaluescomputedfor a badlybalanced matrix. TheactualalgorithmusedisduetoOsborne,asdiscussedin [1]. Itconsistsofa sequenceofsimilaritytransformationsbydiagonalmatrices D. Toavoidintroducing roundingerrors during the balancing process, the elements of Dare restricted to be exactpowersoftheradixbaseemployedforfloating-pointarithmetic(i.e.,2formost machines, but 16 for IBM mainframe architectures). The output is a matrix that is balanced in the norm given by summing the absolute magnitudes of the matrixelements. ThisismoreefficientthanusingtheEuclideannorm,andequallyeffective: A large reduction in one norm implies a large reduction in the other. Note that if the off-diagonal elements of any row or column of a matrix are all zero, then the diagonal element is an eigenvalue. If the eigenvalue happens to 11.5ReductionofaGeneralMatrixtoHessenbergForm 477Sample 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).be ill-conditioned (sensitive to small changes in the matrix elements), it will have relatively large errors when determined by the routine hqr(§11.6). Had we merely inspected the matrix beforehand, we could have determined the isolated eigenvalue exactly and then deleted the corresponding row and column from the matrix. You should consider whether such a pre-inspection might be useful in your application.(For symmetric matrices, the routines we gave will determine isolated eigenvalues accurately in all cases.) The routine balancdoes not keep track of the accumulated similarity trans- formation of the original matrix, since we will only be concerned with finding eigenvalues of nonsymmetric matrices, not eigenvectors. Consult [1-3]if you want to keep track of the transformation. SUBROUTINE balanc(a,n,np) INTEGER n,np REAL a(np,np),RADIX,SQRDXPARAMETER (RADIX=2.,SQRDX=RADIX**2) Given an nbynmatrix astored in an array of physical dimensions npbynp, this routine replaces it by a balanced matrix with identical eigenvalues. A symmetric matrix is already balancedandisunaffectedbythisprocedure. Theparameter RADIXshouldbethemachine’s floating-point radix. INTEGER i,j,last REAL c,f,g,r,s 1 continue last=1 do14i=1,n Calculate row and column norms. c=0.r=0. do 11j=1,n if(j.ne.i)then c=c+abs(a(j,i))r=r+abs(a(i,j)) endif enddo 11 if(c.ne.0..and.r.ne.0.)then If both are nonzero, g=r/RADIX f=1.s=c+r 2 if(c.lt.g)then find the integer power of the machine radix that comes closest to balancing the matrix. f=f*RADIX c=c*SQRDX goto 2endif g=r*RADIX 3 if(c.gt.g)then f=f/RADIX c=c/SQRDX goto 3endifif((c+r)/f.lt.0.95*s)then last=0 g=1./fdo 12j=1,n Apply similarity transformation. a(i,j)=a(i,j)*g enddo 12 do13j=1,n a(j,i)=a(j,i)*f enddo 13 endif endif enddo 14 478 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).if(last.eq.0)goto 1 return END Reduction to Hessenberg Form Thestrategyforfindingtheeigensystemofageneralmatrixparallelsthatofthe symmetric case. First we reducethe matrix to a simpler form, and then we perform aniterativeprocedureonthesimplifiedmatrix. The simplerstructurewe usehereis calledHessenberg form. An upper Hessenberg matrix has zeros everywhere below the diagonal except for the first subdiagonal row. For example, in the 6×6case, the nonzero elements are:  ×××××× ×××××× ××××× ×××× ××× ××  By now you should be able to tell at a glance that such a structure can be achieved by a sequence of Householder transformations, each one zeroing the required elements in a column of the matrix. Householder reduction to Hessenberg form is in fact an accepted technique. An alternative, however, is a procedureanalogous to Gaussian elimination with pivoting. We will use this elimination proceduresinceitis aboutafactorof2moreefficientthantheHouseholdermethod, and also since we want to teach you the method. It is possible to constructmatricesfor which the Householder reduction, being orthogonal, is stable and elimination is not, but such matrices are extremely rare in practice. Straight Gaussian elimination is not a similarity transformation of the matrix. Accordingly, the actual elimination procedure used is slightly different. Before the rth stage, the original matrix A≡A 1has become Ar, which is upper Hessenberg in its first r−1rows and columns. The rth stage then consists of the following sequence of operations: •Find the element of maximum magnitude in the rth column below the diagonal. If it is zero, skip the next two “bullets” and the stage is done. Otherwise, suppose the maximum element was in row r/prime. •Interchange rows r/primeand r+1. This is the pivoting procedure. To make the permutation a similarity transformation, also interchange columns r/prime and r+1. •For i=r+2 ,r+3 ,...,N, compute the multiplier ni,r+1≡air ar+1 ,r Subtract ni,r+1times row r+1from row i. To make the elimination a similaritytransformation,also add ni,r+1timescolumn itocolumn r+1. A total of N−2such stages are required. 11.5Reductionofa GeneralMatrix toHessenbergForm 479Sample 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).Whenthemagnitudesofthematrixelementsvaryovermanyorders,youshould trytorearrangethematrixsothatthelargestelementsareinthetopleft-handcorner.This reduces the roundofferror,since the reductionproceedsfrom left to right. Since we are concerned only with eigenvalues, the routine elmhesdoes not keep track of the accumulated similarity transformation. The operation count isabout 5N 3/6for large N. SUBROUTINE elmhes(a,n,np) INTEGER n,np REAL a(np,np) Reduction to Hessenberg form by the elimination method. The real, nonsymmetric, nby nmatrix a, stored in an array of physical dimensions npbynp, is replaced by an upper Hessenberg matrix with identical eigenvalues. Recommended, but not required, isthat this routine be preceded by balanc. On output, the Hessenberg matrix is in elements a(i,j) with i≤j+1. Elements with i>j+1are to be thought of as zero, but are returned with random values. INTEGER i,j,m REAL x,ydo 17m=2,n-1 m is called r+1in the text. x=0. i=mdo 11j=m,n Find the pivot. if(abs(a(j,m-1)).gt.abs(x))then x=a(j,m-1) i=j endif enddo 11 if(i.ne.m)then Interchange rows and columns. do12j=m-1,n y=a(i,j) a(i,j)=a(m,j) a(m,j)=y enddo 12 do13j=1,n y=a(j,i)a(j,i)=a(j,m)a(j,m)=y enddo 13 endif if(x.ne.0.)then Carry out the elimination. do16i=m+1,n y=a(i,m-1) if(y.ne.0.)then y=y/x a(i,m-1)=y do14j=m,n a(i,j)=a(i,j)-y*a(m,j) enddo 14 do15j=1,n a(j,m)=a(j,m)+y*a(j,i) enddo 15 endif enddo 16 endif enddo 17 returnEND 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 ) The QRtransformationpreservesthe 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 ksand ks+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 )