486 Chapter 11. Eigensystems if(1！=m)[ Interchange rows and columns. for (j=m-1;j<=n;j++)SWAP(a[i][j],a[m][j]) for (j=1;j<=n;j++)SWAP(a[j][i],a[j][m]) 1f(x)[ Carry out the elimination for(1=m+1;1<=n;1+)[ 1f((y=a[1][m-1])1=0.0)C y /=xi a[i][m-1]=y; for (j=m;j<=n;j++) http://www.nr. readable files a[i][j]-=y*a[m][j]; for (j=1;j<=n;j++) a[j][回间]+ey*a[j][i]; .com or call granted for 19881992 2 222 (including this one) 111800-672 /Cambridge NUMERICAL RECIPES IN CITED REFERENCES AND FURTHER READING: (Nort Wilkinson,J.H.,and Reinsch,C.1971,Linear Algebra,vol.Il of Handbook for Automatic Com- America users to make one paper server computer, University Press. THE putation (New York:Springer-Verlag).[1] Smith,B.T.,et al.1976,Matrix Eigensystem Routines-EISPACK Guide,2nd ed.,vol.6 of ART Lecture Notes in Computer Science (New York:Springer-Verlag).[2] Stoer,J.,and Bulirsch,R.1980,Introduction to Numerical Analysis (New York:Springer-Verlag). 9 Programs 66.5.4.[3] OF SCIENTIFIC 11.6 The QR Algorithm for Real Hessenberg Matrices 1920 MPUTING (ISBN Recall the following relations for the QR algorithm with shifts: Numerical 1021 Q·(Ag-kI)=Rg (11.6.1) 43108 where Q is orthogonal and R is upper triangular,and (outside Recipes A+1=R。·Q+k1 North Software. (11.6.2) =Q。·A。·Q7 visit website The QR transformation preserves the upper Hessenberg form of the original matrix A =A1,and the workload on such a matrix is O(n2)per iteration as opposed to O(n3)on a general matrix.As s-oo,As converges to a form where the eigenvalues are either isolated on the diagonal or are eigenvalues of a 2 x 2 submatrix on the diagonal. As we pointed out in 811.3,shifting is essential for rapid convergence.A key difference here is that a nonsymmetric real matrix can have complex eigenvalues.This
486 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). } } if (i != m) { Interchange rows and columns. for (j=m-1;j<=n;j++) SWAP(a[i][j],a[m][j]) for (j=1;j<=n;j++) SWAP(a[j][i],a[j][m]) } if (x) { Carry out the elimination. for (i=m+1;i<=n;i++) { if ((y=a[i][m-1]) != 0.0) { y /= x; a[i][m-1]=y; for (j=m;j<=n;j++) a[i][j] -= y*a[m][j]; for (j=1;j<=n;j++) a[j][m] += y*a[j][i]; } } } } } 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] Stoer, J., and Bulirsch, R. 1980, Introduction to Numerical Analysis (New York: Springer-Verlag), §6.5.4. [3] 11.6 The QR Algorithm for Real Hessenberg Matrices Recall the following relations for the QR algorithm with shifts: Qs · (As − ks1) = Rs (11.6.1) where Q is orthogonal and R is upper triangular, and As+1 = Rs · QT s + ks1 = Qs · As · QT s (11.6.2) The QR transformation preserves the upper Hessenberg form of the original matrix A ≡ A1, and the workload on such a matrix is O(n2) per iteration as opposed to O(n3) on a general matrix. As s → ∞, As converges to a form where the eigenvalues are either isolated on the diagonal or are eigenvalues of a 2×2 submatrix on the diagonal. As we pointed out in §11.3, shifting is essential for rapid convergence. A key difference here is that a nonsymmetric real matrix can have complex eigenvalues. This
11.6 The QR Algorithm for Real Hessenberg Matrices 487 means that good choices for the shifts s may be complex,apparently necessitating complex arithmetic. Complex arithmetic can be avoided.however,by a clever trick.The trick depends on a result analogous to the lemma we used for implicit shifts in $11.3.The lemma we need here states that if B is a nonsingular matrix such that B.Q=Q.H (11.6.3) where Q is orthogonal and H is upper Hessenberg,then Q and H are fully determined by the first column of Q.(The determination is unique if H has positive subdiagonal 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 OR algorithm, 菲 either with two real shifts ks and k+1,or with complex conjugate values ks and s+1=k*.This gives a real matrix As+2,where 会 A+2=Q+1Q。A..QT.oT+1 (11.6.4) The Q's are determined by 9 Ag-k1=Q.R。 (11.6.5) As+1=Q。·A。·Q (11.6.6) As+1-ks+11=Qg+1·Rs+1 (11.6.7) Using (11.6.6),equation (11.6.7)can be rewritten A。-k+11=Qf.Q+1·Rg+1·Q。 (11.6.8) Hence,if we define M=(Ag-k+11)·(Ag-ks1) (11.6.9) equations (11.6.5)and (11.6.8)give Numer R=Q·M (11.6.10) where Q=Qs+1·Q (11.6.11) R=Rg+1·Rs (11.6.12) Equation (11.6.4)can be rewritten A。·Q=Q·A+2 (11.6.13) 品 Thus suppose we can somehow find an upper Hessenberg matrix H such that A。.0=o.H (11.6.14) where Qis orthogonal.Ifhas the same first column as QT(ie.,Qhas the same first row as Q),then Q=Q and As+2 =H
11.6 The QR Algorithm for Real Hessenberg Matrices 487 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). means that good choices for the shifts ks may be complex, apparently necessitating complex arithmetic. Complex arithmetic can be avoided, however, by a clever trick. The trick depends on a result analogous to the lemma we used for implicit shifts in §11.3. The lemma we need here states that if B is a nonsingular matrix such that B · Q = Q · H (11.6.3) where Q is orthogonal and H is upper Hessenberg, then Q and H are fully determined by the first column of Q. (The determination is unique if H has positive subdiagonal 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 QR algorithm, either with two real shifts ks and ks+1, or with complex conjugate values ks and 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) The Q’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) Thus suppose we can somehow find an upper Hessenberg matrix H such that As · QT = QT · H (11.6.14) where Q is orthogonal. If QT has the same first column as QT (i.e., Q has the same first row as Q), then Q = Q and As+2 = H
488 Chapter 11.Eigensystems The first row of Q is found as follows.Equation (11.6.10)shows that Q is 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 Q is the first row of P1 since Pi is an (i-1)x(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 P. The Householder matrix P1 is determined by the first column of M.Since As is upper Hessenberg,equation (11.6.9)shows that the first column of M has the form pi,q1,r1,0,...,0],where p1=a11-a11(kg+kg+1)+kkg+1+a12a21 91=a21(a11+a22-kg-kg+1) (11.6.15) Cam T1=a21a32 RECIPES Hence P1=1-2w1·w (11.6.16) o Press. where wi has only its first 3 elements nonzero(cf.equation 11.2.5).The matrix P1.A.Pf is therefore upper Hessenberg with 3 extra elements: Programs ×1 XXXXXX × IENTIFIC X 6 P1·A1P= × (11.6.17) 千 19200 COMPUTING (ISBN This matrix can be restored to upper Hessenberg form without affecting the first row 10-621 by a sequence of Householder similarity transformations.The first such Householder matrix,P2,acts on elements 2,3,and 4 in the first column,annihilating elements Fuurggoglrion 43106 3 and 4.This produces a matrix of the same form as(11.6.17),with the 3 extra elements appearing one column over: (outside XX XX ×1 North Software. XX + + + 女 (11.6.18) × ××」 Proceeding in this way up to Pn-1,we see that at each stage the Householder matrix Pr has a vector wr that is nonzero only in elements r,r+1,and r+2. These elements are determined by the elements r,r+1,and r+2 in the(r-1)st
488 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). The first row of Q is found as follows. Equation (11.6.10) shows that Q is the orthogonal matrix that triangularizes the real matrix M. Any real matrix can be triangularized by premultiplying it by a sequence of Householder matrices P 1 (acting on the first column), P2 (acting on the second column), ... , Pn−1. Thus Q = Pn−1 ··· P2 · P1, and the first row of Q is the first row of P1 since Pi is 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 P1 is determined by the first column of M. Since As is upper Hessenberg, equation (11.6.9) shows that the first column of M has 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) Hence P1 = 1 − 2w1 · wT 1 (11.6.16) where w1 has only its first 3 elements nonzero (cf. equation 11.2.5). The matrix P1 · As · PT 1 is therefore upper Hessenberg with 3 extra elements: P1 · A1 · PT 1 = ××××××× ××××××× x ×××××× x x ××××× ×××× ××× × × (11.6.17) This matrix can be restored to upper Hessenberg form without affecting the first row by a sequence of Householder similarity transformations. The first such Householder matrix, P2, 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 ××××× x x ×××× ××× × × (11.6.18) Proceeding in this way up to Pn−1, we see that at each stage the Householder matrix Pr has a vector wr that is nonzero only in elements r, r + 1, and r + 2. These elements are determined by the elements r, r + 1, and r + 2 in the (r − 1)st
11.6 The QR Algorithm for Real Hessenberg Matrices 489 column of the current matrix.Note that the preliminary matrix P1 has the same structure as P2,...,Pn-1. The result is that Pn-1…P2P1·AgP.P…Pt-1=H (11.6.19) where H is upper Hessenberg.Thus =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 x 2 matrix in the bottom right-hand corner of the current As.This gives ksks+2 an-1,n-1+ann (11.6.22) ksks+1 an-1,n-1ann -an-1,nan,n-1 RECIPES è Substituting (11.6.22)in (11.6.15),we get Press. 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] T1=a21a32 (11.6.23) We have judiciously grouped terms to reduce possible roundoff when there are IENTIFIC( small off-diagonal elements.Since only the ratios of elements are relevant for a Householder transformation,we can omit the factor a21 from(11.6.23). In summary,to carry out a double OR step we construct the Householder matrices Pr,r=1,...,n-1.For Pi we use p1,g,and ri given by (11.6.23).For the remaining matrices,pr,r,and rr are determined by the (r,r-1),(r+1,r-1), and (r+2,r-1)elements of the current matrix.The number of arithmetic 10621 operations can be reduced by writing the nonzero elements of the 2w.wT part of the Householder matrix in the form 43106 (p±s)/(±s) (outside 2w.w q/(±s) 1 q/(p±s)r/(p±s)] (11.6.24) 腿 r/（土s】 where s2=p2+g2+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 possible ways of terminating the iteration for an eigenvalue.First,if an.n-1 becomes "negligible,"then ann is an eigenvalue.We can then delete the nth row and column
11.6 The QR Algorithm for Real Hessenberg Matrices 489 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). column of the current matrix. Note that the preliminary matrix P1 has 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) where H is 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 A s. This gives ks + ks+2 = an−1,n−1 + ann ksks+1 = an−1,n−1ann − an−1,nan,n−1 (11.6.22) 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 a 21 from (11.6.23). In summary, to carry out a double QR step we construct the Householder matrices Pr, r = 1,...,n − 1. For P1 we use p1, q1, and r1 given by (11.6.23). For the remaining matrices, pr, qr, and rr are determined by the (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 · w T part 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 possible ways of terminating the iteration for an eigenvalue. First, if a n,n−1 becomes “negligible,” then ann is an eigenvalue. We can then delete the nth row and column
490 Chapter 11.Eigensystems of the matrix and look for the next eigenvalue.Alternatively,an-1.n-2 may become negligible.In this case the eigenvalues of the 2 x 2 matrix 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. The test for convergence to an eigenvalue is combined with a test for negligible subdiagonal elements that allows splitting of the matrix into submatrices.We find the largest i such that ai.is negligible.Ifi=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 i to n(i being set to unity if there is no small subdiagonal element). 三 After determining i,the submatrix in rows i to n is examined to see if the product 81 of any two consecutive subdiagonal elements is small enough that we can work with 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 r according to equations(11.6.23) with 1 replaced by m and 2 by m+1.Ifthese were indeed the elements of the special "first"Householder matrix in a double QR step,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 require that the first two of these elements be small compared with the local diagonal elements am-1.m-1,amm and am+1.m+1.A satisfactory approximate criterion is 9 lam,m-1(lg+r) #include "nrutil.h" void hgr(float **a,int n,float wr[],float wi[]) Finds all eigenvalues of an upper Hessenberg matrix a[1..n][1..n].On input a can be exactly as output from elmhes 811.5;on output it is destroyed.The real and imaginary parts of the eigenvalues are returned in wr[1..n]and wi[1..n],respectively
490 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). of the matrix and look for the next eigenvalue. Alternatively, a n−1,n−2 may become negligible. In this case the eigenvalues of the 2 × 2 matrix 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. The test for convergence to an eigenvalue is combined with a test for negligible subdiagonal elements that allows splitting of the matrix into submatrices. We find the largest i such that ai,i−1 is 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 rowsi to n (i being set to unity if there is no small subdiagonal element). After determining i, the submatrix in rowsi to n is examined to see if the product of any two consecutive subdiagonal elements is small enough that we can work with 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 r according to equations (11.6.23) with 1 replaced by m and 2 by m+ 1. If these were indeed the elements of the special “first” Householder matrix in a double QR step, 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 require that the first two of these elements be small compared with the local diagonal elements am−1,m−1, amm and am+1,m+1. A satisfactory approximate criterion is |am,m−1|(|q| + |r|) |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 30 unsuccessful iterations, the routine reports failure. The operation count for the QR algorithm described here is ∼ 5k 2 per iteration, where k is the current size of the matrix. The typical average number of iterations per eigenvalue is ∼ 1.8, so the total operation count for all the eigenvalues is ∼ 3n3. This estimate neglects any possible efficiency due to splitting or sparseness of the matrix. The following routine hqr is based algorithmically on the above description, in turn following the implementations in [1,2]. #include #include "nrutil.h" void hqr(float **a, int n, float wr[], float wi[]) Finds all eigenvalues of an upper Hessenberg matrix a[1..n][1..n]. On input a can be exactly as output from elmhes §11.5; on output it is destroyed. The real and imaginary parts of the eigenvalues are returned in wr[1..n] and wi[1..n], respectively. {
11.6 The QR Algorithm for Real Hessenberg Matrices 491 int nn,m,1,k,j,its,i,mmin; float z,y,x,w,v,u,t,s,r,q,p,anormi anorm=0.0; Compute matrix norm for possible use in lo- for(i=1:i=1){ Begin search for next eigenvalue. 1ts=0; do for(1=nn;1>=2;1-)[ Begin iteration:look for single small subdi- s=fabs(a[1-1][1-1])+fabs(a[1][1]); agonal element. if (s ==0.0)s=anorm; if((f1oat)(fabs(a[1][1-1])+s)=s){ 18881992 a[1][1-1]=0.0; break; 100 (including this one) 2 872 x=a[nn][nn]; if (1 ==nn){ One root found. 7423 from NUMERICAL RECIPES IN wr[nn]=x+t; wi[nn--]=0.0; else y=a[nn-1][nn-1]; w=a[nn][nn-1]*a[nn-1][nn]; (North America to to to make one paper THE if(1==(nn-1))[ Two roots found... /Cambridge University Press.Programs ART p=0.5*(y-x); q-p*p+v; z=sqrt(fabs(q)); x+t; send 1f(q>=0.0){ ...a real pair. z=p+SIGN(z,p); wr [nn-1]=wr [nn]=x+z; Copyright (C) if (z)wr[nn]=x-w/z; wi[nn-1]=wi[nn]=0.0; else ...a complex pair. wr [nn-1]=wr [nn]=x+p; to directcustser 1988-1992 by Numerical Recipes OF SCIENTIFIC COMPUTING(ISBN wi[nn-1]=-(wi[nn]=z); nn-=2; else No roots found.Continue iteration v@cambr 10-621 if (its =30)nrerror("Too many iterations in hgr"); 1f(1ts==10I11ts==20){ Form exceptional shift. t+=X; -43108 for (i=1;i=1;m-){ Form shift and then look for z=a[m][m]; 2 consecutive small sub- 工=X=Z: diagonal elements. 8=y-z; p=(r*s-w)/a[m+1][m]+a[m][m+1]: Equation (11.6.23) qa[m+1][m+1]-z-r-s; r=a[m+2][m+1]; s=fabs(p)+fabs(q)+fabs(r); Scale to prevent overflow or p/=s; underflow. 9/as; r /s; if (m ==1)break;
11.6 The QR Algorithm for Real Hessenberg Matrices 491 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). int nn,m,l,k,j,its,i,mmin; float z,y,x,w,v,u,t,s,r,q,p,anorm; anorm=0.0; Compute matrix norm for possible use in lofor (i=1;i= 1) { Begin search for next eigenvalue. its=0; do { for (l=nn;l>=2;l--) { Begin iteration: look for single small subdis=fabs(a[l-1][l-1])+fabs(a[l][l]); agonal element. if (s == 0.0) s=anorm; if ((float)(fabs(a[l][l-1]) + s) == s) { a[l][l-1]=0.0; break; } } x=a[nn][nn]; if (l == nn) { One root found. wr[nn]=x+t; wi[nn--]=0.0; } else { y=a[nn-1][nn-1]; w=a[nn][nn-1]*a[nn-1][nn]; if (l == (nn-1)) { Two roots found... p=0.5*(y-x); q=p*p+w; z=sqrt(fabs(q)); x += t; if (q >= 0.0) { ...a real pair. z=p+SIGN(z,p); wr[nn-1]=wr[nn]=x+z; if (z) wr[nn]=x-w/z; wi[nn-1]=wi[nn]=0.0; } else { ...a complex pair. wr[nn-1]=wr[nn]=x+p; wi[nn-1]= -(wi[nn]=z); } nn -= 2; } else { No roots found. Continue iteration. if (its == 30) nrerror("Too many iterations in hqr"); if (its == 10 || its == 20) { Form exceptional shift. t += x; for (i=1;i=l;m--) { Form shift and then look for 2 consecutive 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=fabs(p)+fabs(q)+fabs(r); Scale to prevent overflow or p /= s; underflow. q /= s; r /= s; if (m == l) break;
492 Chapter 11.Eigensystems u=fabs(a [m][m-1])*(fabs(q)+fabs(r)); v=fabs(p)*(fabs(a [m-1][m-1])+fabs(z)+fabs(a [m+1][m+1])); if ((float)(u+v)==v)break; Equation (11.6.26). for(i=m+2;1<=nn;i++)[ a[i][i-2]=0.0; 1f(11=(m+2))a[i][i-3]=0.0; for (k=m;k<=nn-1;k++){ Double QR step on rows 1 to nn and columns m to nn. if (k !m){ p=a[k][k-1]: Begin setup of Householder http://www.nr. q=a0k+1][k-1]; vector. r=0.0; 1f(kI=(nn-1))r=a[k+2][k-1]; 83 if ((x=fabs(p)+fabs(g)+fabs(r))!=0.0){ 鱼 198891992 p/=x; Scale to prevent overflow or q/=x; underflow. 1.800 (including this one) r/=x; if ((s=SIGN(sqrt(p*p+q*q+r*r),p))!-0.0){ 7423 to to tusers to make one paper by Cambridge University Press.Programs from NUMERICAL RECIPES IN 1f(k=m){ 1f(11=m) ak][k-1]=-ak]k-1] else a0k][k-1]=-g*x; (North America THE p+=s; Equations (11.6.24). x=p/s; y=q/s: only). z=r/s; q/=p: send r /=pi for (i=k;i<=nn;++) Row modification. strictly prohibited. p=a[k][j]+q*ak+1][j]; if (k !=(nn-1)){ email to directcustsen Copyright(C) p+=rta[k+2][j]; ak+2][j]-=p*z; ak+1][]-=p*y ART OF SCIENTIFIC COMPUTING(ISBN 0-521 a[k][j]-p*x; mmin nn<k+3 nn k+3; for (i=l;i<=mmin;i++){ Column modification. p=x*a[i]k]+y*a[1][k+1]; if (k !(nn-1)){ @cambridge.org(outside North America). To order Numerical Recipes books or 1988-1992 by Numerical Recipes p+=z*a[i][k+2] 1-431085 a[1][k+2]-=p*r; a[1][k+1]-=p*q a[i][k]-=p; Software. ying of machine while (1<nn-1) 22
492 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). u=fabs(a[m][m-1])*(fabs(q)+fabs(r)); v=fabs(p)*(fabs(a[m-1][m-1])+fabs(z)+fabs(a[m+1][m+1])); if ((float)(u+v) == v) break; Equation (11.6.26). } for (i=m+2;i<=nn;i++) { a[i][i-2]=0.0; if (i != (m+2)) a[i][i-3]=0.0; } for (k=m;k<=nn-1;k++) { Double QR step on rows l to nn and columns m to nn. if (k != m) { p=a[k][k-1]; Begin setup of Householder q=a[k+1][k-1]; vector. r=0.0; if (k != (nn-1)) r=a[k+2][k-1]; if ((x=fabs(p)+fabs(q)+fabs(r)) != 0.0) { p /= x; Scale to prevent overflow or q /= x; underflow. r /= x; } } if ((s=SIGN(sqrt(p*p+q*q+r*r),p)) != 0.0) { if (k == m) { if (l != m) a[k][k-1] = -a[k][k-1]; } else a[k][k-1] = -s*x; p += s; Equations (11.6.24). x=p/s; y=q/s; z=r/s; q /= p; r /= p; for (j=k;j<=nn;j++) { Row modification. p=a[k][j]+q*a[k+1][j]; if (k != (nn-1)) { p += r*a[k+2][j]; a[k+2][j] -= p*z; } a[k+1][j] -= p*y; a[k][j] -= p*x; } mmin = nn<k+3 ? nn : k+3; for (i=l;i<=mmin;i++) { Column modification. p=x*a[i][k]+y*a[i][k+1]; if (k != (nn-1)) { p += z*a[i][k+2]; a[i][k+2] -= p*r; } a[i][k+1] -= p*q; a[i][k] -= p; } } } } } } while (l < nn-1); } }
11.7 Eigenvalues or Eigenvectors by Inverse Iteration 493 CITED REFERENCES AND FURTHER READING: Wilkinson,J.H.,and Reinsch,C.1971,Linear Algebra,vol.Il of Handbook for Automatic Com- putation (New York:Springer-Verlag).[1] Golub,G.H.,and Van Loan,C.F.1989,Matrix Computations,2nd ed.(Baltimore:Johns Hopkins University Press),87.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.7 Improving Eigenvalues and/or Finding Eigenvectors by Inverse Iteration The basic idea behind inverse iteration is quite simple.Let y be the solution of the linear system (A-T1).y=b (11.7.1) 令 where b is a random vector and r is close to some eigenvalue A of A.Then the 需 solution y will be close to the eigenvector corresponding to A.The procedure can be iterated:Replace b by y and solve for a new y,which will be even closer to the true eigenvector. Program We can see why this works by expanding both y and b as linear combinations of the eigenvectors xj of A: 色 OF SCIENTIFIC( 61 y=b=∑Bx (11.7.2) Then (11.7.1)gives ∑a-Ts=∑ (11.7.3) Numerica 10.621 so that 、86六 43106 月 =-T (11.7.4) (outside and y-刀月当 (11.7.5) 入-T If r is close to An,say,then provided Bn is not accidentally too small,y will be approximately xn,up to a normalization.Moreover,each iteration of this procedure gives another power ofj-in the denominator of(11.7.5).Thus the convergence 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)
11.7 Eigenvalues or Eigenvectors by Inverse Iteration 493 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). 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] Golub, G.H., and Van Loan, C.F. 1989, Matrix Computations, 2nd ed. (Baltimore: Johns Hopkins 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.7 Improving Eigenvalues and/or Finding Eigenvectors by Inverse Iteration The basic idea behind inverse iteration is quite simple. Let y be the solution of the linear system (A − τ1) · y = b (11.7.1) where b is a random vector and τ is close to some eigenvalue λ of A. Then the solution y will be close to the eigenvector corresponding to λ. The procedure can be iterated: Replace b by y and solve for a new y, which will be even closer to the true eigenvector. We can see why this works by expanding both y and b as linear combinations of the eigenvectors xj of A: y = j αjxj b = j βjxj (11.7.2) Then (11.7.1) gives j αj (λj − τ)xj = j βjxj (11.7.3) so that αj = βj λj − τ (11.7.4) and y = j βjxj λj − τ (11.7.5) If τ is close to λn, say, then provided βn is not accidentally too small, y will be approximately xn, up to a normalization. Moreover, each iteration of this procedure gives another power of λj − τ in the denominator of (11.7.5). Thus the convergence 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)
