Rayleigh quotient iteration

From testwiki
Jump to navigation Jump to search

Rayleigh quotient iteration is an eigenvalue algorithm which extends the idea of the inverse iteration by using the Rayleigh quotient to obtain increasingly accurate eigenvalue estimates.

Rayleigh quotient iteration is an iterative method, that is, it delivers a sequence of approximate solutions that converges to a true solution in the limit. Very rapid convergence is guaranteed and no more than a few iterations are needed in practice to obtain a reasonable approximation. The Rayleigh quotient iteration algorithm converges cubically for Hermitian or symmetric matrices, given an initial vector that is sufficiently close to an eigenvector of the matrix that is being analyzed.

Algorithm

The algorithm is very similar to inverse iteration, but replaces the estimated eigenvalue at the end of each iteration with the Rayleigh quotient. Begin by choosing some value μ0 as an initial eigenvalue guess for the Hermitian matrix A. An initial vector b0 must also be supplied as initial eigenvector guess.

Calculate the next approximation of the eigenvector bi+1 by bi+1=(AμiI)1bi(AμiI)1bi, where I is the identity matrix, and set the next approximation of the eigenvalue to the Rayleigh quotient of the current iteration equal to
μi+1=bi+1*Abi+1bi+1*bi+1.

To compute more than one eigenvalue, the algorithm can be combined with a deflation technique.Template:Citation needed

Note that for very small problems it is beneficial to replace the matrix inverse with the adjugate, which will yield the same iteration because it is equal to the inverse up to an irrelevant scale (the inverse of the determinant, specifically). The adjugate is easier to compute explicitly than the inverse (though the inverse is easier to apply to a vector for problems that aren't small), and is more numerically sound because it remains well defined as the eigenvalue converges.

Example

Consider the matrix

A=[123121321]

for which the exact eigenvalues are λ1=3+5, λ2=35 and λ3=2, with corresponding eigenvectors

v1=[1φ11], v2=[1φ1] and v3=[101].

(where φ=1+52 is the golden ratio).

The largest eigenvalue is λ15.2361 and corresponds to any eigenvector proportional to v1[10.61801].

We begin with an initial eigenvalue guess of

b0=[111],μ0=200.

Then, the first iteration yields

b1[0.579270.573480.57927],μ15.3355

the second iteration,

b2[0.646760.404220.64676],μ25.2418

and the third,

b3[0.647930.400450.64793],μ35.2361

from which the cubic convergence is evident.

Octave implementation

The following is a simple implementation of the algorithm in Octave.

function x = rayleigh(A, epsilon, mu, x)
  x = x / norm(x);
  % the backslash operator in Octave solves a linear system
  y = (A - mu * eye(rows(A))) \ x; 
  lambda = y' * x;
  mu = mu + 1 / lambda
  err = norm(y - lambda * x) / norm(y)

  while err > epsilon
    x = y / norm(y);
    y = (A - mu * eye(rows(A))) \ x;
    lambda = y' * x;
    mu = mu + 1 / lambda
    err = norm(y - lambda * x) / norm(y)
  end

end

See also

References

  • Lloyd N. Trefethen and David Bau, III, Numerical Linear Algebra, Society for Industrial and Applied Mathematics, 1997. Template:Isbn.
  • Rainer Kress, "Numerical Analysis", Springer, 1991. Template:Isbn

Template:Numerical linear algebra