🌀Eigenvalue Problems

Eigenvalues and the Power Method

Eigenpairs and their physical meaning, power iteration, the Rayleigh quotient, PageRank and Gershgorin discs

Puget Sound, 7 November 1940

The Tacoma Narrows Bridge has been open for four months. Its main span, 853 m long and remarkably slender, has always moved in the wind — drivers call it Galloping Gertie — but the motion has been an up-and-down heave that the engineers expected to damp out. That morning, in a steady wind of about 64 km/h, the motion changes. A cable band slips, and the deck begins to twist: one half rotating one way, the other half the opposite way, with the centre line of the span almost still. The twisting grows with every cycle. Shortly after eleven o'clock, the deck tears apart and falls into the water. Nobody is killed. The board of engineers that investigated — Theodore von Kármán among them — concluded that the wind had fed energy into one particular shape of motion, at its own frequency, faster than the structure could dissipate it. Not a resonance with gusts: a self-excited torsional oscillation that modern engineers call flutter. A structure has a handful of such shapes and frequencies, and every load it ever feels can be written as a mixture of them. In the language of linear algebra they are eigenvectors and eigenvalues.

A nonzero vector

\mathbf v

is an eigenvector of a square matrix

A

, with eigenvalue

\lambda

, if

A\mathbf v = \lambda\mathbf v.

A

does not turn

\mathbf v

; it only stretches it by

\lambda

. What that means depends on what

A

is:

for a stiffness matrix (with unit masses), eigenvectors are mode shapes and $\lambda = \omega^2$ gives the natural frequencies;
for a state matrix in $\dot{\mathbf x} = A\mathbf x$ , each eigenvalue is a growth rate: the system is stable when every $\operatorname{Re}\lambda < 0$ , and $\operatorname{Im}\lambda$ is a frequency;
for a stress tensor, eigenvalues are principal stresses and eigenvectors principal directions; for an inertia tensor, principal moments and principal axes of a spacecraft;
for a transition matrix of a Markov chain, the eigenvector for $\lambda = 1$ is the long-run distribution.

Eigenvalues are roots of

\det(A - \lambda I) = 0

, but lesson 1-4's Wilkinson polynomial showed why we never compute them that way.

Derivation

Why eigenvalues are frequencies and growth rates

Growth rates. Look for a solution of

\dot{\mathbf x} = A\mathbf x

that keeps its shape and only changes in size:

\mathbf x(t) = \mathbf v\,e^{\lambda t}

. Then

\dot{\mathbf x} = \lambda\mathbf v\,e^{\lambda t}

and

A\mathbf x = A\mathbf v\,e^{\lambda t}

; cancelling the exponential leaves

A\mathbf v = \lambda\mathbf v.

So every eigenpair is a solution, and any combination of them is too. Writing

\lambda = \sigma + i\omega

e^{\lambda t} = e^{\sigma t}(\cos\omega t + i\sin\omega t)

: the size changes like

e^{\sigma t}

and the phase turns at rate

\omega

. Every mode decays exactly when every

\operatorname{Re}\lambda = \sigma < 0

Natural frequencies. For undamped vibration with unit masses,

\ddot{\mathbf x} + K\mathbf x = \mathbf 0

. Try

\mathbf x = \boldsymbol\phi\cos\omega t

; then

\ddot{\mathbf x} = -\omega^2\boldsymbol\phi\cos\omega t

and

K\boldsymbol\phi = \omega^2\boldsymbol\phi:

the eigenvalues of the stiffness matrix are the squared natural frequencies, the eigenvectors the mode shapes. Lesson 6-4 puts the mass matrix back in.

💡The power method

Start from almost any vector and multiply by

A

repeatedly, normalising each time:

\mathbf x_{k+1} = \dfrac{A\mathbf x_k}{\|A\mathbf x_k\|}.

Write

\mathbf x_0 = c_1\mathbf v_1 + c_2\mathbf v_2 + \dots

in eigenvectors with

|\lambda_1| > |\lambda_2| \ge \dots

. Then

A^k\mathbf x_0 = c_1\lambda_1^k\mathbf v_1 + c_2\lambda_2^k\mathbf v_2 + \dots

, and every other component shrinks relative to the first by

|\lambda_2/\lambda_1|^k

. The iterate turns towards the dominant eigenvector, linearly, at rate

|\lambda_2/\lambda_1|

. Each step costs one matrix–vector product — nothing else — which is why it works on matrices with billions of rows.

\rho(\mathbf x) = \frac{\mathbf x^{\mathsf T}A\mathbf x}{\mathbf x^{\mathsf T}\mathbf x}

The Rayleigh quotient: the best eigenvalue estimate for a given approximate eigenvector. For symmetric A, an eigenvector error ε gives an eigenvalue error of order ε², so the estimate converges twice as fast as the vector.

Derivation

Why the Rayleigh quotient's error is squared

Let

\mathbf v

be a unit eigenvector,

A\mathbf v = \lambda\mathbf v

, and let the approximation be

\mathbf x = \mathbf v + \varepsilon\mathbf w

, with

