f11-4

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] 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