c11-6 (Numerical Recipes in C)

486Chapter 11.Eigensystems}}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 ofLecture 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 HessenbergMatricesRecall the following relations for the QR algorithm with shifts:Qs · (As − ks 1) = Rs(11.6.1)where Q is orthogonal and R is upper triangular, andAs+1 = Rs · QTs + ks 1= Qs · As · QTs(11.6.2)The QR transformation preserves the upper Hessenberg form of the original matrixA ≡ A1 , and the workload on such a matrix is O(n2 ) per iteration as opposedto O(n3 ) on a general matrix. As s → ∞, As converges to a form wherethe eigenvalues are either isolated on the diagonal or are eigenvalues of a 2 × 2submatrix on the diagonal.As we pointed out in §11.3, shifting is essential for rapid convergence.

A keydifference here is that a nonsymmetric real matrix can have complex eigenvalues.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (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];}}}11.6 The QR Algorithm for Real Hessenberg Matrices487This means that good choices for the shifts ks may be complex, apparentlynecessitating complex arithmetic.Complex arithmetic can be avoided, however, by a clever trick.

The trickdepends on a result analogous to the lemma we used for implicit shifts in §11.3. Thelemma we need here states that if B is a nonsingular matrix such that(11.6.3)where Q is orthogonal and H is upper Hessenberg, then Q and H are fully determinedby the first column of Q. (The determination is unique if H has positive subdiagonalelements.) The lemma can be proved by induction analogously to the proof givenfor 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 andks+1 = ks *. This gives a real matrix As+2 , whereAs+2 = Qs+1 · Qs · As · QTs · QTs+1 ·(11.6.4)The Q’s are determined byAs − ks 1 = QTs · Rs(11.6.5)As+1 = Qs · As ·As+1 − ks+1 1 = QTs+1 · Rs+1QTs(11.6.6)(11.6.7)Using (11.6.6), equation (11.6.7) can be rewrittenAs − ks+1 1 = QTs · QTs+1 · Rs+1 · Qs(11.6.8)M = (As − ks+1 1) · (As − ks 1)(11.6.9)Hence, if we defineequations (11.6.5) and (11.6.8) giveR= Q·M(11.6.10)Q = Qs+1 · Qs(11.6.11)R = Rs+1 · Rs(11.6.12)whereEquation (11.6.4) can be rewrittenAs · QT = QT · As+2(11.6.13)Thus suppose we can somehow find an upper Hessenberg matrix H such thatTTAs · Q = Q · HT(11.6.14)where Q is orthogonal.

If Q has the same first column as QT (i.e., Q has the samefirst row as Q), then Q = Q and As+2 = H.Sample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).B·Q = Q·H488Chapter 11.Eigensystemsp1 = a211 − a11 (ks + ks+1 ) + ks ks+1 + a12 a21q1 = a21 (a11 + a22 − ks − ks+1 )(11.6.15)r1 = a21 a32HenceP1 = 1 − 2w1 · wT1(11.6.16)where w1 has only its first 3 elements nonzero (cf.

equation 11.2.5). The matrixP1 · As · PT1 is therefore upper Hessenberg with 3 extra elements:××xTP1 · A1 · P1 =  x×××x×××××××××××××××××××××××××××××(11.6.17)This matrix can be restored to upper Hessenberg form without affecting the first rowby a sequence of Householder similarity transformations.

The first such Householdermatrix, P2 , acts on elements 2, 3, and 4 in the first column, annihilating elements3 and 4. This produces a matrix of the same form as (11.6.17), with the 3 extraelements appearing one column over:×××××xx××××x×××××××××××××××××××××××××(11.6.18)Proceeding in this way up to Pn−1 , we see that at each stage the Householdermatrix 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)stSample page from NUMERICAL RECIPES IN C: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43108-5)Copyright (C) 1988-1992 by Cambridge University Press.Programs Copyright (C) 1988-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 copying of machinereadable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes books,diskettes, or CDROMsvisit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to trade@cup.cam.ac.uk (outside North America).The first row of Q is found as follows.

Equation (11.6.10) shows that Q isthe orthogonal matrix that triangularizes the real matrix M. Any real matrix canbe triangularized by premultiplying it by a sequence of Householder matrices P1(acting on the first column), P2 (acting on the second column), . . . , Pn−1 . ThusQ = 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 Qsatisfying (11.6.14) whose first row is that of P1 .The Householder matrix P1 is determined by the first column of M. Since Asis upper Hessenberg, equation (11.6.9) shows that the first column of M has theform [p1 , q1, r1 , 0, ..., 0]T , where11.6 The QR Algorithm for Real Hessenberg Matrices489column of the current matrix.

Note that the preliminary matrix P1 has the samestructure as P2 , . . . , Pn−1 .The result is thatPn−1 · · · P2 · P1 · As · PT1 · PT2 · · · PTn−1 = H(11.6.19)Q = Q = Pn−1 · · · P2 · P1(11.6.20)As+2 = H(11.6.21)andThe shifts of origin at each stage are taken to be the eigenvalues of the 2 × 2matrix in the bottom right-hand corner of the current As . This givesks + ks+2 = an−1,n−1 + annks ks+1 = an−1,n−1ann − an−1,n an,n−1(11.6.22)Substituting (11.6.22) in (11.6.15), we getp1 = 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 = a21 a32(11.6.23)We have judiciously grouped terms to reduce possible roundoff when there aresmall off-diagonal elements.

Since only the ratios of elements are relevant for aHouseholder transformation, we can omit the factor a21 from (11.6.23).In summary, to carry out a double QR step we construct the Householdermatrices Pr , r = 1, . . . , n − 1. For P1 we use p1 , q1 , and r1 given by (11.6.23). Forthe 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 arithmeticoperations can be reduced by writing the nonzero elements of the 2w · wT part ofthe Householder matrix in the form(p ± s)/(±s)2w · wT =  q/(±s)  · [ 1 q/(p ± s) r/(p ± s) ](11.6.24)r/(±s)wheres2 = p2 + q 2 + r 2(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.