f11-7

11.7EigenvaluesorEigenvectorsbyInverseIteration 487Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).11.7 Improving Eigenvalues and/or Finding Eigenvectors by Inverse Iteration The basic idea behind inverse iteration is quite simple. Let ybe the solution of the linear system (A−τ1)·y=b (11.7.1 ) wherebis a random vector and τis close to some eigenvalue λofA. Then the solutionywill be close to the eigenvector corresponding to λ. The procedure can be iterated: Replace bbyyand solve for a new y, which will be even closer to the true eigenvector. We can see why this works by expanding both yandbas linear combinations of the eigenvectors xjofA: y=/summationdisplay jαjxjb=/summationdisplay jβjxj (11.7.2 ) Then (11.7.1) gives /summationdisplay jαj(λj−τ)xj=/summationdisplay jβjxj (11.7.3 ) so that αj=βj λj−τ(11.7.4 ) and y=/summationdisplay jβjxj λj−τ(11.7.5 ) Ifτis close to λn, say, then provided βnis not accidentally too small, ywill be approximately xn, upto a normalization. Moreover,eachiterationofthis procedure givesanotherpowerof λj−τinthedenominatorof(11.7.5). Thustheconvergence 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 ) wherebkandτkare our current guesses for some eigenvector and eigenvalue of interest (let’s say, xnandλn). Normalize bkso thatbk·bk=1. The exact eigenvector and eigenvalue satisfy A·xn=λnxn (11.7.7 ) so (A−τk1)·xn=(λn−τk)xn (11.7.8 ) Sinceyof (11.7.6)is an improvedapproximationto xn, we normalizeit and set bk+1 =y |y|(11.7.9 ) 488 Chapter11. EigensystemsSample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).We get an improved estimate of the eigenvalue by substituting our improved guess yforxnin (11.7.8). By (11.7.6), the left-hand side is bk, so calling λnour new value τk+1, we find τk+1 =τk+1 bk·y(11.7.10 ) While the above formulas look simple enough, in practice the implementation canbequitetricky. Thefirst questiontoberesolvedis whento useinverseiteration. Most of the computational load occurs in solving the linear system (11.7.6). Thusa possible strategy is first to reduce the matrix Ato a special form that allows easy solution of (11.7.6). Tridiagonal form for symmetric matrices or Hessenberg for nonsymmetric are the obvious choices. Then apply inverse iteration to generateall the eigenvectors. While this is an O(N 3)method for symmetric matrices, it is many times less efficient than the QLmethod given earlier. In fact, even the best inverse iteration packages are less efficient than the QLmethod as soon as more than about 25 percent of the eigenvectors are required. Accordingly, inverse iteration is generally used when one already has good eigenvalues and wants onlya few selected eigenvectors. You can write a simple inverse iteration routine yourself using LUdecompo- sition to solve (11.7.6). You can decide whether to use the general LUalgorithm we gave in Chapter 2 or whether to take advantage of tridiagonal or Hessenberg form. Note that, since the linear system (11.7.6) is nearly singular, you must be carefulto use a versionof LUdecompositionlike that in §2.3whichreplacesa zero pivot with a very small number. We have chosen not to give a general inverse iteration routine in this book, because it is quite cumbersome to take account of all the cases that can arise. Routinesaregiven,forexample,in [1,2]. Ifyouusethese,orwriteyourownroutine, you may appreciate the following pointers. One starts by supplying an initial value τ0for the eigenvalue λnof interest. Choose a random normalized vector b0as the initial guess for the eigenvector xn, and solve (11.7.6). The new vector yis bigger than b0by a “growth factor” |y|, which ideally shouldbe large. Equivalently,the changein the eigenvalue,which by (11.7.10)is essentially 1/|y|, should be small. The followingcases can arise: If the growth factor is too small initially, then we assume we have made a “bad” choice of random vector. This can happen not just because of a small βnin (11.7.5), but also in the case of a defective matrix, when (11.7.5) does not even apply (see, e.g., [1]or[3]for details). We go back to the beginning and choose a new initial vector. Thechange |b1−b0|mightbelessthansometolerance /epsilon1. Wecanusethis as a criterion for stopping, iterating until it is satisfied, with a maximum of 5 – 10 iterations, say. After a few iterations, if |bk+1−bk|is not decreasing rapidly enough, we can try updating the eigenvalue according to (11.7.10). If τk+1 =τk to machine accuracy, we are not going to improve the eigenvector much more and can quit. Otherwise start another cycle of iterations with the new eigenvalue. Thereasonwe donotupdatetheeigenvalueateverystepis that whenwe solve the linear system (11.7.6) by LUdecomposition, we can save the decomposition 11.7EigenvaluesorEigenvectorsbyInverseIteration 489Sample page from NUMERICAL RECIPES IN FORTRAN 77: THE ART OF SCIENTIFIC COMPUTING (ISBN 0-521-43064-X) Copyright (C) 1986-1992 by Cambridge University Press.Programs Copyright (C) 1986-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 copyin g of machine- readable files (including this one) to any servercomputer, is strictly prohibited. To order Numerical Recipes booksor CDROMs, v isit website http://www.nr.com or call 1-800-872-7423 (North America only),or send email to [email protected] (outside North Amer ica).ifτkis fixed. We only need do the backsubstitution step each time we update bk. The number of iterations we decide to do with a fixed τkis a trade-off between the quadratic convergence but O(N3)workload for updating τkat each step and the linearconvergencebut O(N2)loadforkeeping τkfixed. Ifyouhavedeterminedthe eigenvalue by one of the routines given earlier in the chapter, it is probably correctto machine accuracy anyway, and you can omit updating it. There are two differentpathologies that can arise duringinverse iteration. The firstismultipleorcloselyspacedroots. Thisismoreoftenaproblemwithsymmetric matrices. Inverseiterationwill findonlyoneeigenvectorfora giveninitialguess τ 0. A good strategy is to perturb the last few significant digits in τ0and then repeat the iteration. Usually this provides an independenteigenvector. Special steps generally have to be taken to ensure orthogonality of the linearly independent eigenvectors, whereas the Jacobi and QLalgorithms automatically yield orthogonal eigenvectors even in the case of multiple eigenvalues. The second problem, peculiar to nonsymmetric matrices, is the defective case. Unless one makes a “good” initial guess, the growth factor is small. Moreover,iteration does not improve matters. In this case, the remedy is to choose random initialvectors,solve(11.7.6)once,andquitassoonas anyvectorgivesanacceptably large growth factor. Typically only a few trials are necessary. One further complication in the nonsymmetric case is that a real matrix can have complex-conjugate pairs of eigenvalues. You will then have to use complexarithmetic to solve (11.7.6) for the complex eigenvectors. For any moderate-sized (or larger) nonsymmetric matrix, our recommendation is to avoid inverse iteration i nf a v o ro fa QRmethod that includes the eigenvector computation in complex arithmetic. You will find routines for this in [1,2]and other places. CITED REFERENCES AND FURTHER READING: Acton, F.S. 1970, Numerical Methods That Work ; 1990, corrected edition (Washington: Mathe- matical Association of America). Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra , vol. II of Handbook for Automatic Com- putation(New York: Springer-Verlag), p. 418. [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.,andBulirsch,R.1980, IntroductiontoNumericalAnalysis (NewYork:Springer-Verlag), p. 356. [3]