🌀Eigenvalue Problems

Inverse Iteration and Shifts

Inverse iteration, shifted inverse iteration, Rayleigh quotient iteration and deflation

The first number a structural engineer asks for

Every seismic design code — Eurocode 8, ASCE 7, IS 1893 — starts the same way: estimate the fundamental period

T_1

of the structure, the period of its slowest natural mode. The design spectrum then gives the acceleration the structure must resist for that period, and the base shear follows. The codes even offer rules of thumb for a first guess — a moment-frame building of

N

storeys has

T_1

roughly

0.1N

seconds — but a computed value is required for anything important. That period corresponds to the smallest eigenvalue of

K\phi = \omega^2 M\phi

. A finite-element model of a tall building has hundreds of thousands of degrees of freedom. The power method converges to the largest eigenvalue, some local jiggle of one element at thousands of hertz. What we need is at the opposite end of the spectrum, and we need it to four digits.

💡Inverse iteration

A\mathbf v = \lambda\mathbf v

then

A^{-1}\mathbf v = \lambda^{-1}\mathbf v

: same eigenvectors, reciprocal eigenvalues. The smallest eigenvalue of

A

is the largest of

A^{-1}

. So run the power method on

A^{-1}

\text{solve } A\mathbf y = \mathbf x_k, \qquad \mathbf x_{k+1} = \mathbf y/\|\mathbf y\|.

Never form

A^{-1}

. Factorise

A = LU

(or Cholesky) once; every iteration is then two triangular solves (Unit 3). The rate is

|\lambda_1/\lambda_2|

— for the vibrating string,

1/4

per step.

📐The fundamental mode of a string

The finite-difference matrix of

-u'' = \lambda u

(0,1)

with 99 interior points has smallest eigenvalue

4\sin^2(\pi h/2)/h^2 = 9.8687926854

, close to the continuous

\pi^2 = 9.8696

. Inverse iteration from the all-ones vector gives

9.99900

9.87016

9.86881

, … and has all eleven digits after 8 iterations. It converges even faster than the predicted

1/4

because the symmetric starting vector contains no antisymmetric second mode: the effective ratio is

\lambda_1/\lambda_3 = 1/9

Derivation

The string's exact discrete eigenvalues

The matrix here is lesson 3-3's second-difference matrix divided by

h^2

, and lesson 4-1 found the eigenvalues of

\operatorname{tridiag}(-1, 2, -1)

exactly:

4\sin^2(j\pi h/2)

, with sampled sine waves as eigenvectors. Dividing by

h^2

\lambda_j = \frac{4}{h^2}\sin^2\frac{j\pi h}{2}, \qquad j = 1, \dots, n.

For small

h

\sin x \approx x - x^3/6

, so

\lambda_j \approx j^2\pi^2\big(1 - j^2\pi^2h^2/12\big)

: close to the continuous eigenvalues

j^2\pi^2

-u'' = \lambda u

, slightly low, with an

O(h^2)

error. With

h = 0.01

\lambda_1 = 9.86879

against

\pi^2 = 9.86960

. The ratios

\lambda_1 : \lambda_2 : \lambda_3 \approx 1 : 4 : 9

give the rates

1/4

and

1/9

quoted above.

💡Shifts: the eigenvalue nearest σ

Apply inverse iteration to

A - \sigma I

. Its eigenvalues are

\lambda_i - \sigma

, so the iteration converges to the eigenvector whose eigenvalue is nearest the shift

\sigma

, at rate

\dfrac{|\lambda_{\text{nearest}} - \sigma|}{|\lambda_{\text{next nearest}} - \sigma|}.

A good shift makes this tiny. With

\sigma = 85

, aimed at the string's third eigenvalue

88.7607

, five iterations give

88.760707938

. The closer

\sigma

is, the faster — and when

\sigma

is almost exactly an eigenvalue,

A - \sigma I

is almost singular. The solve is then badly conditioned, but harmlessly so: the large error points almost entirely along the eigenvector we want, and normalisation removes its size.

Rayleigh quotient iteration. Why keep the shift fixed? Use the best current estimate — the Rayleigh quotient — as the shift at every step:

\mu_k = \rho(\mathbf x_k), \qquad (A - \mu_k I)\mathbf y = \mathbf x_k, \qquad \mathbf x_{k+1} = \mathbf y/\|\mathbf y\|.

