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

团购合买资源类别：文库，文档格式：PDF，文档页数：5，文件大小：58.55KB

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.

11.5 Reduction of a General Matrix to Hessenberg Form 483 Apart from referring you to the literature,and to the collected routines in [1,2].we are going to sidestep the problem of eigenvectors,giving algorithms for eigenvalues only.If you require just a few eigenvectors,you can read 811.7 and consider finding them by inverse iteration.We consider the problem of finding all eigenvectors of a nonsymmetric matrix as lying beyond the scope of this book. Balancing The sensitivity of eigenvalues to rounding errors during the execution of some algorithms can be reduced by the procedure of balancing.The errors in the eigensystem found by a numerical procedure are generally proportional to the g Euclidean norm of the matrix,that is,to the square root of the sum of the squares of the elements.The idea of balancing is to use similarity transformations to make corresponding rows and columns of the matrix have comparable norms,thus reducing the overall norm of the matrix while leaving the eigenvalues unchanged. A symmetric matrix is already balanced. 之物 Balancing is a procedure with of order N2 operations.Thus,the time taken by the procedure balanc,given below,should never be more than a few percent of the total time required to find the eigenvalues.It is therefore recommended that you ahways balance nonsymmetric matrices.It never hurts,and it can substantially improve the accuracy of the eigenvalues computed for a badly balanced matrix. 当寻%n 9 The actual algorithm used is due to Osborne.as discussed in 111.It consists of a sequence of similarity transformations by diagonal matrices D.To avoid introducing rounding errors during the balancing process,the elements of D are restricted to be exact powers of the radix base employed for floating-point arithmetic(i.e.,2 for most OF SCIENTIFIC machines,but 16 for IBM mainframe architectures).The output is a matrix that is balanced in the norm given by summing the absolute magnitudes of the matrix elements.This is more efficient than using the Euclidean norm,and equally effective: A large reduction in one norm implies a large reduction in the other. Note that if the off-diagonal elements of any row or column of a matrix are all zero,then the diagonal element is an eigenvalue.If the eigenvalue happens to be ill-conditioned(sensitive to small changes in the matrix elements),it will have relatively large errors when determined by the routine hgr(811.6).Had we merely Numerica 10621 inspected the matrix beforehand,we could have determined the isolated eigenvalue 43106 exactly and then deleted the corresponding row and column from the matrix.You should consider whether such a pre-inspection might be useful in your application. (For symmetric matrices,the routines we gave will determine isolated eigenvalues (outside Recipes 腿 accurately in all cases.) The routine balanc does not keep track of the accumulated similarity trans- North formation of the original matrix,since we will only be concerned with finding eigenvalues of nonsymmetric matrices,not eigenvectors.Consult [1-3]if you want to keep track of the transformation. #include #define RADIX 2.0 void balanc(float *a ,int n) Given a matrix a[1..n][1..n],this routine replaces it by a balanced matrix with identical eigenvalues.A symmetric matrix is already balanced and is unaffected by this procedure.The parameter RADIX should be the machine's floating-point radix

11.5 Reduction of a General Matrix to Hessenberg Form 483 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). Apart from referring you to the literature, and to the collected routines in [1,2], we are going to sidestep the problem of eigenvectors, giving algorithms for eigenvalues only. If you require just a few eigenvectors, you can read §11.7 and consider finding them by inverse iteration. We consider the problem of finding all eigenvectors of a nonsymmetric matrix as lying beyond the scope of this book. Balancing The sensitivity of eigenvalues to rounding errors during the execution of some algorithms can be reduced by the procedure of balancing. The errors in the eigensystem found by a numerical procedure are generally proportional to the Euclidean norm of the matrix, that is, to the square root of the sum of the squares of the elements. The idea of balancing is to use similarity transformations to make corresponding rows and columns of the matrix have comparable norms, thus reducing the overall norm of the matrix while leaving the eigenvalues unchanged. A symmetric matrix is already balanced. Balancing is a procedure with of order N 2 operations. Thus, the time taken by the procedure balanc, given below, should never be more than a few percent of the total time required to find the eigenvalues. It is therefore recommended that you always balance nonsymmetric matrices. It never hurts, and it can substantially improve the accuracy of the eigenvalues computed for a badly balanced matrix. The actual algorithm used is due to Osborne, as discussed in [1]. It consists of a sequence of similarity transformations by diagonal matrices D. To avoid introducing rounding errors during the balancing process, the elements of D are restricted to be exact powers of the radix base employed for floating-point arithmetic (i.e., 2 for most machines, but 16 for IBM mainframe architectures). The output is a matrix that is balanced in the norm given by summing the absolute magnitudes of the matrix elements. This is more efficient than using the Euclidean norm, and equally effective: A large reduction in one norm implies a large reduction in the other. Note that if the off-diagonal elements of any row or column of a matrix are all zero, then the diagonal element is an eigenvalue. If the eigenvalue happens to be ill-conditioned (sensitive to small changes in the matrix elements), it will have relatively large errors when determined by the routine hqr (§11.6). Had we merely inspected the matrix beforehand, we could have determined the isolated eigenvalue exactly and then deleted the corresponding row and column from the matrix. You should consider whether such a pre-inspection might be useful in your application. (For symmetric matrices, the routines we gave will determine isolated eigenvalues accurately in all cases.) The routine balanc does not keep track of the accumulated similarity transformation of the original matrix, since we will only be concerned with finding eigenvalues of nonsymmetric matrices, not eigenvectors. Consult [1-3] if you want to keep track of the transformation. #include #define RADIX 2.0 void balanc(float **a, int n) Given a matrix a[1..n][1..n], this routine replaces it by a balanced matrix with identical eigenvalues. A symmetric matrix is already balanced and is unaffected by this procedure. The parameter RADIX should be the machine’s floating-point radix

