f11-0

Sample 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).Chapter 11. Eigensystems 11.0 Introduction AnN×NmatrixAis said to have an eigenvector xand corresponding eigenvalue λif A·x=λx (11.0.1 ) Obviously any multiple of an eigenvector xwill also be an eigenvector, but we won’tconsidersuchmultiplesasbeingdistincteigenvectors. (Thezerovectorisnotconsidered to be an eigenvector at all.) Evidently (11.0.1) can hold only if det|A−λ1|=0 ( 11.0.2 ) which,ifexpandedout,isan Nthdegreepolynomialin λwhoserootsaretheeigen- values. This proves that there are always N(not necessarily distinct) eigenvalues. Equaleigenvaluescomingfrommultiplerootsarecalled degenerate . Root-searching in the characteristic equation (11.0.2) is usually a very poor computational methodfor finding eigenvalues. We will learn much better ways in this chapter, as well as efficient ways for finding corresponding eigenvectors. The above two equations also prove that every one of the Neigenvalues has a (not necessarily distinct) corresponding eigenvector: If λis set to an eigenvalue, then the matrix A−λ1is singular, and we know that every singular matrix has at least onenonzerovectorinits nullspace(see §2.6onsingularvaluedecomposition). If you add τxto both sides of (11.0.1),you will easily see that the eigenvalues of any matrix can be changed or shiftedby an additive constant τby adding to the matrix that constant times the identity matrix. The eigenvectors are unchanged by this shift. Shifting, as we will see, is an important part of many algorithms for computing eigenvalues. We see also that there is no special significance to a zeroeigenvalue. Any eigenvalue can be shifted to zero, or any zero eigenvalue can be shifted away from zero. 449 450 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).Definitionsand Basic Facts A matrix is called symmetric if it is equal to its transpose, A=ATor aij=aji (11.0.3 ) Itiscalled Hermitian orself-adjoint ifitequalsthecomplex-conjugateofitstranspose (itsHermitian conjugate , denoted by “ †”) A=A†or aij=aji*( 11.0.4 ) It is termed orthogonal if its transpose equals its inverse, AT·A=A·AT=1 (11.0.5 ) andunitaryif its Hermitian conjugate equals its inverse. Finally, a matrix is called normalif itcommutes with its Hermitian conjugate, A·A†=A†·A (11.0.6 ) For real matrices, Hermitian means the same as symmetric, unitary means the same as orthogonal, and bothof these distinct classes are normal. Thereasonthat“Hermitian”isanimportantconcepthastodowitheigenvalues. The eigenvalues of a Hermitian matrix are all real. In particular, the eigenvalues of a real symmetric matrix are all real. Contrariwise, the eigenvalues of a realnonsymmetricmatrixmayincluderealvalues,butmayalsoincludepairsofcomplex conjugate values; and the eigenvalues of a complex matrix that is not Hermitian will in general be complex. The reason that “normal” is an important concept has to do with the eigen- vectors. The eigenvectors of a normal matrix with nondegenerate (i.e., distinct)eigenvaluesarecompleteandorthogonal,spanningthe N-dimensionalvectorspace. Foranormalmatrixwithdegenerateeigenvalues,wehavetheadditionalfreedomof replacing the eigenvectorscorrespondingto a degenerate eigenvalueby linear com-binationsofthemselves. Usingthisfreedom,wecanalwaysperformGram-Schmidt orthogonalization(consultanylinearalgebratext)and finda set ofeigenvectorsthat are complete and orthogonal, just as in the nondegenerate case. The matrix whose columns are an orthonormalset of eigenvectors is evidently unitary. A special case is that the matrix of eigenvectors of a real, symmetric matrix is orthogonal, sincethe eigenvectors of that matrix are all real. When a matrix is not normal, as typified by any random, nonsymmetric, real matrix,theningeneralwecannotfind anyorthonormalsetofeigenvectors,noreven anypairsofeigenvectorsthatareorthogonal(exceptperhapsbyrarechance). While theNnon-orthonormaleigenvectors will “usually” span the N-dimensional vector space, they do not always do so; that is, the eigenvectors are not always complete. Such a matrix is said to be defective. 11.0 Introduction 451Sample 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).Left and RightEigenvectors While the eigenvectors of a non-normal matrix are not particularly orthogonal among themselves, they dohave an orthogonality relation with a different set of vectors, which we must now define. Up to now our eigenvectorshave been column vectors that are multiplied to the right of a matrix A, as in (11.0.1). These, more explicitly, are termed right eigenvectors . We could also, however, try to find row vectors, which multiply Ato the left and satisfy x·A=λx (11.0.7 ) These are called left eigenvectors . By taking the transpose of equation (11.0.7), we seethateverylefteigenvectoris thetransposeofarighteigenvector ofthetranspose ofA. Now by comparing to (11.0.2), and using the fact that the determinant of a matrix equals the determinant of its transpose, we also see that the left and righteigenvaluesofAare identical. If the matrix Ais symmetric, then the left and right eigenvectors are just transposes of each other, that is, have the same numerical values as components. Likewise, if the matrix is self-adjoint, the left and right eigenvectors are Hermitian conjugates of each other. For the general nonnormal case, however, we have thefollowing calculation: Let X Rbe the matrix formed by columns from the right eigenvectors,and XLbethematrixformedbyrowsfromthelefteigenvectors. Then (11.0.1) and (11.0.7) can be rewritten as A·XR=XR·diag (λ1...λ N)XL·A=diag (λ1...λ N)·XL (11.0.8 ) Multiplying the first of these equations on the left by XL, the second on the right byXR, and subtracting the two, gives (XL·XR)·diag (λ1...λ N)=diag (λ1...λ N)·(XL·XR)( 11.0.9 ) Thissaysthatthematrixofdotproductsoftheleftandrighteigenvectorscommutes with the diagonalmatrix ofeigenvalues. But the only matricesthat commutewith adiagonalmatrix ofdistinctelements arethemselvesdiagonal. Thus,iftheeigenvalues arenondegenerate,eachlefteigenvectorisorthogonaltoallrighteigenvectorsexcept its correspondingone, and vice versa. By choice of normalization,the dot productsofcorrespondingleftandrighteigenvectorscanalwaysbemadeunityforanymatrix with nondegenerate eigenvalues. If some eigenvalues are degenerate, then either the left or the right eigenvec- tors corresponding to a degenerate eigenvalue must be linearly combined among themselves to achieve orthogonality with the right or left ones, respectively. Thiscan always be done by a procedure akin to Gram-Schmidt orthogonalization. The normalizationcanthenbeadjustedtogiveunityforthenonzerodotproductsbetween correspondingleftandrighteigenvectors. Ifthedotproductofcorrespondingleftandright eigenvectors is zero at this stage, then you have a case where the eigenvectors are incomplete! Note that incomplete eigenvectors can occur only where there are degenerate eigenvalues, but do not always occur in such cases (in fact, never occur for the class of “normal” matrices). See [1]for a clear discussion. 452 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).In both the degenerate and nondegenerate cases, the final normalization to unity of all nonzero dot products produces the result: The matrix whose rowsare left eigenvectors is the inverse matrix of the matrix whose columns are right eigenvectors, if the inverse exists . Diagonalizationofa Matrix Multiplying the first equation in (11.0.8) by XL, and using the fact that XL andXRare matrix inverses, we get X−1 R·A·XR=diag (λ1...λ N)( 11.0.10 ) This is a particular case of a similarity transform of the matrix A, A→Z−1·A·Z (11.0.11 ) for some transformation matrix Z. Similarity transformations play a crucial role in the computation of eigenvalues, because they leave the eigenvalues of a matrix unchanged. This is easily seen from det/vextendsingle/vextendsingleZ−1·A·Z−λ1/vextendsingle/vextendsingle=det/vextendsingle/vextendsingleZ−1·(A−λ1)·Z/vextendsingle/vextendsingle =det|Z|det|A−λ1|det/vextendsingle/vextendsingleZ−1/vextendsingle/vextendsingle =det|A−λ1|(11.0.12 ) Equation(11.0.10)showsthatanymatrixwithcompleteeigenvectors(whichincludes all normal matrices and “most” random nonnormal ones) can be diagonalized by a similarity transformation,that the columns of the transformationmatrix that effects the diagonalization are the right eigenvectors, and that the rows of its inverse are the left eigenvectors. For real, symmetricmatrices, the eigenvectorsare real and orthonormal,so the transformation matrix is orthogonal. The similarity transformation is then also an orthogonal transformation of the form A→ZT·A·Z (11.0.13 ) Whilerealnonsymmetricmatricescanbediagonalizedintheirusualcaseofcomplete eigenvectors,thetransformationmatrixisnotnecessarilyreal. Itturnsout,however, thatarealsimilaritytransformationcan“almost”dothejob. Itcanreducethematrixdown to a form with little two-by-twoblocks along the diagonal, all other elements zero. Each two-by-two block corresponds to a complex-conjugatepair of complex eigenvalues. Wewillseethisideaexploitedinsomeroutinesgivenlaterinthechapter. The “grand strategy” of virtually all modern eigensystem routines is to nudge the matrix Atowards diagonalformby a sequenceof similarity transformations, A→P −1 1·A·P1→P−1 2·P−1 1·A·P1·P2 →P−1 3·P−1 2·P−1 1·A·P1·P2·P3→etc.(11.0.14 ) 11.0 Introduction 453Sample 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 we get all the way to diagonal form, then the eigenvectors are the columns of the accumulated transformation XR=P1·P2·P3·... (11.0.15 ) Sometimeswedonotwanttogoallthewaytodiagonalform. Forexample,ifweare interestedonlyineigenvalues,noteigenvectors,it is enoughtotransformthe matrix Ato be triangular, with all elements below (or above) the diagonal zero. In this case the diagonal elements are already the eigenvalues, as you can see by mentally evaluating (11.0.2) using expansion by minors. There are two rather different sets of techniques for implementing the grand strategy (11.0.14). It turns out that they work rather well in combination, so most moderneigensystemroutinesuseboth. Thefirstsetoftechniquesconstructsindivid- ualPi’s as explicit “atomic” transformationsdesigned to performspecific tasks, for examplezeroinga particularoff-diagonalelement(Jacobitransformation, §11.1),or a whole particular row or column (Householder transformation, §11.2; elimination method, §11.5). Ingeneral,a finitesequenceofthese simpletransformationscannot completely diagonalize a matrix. There are then two choices: either use the finite sequence of transformations to go most of the way (e.g., to some special form liketridiagonal orHessenberg ,see§11.2and §11.5below)andfollowupwiththesecond set oftechniquesaboutto be mentioned;or else iterate the finite sequenceofsimple transformations over and over until the deviation of the matrix from diagonal is negligibly small. This latter approach is conceptually simplest, so we will discuss it in the next section; however, for Ngreater than ∼10, it is computationally inefficient by a roughly constant factor ∼5. The second set of techniques, called factorization methods , is more subtle. Suppose that the matrix Acan be factored into a left factor F Land a right factor FR. Then A=FL·FRorequivalently F−1 L·A=FR (11.0.16 ) Ifwenowmultiplybacktogetherthefactorsinthereverseorder,andusethesecond equation in (11.0.16) we get FR·FL=F−1 L·A·FL (11.0.17 ) which we recognize as having effected a similarity transformation on Awith the transformationmatrixbeing FL!I n§11.3and §11.6we will discuss the QR method which exploits this idea. Factorization methods also do not converge exactly in a finite number of transformations. But the better ones do converge rapidly and reliably, and, when followingan appropriateinitial reductionby simple similarity transformations,they are the methods of choice. 454 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).“Eigenpackages ofCanned Eigenroutines” Youhaveprobablygatheredbynowthatthesolutionofeigensystemsisafairly complicated business. It is. It is one of the few subjects covered in this book forwhich we do notrecommend that you avoid canned routines. On the contrary, the purpose of this chapter is precisely to give you some appreciation of what is going oninsidesuchcannedroutines,sothatyoucanmakeintelligentchoicesaboutusing them, and intelligent diagnoses when something goes wrong. Youwillfindthatalmostallcannedroutinesinusenowadaystracetheirancestry back to routines published in Wilkinson and Reinsch’s Handbook for Automatic Computation,Vol. II,LinearAlgebra [2]. Thisexcellentreference,containingpapers by a number of authors, is the Bible of the field. A public-domain implementationof theHandbook routines in FORTRAN is the EISPACK set of programs [3]. The routinesinthischapteraretranslationsofeitherthe Handbook orEISPACKroutines, so understanding these will take you a lot of the way towards understanding those canonical packages. IMSL[4]and NAG [5]each provide proprietary implementations, in FORTRAN, of what are essentially the Handbook routines. Agood“eigenpackage”willprovideseparateroutines,orseparatepathsthrough sequences of routines, for the following desired calculations: •all eigenvalues and no eigenvectors •all eigenvalues and some corresponding eigenvectors •all eigenvalues and all corresponding eigenvectors The purposeof these distinctions is to save computetime and storage; it is wasteful to calculate eigenvectors that you don’t need. Often one is interested only inthe eigenvectors corresponding to the largest few eigenvalues, or largest few in magnitude, or few that are negative. The method usually used to calculate “some” eigenvectorsistypicallymoreefficientthancalculatingalleigenvectorsifyoudesirefewer than about a quarter of the eigenvectors. A good eigenpackage also provides separate paths for each of the above calculations for each of the following special forms of the matrix: •real, symmetric, tridiagonal •real, symmetric,banded(onlya small numberof sub-and superdiagonals are nonzero) •real, symmetric •real, nonsymmetric •complex, Hermitian •complex, non-Hermitian Again,thepurposeofthesedistinctionsistosavetimeandstoragebyusingthe least general routine that will serve in any particular application. In this chapter, as a bare introduction, we give good routines for the following paths: •all eigenvalues and eigenvectors of a real, symmetric, tridiagonal matrix (§11.3) •alleigenvaluesandeigenvectorsofareal,symmetric,matrix( §11.1– §11.3) •all eigenvalues and eigenvectors of a complex, Hermitian matrix (§11.4) •alleigenvaluesandnoeigenvectorsofareal,nonsymmetricmatrix( §11.5– 11.0 Introduction 455Sample 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.6) We also discuss, in §11.7, how to obtain some eigenvectors of nonsymmetric matrices by the method of inverse iteration. Generalizedand NonlinearEigenvalue Problems Manyeigenpackagesalso dealwiththe so-called generalizedeigenproblem ,[6] A·x=λB·x (11.0.18 ) whereAandBare both matrices. Most such problems, where Bis nonsingular, can be handled by the equivalent (B−1·A)·x=λx (11.0.19 ) OftenAandBare symmetric and Bis positive definite. The matrix B−1·Ain (11.0.19) is not symmetric, but we can recover a symmetric eigenvalue problemby using the Cholesky decomposition B=L·L Tof§2.9. Multiplying equation (11.0.18) by L−1, we get C·(LT·x)=λ(LT·x)( 11.0.20 ) where C=L−1·A·(L−1)T(11.0.21 ) The matrix Cis symmetric and its eigenvaluesare the same as those of the original problem (11.0.18); its eigenfunctions are LT·x. The efficient way to form Cis first to solve the equation Y·LT=A (11.0.22 ) for the lower triangle of the matrix Y. Then solve L·C=Y (11.0.23 ) for the lower triangle of the symmetric matrix C. Another generalization of the standard eigenvalue problem is to problems nonlinear in the eigenvalue λ, for example, (Aλ2+Bλ+C)·x=0 ( 11.0.24 ) This can be turned into a linear problem by introducing an additional unknown eigenvector yand solving the 2N×2Neigensystem, /parenleftbigg 0 1 −A−1·C−A−1·B/parenrightbigg ·/parenleftbigg x y/parenrightbigg =λ/parenleftbigg x y/parenrightbigg (11.0.25 ) Thistechniquegeneralizestohigher-orderpolynomialsin λ. Apolynomialofdegree Mproduces a linear MN ×MNeigensystem (see [7]). 456 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: Stoer,J.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), Chapter 6. [1] Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra , vol. II of Handbook for Automatic Com- putation(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] IMSL Math/Library UsersManual (IMSL Inc., 2500CityWest Boulevard, HoustonTX77042).[4] NAG Fortran Library (Numerical Algorithms Group, 256 Banbury Road, Oxford OX27DE, U.K.), Chapter F02. [5] Golub,G.H.,andVanLoan,C.F.1989, MatrixComputations ,2nded.(Baltimore:JohnsHopkins University Press), §7.7. [6] Wilkinson,J.H.1965, TheAlgebraicEigenvalueProblem (NewYork:OxfordUniversityPress).[7] Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America), Chapter 13. Horn,R.A.,andJohnson,C.R.1985, MatrixAnalysis (Cambridge:CambridgeUniversityPress). 11.1 Jacobi Transformations of a Symmetric Matrix The Jacobi method consists of a sequenceof orthogonalsimilarity transforma- tions of the form of equation (11.0.14). Each transformation (a Jacobi rotation )i s just a plane rotation designed to annihilate one of the off-diagonalmatrix elements. Successive transformationsundopreviouslyset zeros, but the off-diagonalelements nevertheless get smaller and smaller, until the matrix is diagonal to machine preci-sion. Accumulating the product of the transformations as you go gives the matrix of eigenvectors, equation (11.0.15), while the elements of the final diagonal matrix are the eigenvalues. The Jacobi methodis absolutelyfoolproofforall real symmetricmatrices. For matricesofordergreaterthanabout10,say,the algorithmis slower,bya significant constant factor, than the QRmethod we shall give in §11.3. However, the Jacobi algorithm is much simpler than the more efficient methods. We thus recommend it for matrices of moderate order, where expense is not a major consideration. The basic Jacobi rotation P pqis a matrix of the form Ppq= 1 ··· c···s... 1... −s···c ··· 1 (11.1.1 ) Here all the diagonal elements are unity except for the two elements cin rows (and columns) pandq. All off-diagonalelementsare zeroexceptthe two elements sand −s. Thenumbers candsarethecosineandsineofarotationangle φ,soc 2+s2=1.