For symmetric matrices the convergence is cubic: errors went

3.7 \times 10^{-2} \to 3.9 \times 10^{-9} \to 8.5 \times 10^{-14}

in the lab. The price: a new factorisation every step, and less control over which eigenvalue you land on — started from a rough guess it converges to whichever eigenpair it is nearest, which may not be the one you wanted.

Derivation

Why Rayleigh quotient iteration is cubic

Suppose the current vector has error

\varepsilon

(its component along the other eigenvectors). Two facts from this unit combine:

The shift is the Rayleigh quotient, so its error is of order $\varepsilon^2$ (lesson 6-1): $|\lambda - \mu_k| \approx c\,\varepsilon^2$ .
One step of shifted inverse iteration multiplies the unwanted components by $|\lambda - \mu_k|/|\lambda_{\text{next}} - \mu_k|$ — and the denominator stays roughly the gap $g$ to the next eigenvalue.

Together:

\varepsilon_{\text{new}} \approx \varepsilon\,\frac{c\,\varepsilon^2}{g} = \frac{c}{g}\,\varepsilon^3.

Since the eigenvalue error is the square of the vector error, it is cubed too, and the number of correct digits roughly triples each step. In the lab the eigenvalue's relative error goes from

1.5 \times 10^{-4}

(about 4 digits) to

1.6 \times 10^{-11}

(about 11) in one step, and the next step reaches rounding level. For nonsymmetric matrices the Rayleigh quotient is only first-order accurate and the convergence drops to quadratic.

nothing whatsoever takes place in the universe in which some relation of maximum and minimum does not appear.

– Leonhard Euler, De Curvis Elasticis (1744), tr. Oldfather et al. (1933)

💡Deflation: one mode after another

Once mode 1 is found, how do we get mode 2? For symmetric matrices the eigenvectors are orthogonal, so project out the converged vectors at every step:

\mathbf x \leftarrow \mathbf x - \sum_{j}(\mathbf v_j^{\mathsf T}\mathbf x)\mathbf v_j

. Inverse iteration on what remains finds the next-smallest eigenvalue. For a five-storey shear building this gives periods of 0.7805, 0.2674 and 0.1696 s — identical to scipy.linalg.eigh. (With a mass matrix, 'orthogonal' means

M

-orthogonal: lesson 6-4.) Deflation one vector at a time accumulates error and is slow for many modes. Block versions — subspace iteration, and the Lanczos method (Cornelius Lanczos, 1950), which is to eigenvalues what conjugate gradients are to linear systems — compute dozens of modes together. scipy.sparse.linalg.eigsh(K, M=M, sigma=0) is shift-and-invert Lanczos: exactly this lesson, industrialised.

💡Where you will use this

Structures and Earthquake Engineering: fundamental periods and the first few modes for response-spectrum analysis — every commercial FE package uses shift-and-invert Lanczos. Buckling: the smallest buckling load factor is the smallest eigenvalue of a generalised problem. Spacecraft Attitude Control: the lowest flexible-appendage frequency must stay well above the controller bandwidth. Rocket Propulsion: the lowest acoustic modes of a combustion chamber, where instability starts.

📚Further reading

Beresford Parlett, The Symmetric Eigenvalue Problem (Prentice-Hall, 1980; SIAM Classics 1998) — inverse iteration, Rayleigh quotient iteration and Lanczos, with the cubic-convergence proof.
Lloyd N. Trefethen and David Bau, Numerical Linear Algebra (SIAM, 1997), Lecture 27, 'Rayleigh Quotient, Inverse Iteration'.
Klaus-Jürgen Bathe, Finite Element Procedures (Prentice-Hall, 1996), chapter 11, 'Solution Methods for Eigenproblems' — subspace iteration and Lanczos as structural codes use them.

↓But there is a problem

Inverse iteration, shifts and deflation find eigenpairs one at a time, each needing a factorisation, a starting vector and a target. To get all eigenvalues of a dense

n \times n

matrix that way is slow and fragile. In 1959 a 25-year-old programmer in London, John Francis, without a university degree, found an algorithm that finds them all at once — by factorising the matrix and multiplying the factors back together in the wrong order, over and over. It was later named one of the ten most important algorithms of the twentieth century. Francis left the field in 1962, and did not find out for 45 years.