\mathbf w

a unit vector perpendicular to

\mathbf v

and

\varepsilon

small. For symmetric

A

the cross terms vanish:

\mathbf v^{\mathsf T}A\mathbf w = (A\mathbf v)^{\mathsf T}\mathbf w = \lambda\,\mathbf v^{\mathsf T}\mathbf w = 0

. So

\rho(\mathbf x) = \frac{(\mathbf v + \varepsilon\mathbf w)^{\mathsf T}A(\mathbf v + \varepsilon\mathbf w)}{(\mathbf v + \varepsilon\mathbf w)^{\mathsf T}(\mathbf v + \varepsilon\mathbf w)} = \frac{\lambda + \varepsilon^2\,\mathbf w^{\mathsf T}A\mathbf w}{1 + \varepsilon^2} = \lambda + \varepsilon^2\,\frac{\mathbf w^{\mathsf T}A\mathbf w - \lambda}{1 + \varepsilon^2}.

The error in the eigenvalue is proportional to

\varepsilon^2

: no term in

\varepsilon

survives. A vector good to 4 digits gives an eigenvalue good to about 8. For the spring chain, the vector error shrinks by

0.786

per step, so the Rayleigh quotient's shrinks by

0.786^2 = 0.62

— the factor

0.008

every 10 steps quoted below.

📐A chain of springs

For the symmetric

4 \times 4

stiffness matrix with diagonal

(4, 4, 4, 3)

and off-diagonals

-1

, the eigenvalues are

2.1206,\ 3,\ 4.3473,\ 5.5321

. The ratio

|\lambda_2/\lambda_1| = 0.786

is close to 1, so the power method is slow: the Rayleigh quotient has error

4.9 \times 10^{-2}

after 10 steps,

4.1 \times 10^{-4}

after 20,

2.7 \times 10^{-8}

after 40 — each 10 steps buy a factor

0.786^{20} \approx 0.008

, the square of the vector's

0.786^{10}

, exactly as the Rayleigh-quotient caption predicts.

A scientist studies what is, whereas an engineer creates what never was.

– Theodore von Kármán, as quoted in Theoretical Foundations for Decision Making in Engineering Design (2001)

🛰The $25,000,000,000 eigenvector

In 1998 Sergey Brin and Lawrence Page, graduate students at Stanford, described a search engine that ranked web pages by importance: a page is important if important pages link to it. Model a surfer who, with probability

d = 0.85

, follows a random link from the current page, and otherwise jumps to any page at random. The long-run fraction of time spent on each page — PageRank — is the dominant eigenvector, eigenvalue 1, of the 'Google matrix'

G = dP + (1-d)\tfrac1n\mathbf 1\mathbf 1^{\mathsf T}

. The random jump guarantees

|\lambda_2| \le d

, so the power method gains a digit about every 14 iterations, however many billions of pages there are. A 2006 SIAM Review article about it carried the title The $25,000,000,000 Eigenvector.

💡Gershgorin discs

Every eigenvalue of

A

lies in at least one disc in the complex plane centred at a diagonal entry

a_{ii}

with radius

R_i = \sum_{j \ne i}|a_{ij}|

(Semyon Gershgorin, 1931). Proof in one line: if

A\mathbf v = \lambda\mathbf v

and

v_i

is the largest component, row

i

gives

|\lambda - a_{ii}| \le \sum_{j\ne i}|a_{ij}||v_j|/|v_i| \le R_i

. It costs nothing and is often enough. For a linearised system

\dot{\mathbf x} = A\mathbf x

whose discs all lie in the left half-plane, stability is proved without computing a single eigenvalue. A strictly diagonally dominant matrix (lesson 4-1) has no disc containing 0, so it is nonsingular.

💡Where you will use this

Spacecraft Attitude Control: principal axes of the inertia tensor; the eigenvalues of the linearised attitude dynamics decide whether spin about an axis is stable. Control Systems: closed-loop poles are eigenvalues; pole placement moves them. Structures: buckling loads are eigenvalues of

K - P K_G

. Geotechnics: principal stresses for failure criteria. Markov models: reliability and degradation chains have their steady state as an eigenvector.

📚Further reading

Lloyd N. Trefethen and David Bau, Numerical Linear Algebra (SIAM, 1997), Lectures 24 and 27, 'Eigenvalue Problems' and 'Rayleigh Quotient, Inverse Iteration'.
Michael Heath, Scientific Computing: An Introductory Survey (2nd ed.; SIAM Classics 2018), chapter 4, 'Eigenvalue Problems'.
Kurt Bryan and Tanya Leise, 'The $25,000,000,000 eigenvector: the linear algebra behind Google', SIAM Review 48(3), 569–581 (2006).

↓But there is a problem

The power method finds the eigenvalue of largest magnitude. For a structure, that is the highest natural frequency of the model — usually a meaningless jiggle of the smallest elements in the mesh. Engineers almost always need the other end. Every building code in the world asks first for the fundamental period of a structure — the period of its slowest mode, the smallest eigenvalue — because that is what an earthquake shakes hardest. The power method never sees it.