c11-7 (Numerical Recipes in C)

49311.7 Eigenvalues or Eigenvectors by Inverse IterationCITED 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 HopkinsUniversity Press), §7.5.11.7 Improving Eigenvalues and/or FindingEigenvectors by Inverse IterationThe basic idea behind inverse iteration is quite simple.

Let y be the solutionof 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 thesolution y will be close to the eigenvector corresponding to λ. The procedure canbe iterated: Replace b by y and solve for a new y, which will be even closer tothe true eigenvector.We can see why this works by expanding both y and b as linear combinationsof the eigenvectors xj of A:y=Xα j xjb=Xjβj xj(11.7.2)jThen (11.7.1) givesXαj (λj − τ )xj =Xjβj xj(11.7.3)jso thatαj =andy=βjλj − τX βj xjλj − τ(11.7.4)(11.7.5)jIf τ is close to λn , say, then provided βn is not accidentally too small, y will beapproximately xn , up to a normalization.

Moreover, each iteration of this proceduregives another power of λj − τ in the denominator of (11.7.5). Thus the convergenceis rapid for well-separated eigenvalues.Suppose at the kth stage of iteration we are solving the equation(A − τk 1) · y = bk(11.7.6)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).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]494Chapter 11.Eigensystemswhere bk and τk are our current guesses for some eigenvector and eigenvalue ofinterest (let’s say, xn and λn ). Normalize bk so that bk · bk = 1. The exacteigenvector and eigenvalue satisfyA · xn = λn xn(11.7.7)(A − τk 1) · xn = (λn − τk )xn(11.7.8)soWe get an improved estimate of the eigenvalue by substituting our improved guessy for xn in (11.7.8). By (11.7.6), the left-hand side is bk , so calling λn our newvalue τk+1 , we find1(11.7.10)τk+1 = τk +bk · yWhile the above formulas look simple enough, in practice the implementationcan be quite tricky.

The first question to be resolved is when to use inverse iteration.Most of the computational load occurs in solving the linear system (11.7.6). Thusa possible strategy is first to reduce the matrix A to a special form that allows easysolution of (11.7.6). Tridiagonal form for symmetric matrices or Hessenberg fornonsymmetric are the obvious choices. Then apply inverse iteration to generateall the eigenvectors.

While this is an O(N 3 ) method for symmetric matrices, itis many times less efficient than the QL method given earlier. In fact, even thebest inverse iteration packages are less efficient than the QL method as soon asmore than about 25 percent of the eigenvectors are required. Accordingly, inverseiteration 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 LU decomposition to solve (11.7.6). You can decide whether to use the general LU algorithmwe gave in Chapter 2 or whether to take advantage of tridiagonal or Hessenbergform. Note that, since the linear system (11.7.6) is nearly singular, you must becareful to use a version of LU decomposition like that in §2.3 which replaces a zeropivot 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.Routines are given, for example, in [1,2] .

If you use these, or write your own routine,you may appreciate the following pointers.One starts by supplying an initial value τ0 for the eigenvalue λn of interest.Choose a random normalized vector b0 as the initial guess for the eigenvector xn ,and solve (11.7.6). The new vector y is bigger than b0 by a “growth factor” |y|,which ideally should be large. Equivalently, the change in the eigenvalue, which by(11.7.10) is essentially 1/ |y|, should be small. The following cases can arise:• If the growth factor is too small initially, then we assume we have madea “bad” choice of random vector. This can happen not just because ofa small βn in (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 backto the beginning and choose a new initial vector.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).Since y of (11.7.6) is an improved approximation to xn , we normalize it and sety(11.7.9)bk+1 =|y|11.7 Eigenvalues or Eigenvectors by Inverse Iteration495CITED REFERENCES AND FURTHER READING:Acton, F.S. 1970, Numerical Methods That Work; 1990, corrected edition (Washington: Mathematical Association of America).Wilkinson, J.H., and Reinsch, C. 1971, Linear Algebra, vol. II of Handbook for Automatic Computation (New York: Springer-Verlag), p.

418. [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),p. 356. [3]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).• The change |b1 − b0 | might be less than some tolerance . We can use thisas a criterion for stopping, iterating until it is satisfied, with a maximumof 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 = τkto machine accuracy, we are not going to improve the eigenvector muchmore and can quit. Otherwise start another cycle of iterations with thenew eigenvalue.The reason we do not update the eigenvalue at every step is that when we solvethe linear system (11.7.6) by LU decomposition, we can save the decompositionif τk is fixed. We only need do the backsubstitution step each time we update bk .The number of iterations we decide to do with a fixed τk is a trade-off between thequadratic convergence but O(N 3 ) workload for updating τk at each step and thelinear convergence but O(N 2 ) load for keeping τk fixed.

If you have determined theeigenvalue 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 different pathologies that can arise during inverse iteration. Thefirst is multiple or closely spaced roots. This is more often a problem with symmetricmatrices. Inverse iteration will find only one eigenvector for a given initial guess τ0 .A good strategy is to perturb the last few significant digits in τ0 and then repeat theiteration. Usually this provides an independent eigenvector. Special steps generallyhave to be taken to ensure orthogonality of the linearly independent eigenvectors,whereas the Jacobi and QL algorithms automatically yield orthogonal eigenvectorseven 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 randominitial vectors, solve (11.7.6) once, and quit as soon as any vector gives an acceptablylarge growth factor. Typically only a few trials are necessary.One further complication in the nonsymmetric case is that a real matrix canhave 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 iterationin favor of a QR method that includes the eigenvector computation in complexarithmetic.

You will find routines for this in [1,2] and other places..