《数字信号处理》教学参考资料（Numerical Recipes in C，The Art of Scientific Computing Second Edition）Chapter 11.6.pdf

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 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 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 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)<pl(lam+1.m+1+lamml+lam-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 9起N 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 kg+kg+1=1.5×(0an,n-1+lan-1,n-2l0 (11.6.27) ksks+1=(lan,n-1+lan-1.n-2)2 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. 61 Numerica 10621 4310 The operation count for the QR algorithm described here is~5k2 per iteration, where k is the current size of the matrix.The typical average number of iterations per (outside Recipes 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. North Software. The following routine hor is based algorithmically on the above description, in turn following the implementations in [1,2]. #include <math.h> #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 <math.h> #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. {