11.4 Hermitian Matrices 481 for(i=m-1;1>=1;1--)[ A plane rotation as in the origi- f=s*e[i]; nal OL.followed by Givens b=c*e[i]; rotations to restore tridiag- e[i+1]=(r=pythag(f,g)); onal form. 1f(r==0.0){ Recover from underflow. d[i+1]-=p; e[m]=0.0; break; s-f/r; c=g/r; gd[1+1]-p; r=(d[i]-g)*s+2.0*c*b; d[i+1]=g+(p=s*r); g=c*r-b; 83g /Next loop can be omitted if eigenvectors not wanted*/ 19881992 for(k=1;k=1)continue; d[1]-=p; e[1]=g; e[m]=0.0; (North America 州bMe se make one paper University Press. THE while (m !1); ART 22 ictly proh Programs CITED REFERENCES AND FURTHER READING: Acton,F.S.1970,Numerica/Methods That Work;1990,corrected edition (Washington:Mathe- matical Association of America).pp.331-335.[1] Wilkinson,J.H.,and Reinsch,C.1971,Linear Algebra,vol.ll of Handbook for Automatic Com- OF SCIENTIFIC COMPUTING(ISBN putation (New York:Springer-Verlag).[2] Smith,B.T..et al.1976,Matrix Eigensystem Routines-EISPACK Guide,2nd ed.,vol.6 of 1988-1992 Lecture Notes in Computer Science (New York:Springer-Verlag).[3] Stoer,J.,and Bulirsch,R.1980,Introduction to Numerical Analysis(New York:Springer-Verlag). Fuurggoglrion Numerica 10621 56.6.6.[4 43106 11.4 Hermitian Matrices (outside Recipes North Software. The complex analog of a real,symmetric matrix is a Hermitian matrix, satisfying equation(11.0.4).Jacobi transformations can be used to find eigenvalues and eigenvectors,as also can Householder reduction to tridiagonal form followed by visit website machine OL iteration.Complex versions of the previous routines jacobi,tred2,and tqli are quite analogous to their real counterparts.For working routines,consult [1.21. An alternative,using the routines in this book,is to convert the Hermitian problem to a real,symmetric one:If C=A+iB is a Hermitian matrix,then the n x n complex eigenvalue problem (A+iB)·(u+iv)=λ(u+iv) (11.4.1)
11.4 Hermitian Matrices 481 Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). for (i=m-1;i>=l;i--) { A plane rotation as in the original QL, followed by Givens rotations to restore tridiagonal form. f=s*e[i]; b=c*e[i]; e[i+1]=(r=pythag(f,g)); if (r == 0.0) { Recover from underflow. d[i+1] -= p; e[m]=0.0; break; } s=f/r; c=g/r; g=d[i+1]-p; r=(d[i]-g)*s+2.0*c*b; d[i+1]=g+(p=s*r); g=c*r-b; /* Next loop can be omitted if eigenvectors not wanted*/ for (k=1;k= l) continue; d[l] -= p; e[l]=g; e[m]=0.0; } } while (m != l); } } CITED REFERENCES AND FURTHER READING: Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America), pp. 331–335. [1] Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol. II of Handbook for Automatic Computation (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] Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: 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 transformations can be used to find eigenvalues and eigenvectors, as also can Householder reduction to tridiagonal form followed by QL iteration. Complex versions of the previous routines jacobi, tred2, and tqli are quite analogous to their real counterparts. For working routines, 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 + iB is a Hermitian matrix, then the n × n complex eigenvalue problem (A + iB) · (u + iv) = λ(u + iv) (11.4.1)
482 Chapter 11. Eigensystems is equivalent to the 2n x 2n real problem [&][日=8 (11.4.2) Note that the 2n x 2n matrix in (11.4.2)is symmetric:AT =A and BT=-B if C is Hermitian. Corresponding to a given eigenvalue A,the vector [ (11.4.3) 超君 is also an eigenvector,as you can verify by writing out the two matrix equa- tions implied by(ll.4.2).Thus ifλ1,2,..,入are the eigenvalues of C,then the 2n eigenvalues of the augmented problem (11.4.2)are A1,A1,A2,A2,..., ICAL An,An;each,in other words,is repeated twice.The eigenvectors are pairs of the form u+iv and i(u+iv);that is,they are the same up to an inessential phase.Thus RECIPES we solve the augmented problem(11.4.2),and choose one eigenvalue and eigenvector from each pair.These give the eigenvalues and eigenvectors of the original matrix C 9 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. 9 CITED REFERENCES AND FURTHER READING: 2色g么 Wilkinson,J.H.,and Reinsch,C.1971,Linear Algebra,vol.Il 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 6 Lecture Notes in Computer Science (New York:Springer-Verlag).[2] (ISBN 11.5 Reduction of a General Matrix to Numerica 10.621 Hessenberg Form Recipes 43108 The algorithms for symmetric matrices,given in the preceding sections,are highly satisfactory in practice.By contrast,it is impossible to design equally (outside satisfactory algorithms for the nonsymmetric case.There are two reasons for this North Software. First,the eigenvalues ofa nonsymmetric matrix can be very sensitive to small changes in the matrix elements.Second,the matrix itself can be defective,so that there is no complete set of eigenvectors.We emphasize that these difficulties are intrinsic properties of certain nonsymmetric matrices,and no numerical procedure can"cure" them.The best we can hope for are procedures that don't exacerbate such problems. The presence of rounding error can only make the situation worse.With finite- precision arithmetic,one cannot even design a foolproof algorithm to determine whether a given matrix is defective or not.Thus current algorithms generally try to find a complete set of eigenvectors,and rely on the user to inspect the results.If any eigenvectors are almost parallel,the matrix is probably defective
482 Chapter 11. Eigensystems Permission is granted for internet users to make one paper copy for their own personal use. Further reproduction, or any copyin Copyright (C) 1988-1992 by Cambridge University Press. Programs Copyright (C) 1988-1992 by Numerical Recipes Software. Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5) g of machinereadable files (including this one) to any server computer, is strictly prohibited. To order Numerical Recipes books or CDROMs, visit website http://www.nr.com or call 1-800-872-7423 (North America only), or send email to directcustserv@cambridge.org (outside North America). is equivalent to the 2n × 2n real problem A −B B A · u v = λ u v (11.4.2) Note that the 2n × 2n matrix in (11.4.2) is symmetric: AT = A and BT = −B if C is Hermitian. Corresponding to a given eigenvalue λ, the vector −v u (11.4.3) is also an eigenvector, as you can verify by writing out the two matrix equations implied by (11.4.2). Thus if λ1, λ2,...,λn are the eigenvalues of C, then the 2n eigenvalues 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 form u + iv and i(u + iv); that is, they are the same up to an inessential phase. Thus we solve the augmented problem (11.4.2), and choose one eigenvalue and eigenvector from each pair. These give the eigenvalues and eigenvectors of the original matrix 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. CITED REFERENCES AND FURTHER READING: Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol. II of Handbook for Automatic Computation (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] 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 equally satisfactory algorithms for the nonsymmetric case. There are two reasons for this. First, the eigenvalues of a nonsymmetric matrix can be very sensitive to small changes in the matrix elements. Second, the matrix itself can be defective, so that there is no complete set of eigenvectors. We emphasize that these difficulties are intrinsic properties of certain nonsymmetric matrices, and no numerical procedure can “cure” them. The best we can hope for are procedures that don’t exacerbate such problems. The presence of rounding error can only make the situation worse. With finiteprecision arithmetic, one cannot even design a foolproof algorithm to determine whether a given matrix is defective or not. Thus current algorithms generally try to find a complete set of eigenvectors, and rely on the user to inspect the results. If any eigenvectors are almost parallel, the matrix is probably defective.