11.5 Reduction of a General Matrix to Hessenberg Form 485 required elements in a column of the matrix.Householder reduction to Hessenberg form is in fact an accepted technique.An alternative,however,is a procedure analogous to Gaussian elimination with pivoting.We will use this elimination procedure since it is about a factor of 2 more efficient than the Householder method, and also since we want to teach you the method.It is possible to construct matrices for which the Householder reduction,being orthogonal,is stable and elimination is not,but such matrices are extremely rare in practice. Straight Gaussian elimination is not a similarity transformation of the matrix Accordingly,the actual elimination procedure used is slightly different.Before the rth stage,the original matrix A =A1 has become Ar,which is upper Hessenberg in its first r-1 rows and columns.The rth stage then consists of the following sequence of operations: Find the element of maximum magnitude in the rth column below the diagonal.If it is zero,skip the next two"bullets"and the stage is done. Cam Otherwise,suppose the maximum element was in row r'. Interchange rows r'and r+1.This is the pivoting procedure.To make the permutation a similarity transformation,also interchange columns r' and r+1. 三d 令 For i=r+2,r+3....,N,compute the multiplier dir Press. ni,r+1三 0r+1, Subtract ni.+times row r+1 from row i.To make the elimination a similarity transformation,also add ni+times column i to columnr+1. A total of N-2 such stages are required. IENTIFIC When the magnitudes of the matrix elements vary over many orders,you should try to rearrange the matrix so that the largest elements are in the top left-hand corner. This reduces the roundofferror,since the reduction proceeds from left to right. Since we are concerned only with eigenvalues,the routine elmhes does not 最管 keep track of the accumulated similarity transformation.The operation count is about 5N3/6 for large N. 10621 #include #define SWAP(g,h){y=(g);(g)=(h);(h)=y;} Numerical Recipes 43108 void elmhes(float a,int n) Reduction to Hessenberg form by the elimination method.The real,nonsymmetric matrix (outside a[1..n][1..n]is replaced by an upper Hessenberg matrix with identical eigenvalues.Rec- ommended,but not required,is that this routine be preceded by balanc.On output,the North Software. Hessenberg matrix is in elements a[i][j]with ij+1 are to be thought of as zero,but are returned with random values. Ame int m,j,i; float y,xi for (m=2;mfabs(x)){ x=a[j][m-1]; 1=j;

11.5 Reduction of a General Matrix to Hessenberg Form 485 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). required elements in a column of the matrix. Householder reduction to Hessenberg form is in fact an accepted technique. An alternative, however, is a procedure analogous to Gaussian elimination with pivoting. We will use this elimination procedure since it is about a factor of 2 more efficient than the Householder method, and also since we want to teach you the method. It is possible to construct matrices for which the Householder reduction, being orthogonal, is stable and elimination is not, but such matrices are extremely rare in practice. Straight Gaussian elimination is not a similarity transformation of the matrix. Accordingly, the actual elimination procedure used is slightly different. Before the rth stage, the original matrix A ≡ A1 has become Ar, which is upper Hessenberg in its first r − 1 rows and columns. The rth stage then consists of the following sequence of operations: • Find the element of maximum magnitude in the rth column below the diagonal. If it is zero, skip the next two “bullets” and the stage is done. Otherwise, suppose the maximum element was in row r . • Interchange rows r and r + 1. This is the pivoting procedure. To make the permutation a similarity transformation, also interchange columns r and r + 1. • For i = r + 2, r + 3,...,N, compute the multiplier ni,r+1 ≡ air ar+1,r Subtract ni,r+1 times row r + 1 from row i. To make the elimination a similarity transformation, also add ni,r+1 times column i to column r + 1. A total of N − 2 such stages are required. When the magnitudes of the matrix elements vary over many orders, you should try to rearrange the matrix so that the largest elements are in the top left-hand corner. This reduces the roundoff error, since the reduction proceeds from left to right. Since we are concerned only with eigenvalues, the routine elmhes does not keep track of the accumulated similarity transformation. The operation count is about 5N3/6 for large N. #include #define SWAP(g,h) {y=(g);(g)=(h);(h)=y;} void elmhes(float **a, int n) Reduction to Hessenberg form by the elimination method. The real, nonsymmetric matrix a[1..n][1..n] is replaced by an upper Hessenberg matrix with identical eigenvalues. Recommended, but not required, is that this routine be preceded by balanc. On output, the Hessenberg matrix is in elements a[i][j] with i ≤ j+1. Elements with i > j+1 are to be thought of as zero, but are returned with random values. { int m,j,i; float y,x; for (m=2;m fabs(x)) { x=a[j][m-1]; i=j;

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

点击进入文档下载页（PDF格式）

已到末页，全文结束

点击下载（PDF格式）

浏览记录