Iterative Estimates for Eigenvalues

In applications, eigenvalues are rarely needed to full precision, and the characteristic polynomial is the wrong tool for computing them: no finite formula solves a general polynomial of degree $5$ or more, and forming then rooting the polynomial is numerically unstable. Iterative methods avoid it, producing accurate approximations from repeated matrix-vector products. The power method and its relatives below also motivate the numerical eigenvalue algorithms.

The dominant eigenvalue

Assume $A$ is diagonalizable with a basis of eigenvectors $v_{1}, \dots, v_{n}$ and a strictly dominant eigenvalue, meaning

∣ λ_{1} ∣ > ∣ λ_{2} ∣ \geq ∣ λ_{3} ∣ \geq \dots \geq ∣ λ_{n} ∣.

Write a starting vector as $x = c_{1} v_{1} + \dots + c_{n} v_{n}$ with $c_{1} \neq = 0$ . Then

A^{k} x = c_{1} (λ_{1})^{k} v_{1} + c_{2} (λ_{2})^{k} v_{2} + \dots + c_{n} (λ_{n})^{k} v_{n} .

Dividing by $(λ_{1})^{k}$ isolates the dominant term:

\frac{1}{( λ _{1} ) ^{k}} A^{k} x = c_{1} v_{1} + c_{2} (\frac{λ _{2}}{λ _{1}})^{k} v_{2} + \dots + c_{n} (\frac{λ _{n}}{λ _{1}})^{k} v_{n} .

Every ratio $λ_{j} / λ_{1}$ has magnitude below $1$ , so its powers vanish. Hence

\frac{1}{( λ _{1} ) ^{k}} A^{k} x ⟶ c_{1} v_{1} as k \to \infty.

The vectors $A^{k} x$ line up with the dominant eigenvector: the direction of $A^{k} x$ converges to the eigenspace of $λ_{1}$ .

Successive iterates

x, A x, A^{2} x, \dots

swing toward the dominant eigenspace; each multiplication by

A

suppresses the subdominant component.

The algorithm

The scaling factor $(λ_{1})^{k}$ is unknown, so instead each iterate is rescaled so its largest entry is $1$ . The largest entry of $A x_{k}$ then converges to $λ_{1}$ , because $A x_{k} \approx λ_{1} x_{k}$ once $x_{k}$ is close to an eigenvector.

Algorithm:

\textsc{PowerMethod}(A, x_0)

— estimate the dominant eigenvalue

1
require x_0 has largest entry equal to 1
2
for k = 0, 1, 2, ... do
3
y <- A x_k
4
mu <- an entry of y with the largest absolute value
eigenvalue estimate
5
x_{k+1} <- (1 / mu) y
rescale largest entry to 1
6
until mu and x_k stop changing
7
return (mu, x_k)
dominant eigenvalue and its eigenvector

$k$	$x_{k}$	$A x_{k}$	$μ_{k}$
$0$	$(0, 1)$	$(5, 2)$	$5$
$1$	$(1, 0.4)$	$(8, 1.8)$	$8$
$2$	$(1, 0.225)$	$(7.125, 1.45)$	$7.125$
$3$	$(1, 0.2035)$	$(7.0175, 1.407)$	$7.0175$
$4$	$(1, 0.2005)$	$(7.0025, 1.401)$	$7.0025$
$5$	$(1, 0.20007)$	$(7.00036, 1.40014)$	$7.00036$

The estimates $μ_{k} \to 7$ and $x_{k} \to (1, 0.2)$ , and indeed $A (1, 0.2)^{⊤} = (7, 1.4)^{⊤} = 7 (1, 0.2)^{⊤}$ : the dominant eigenvalue is $7$ with eigenvector $(1, 0.2)$ .

Rate of convergence

The leftover error in using a scaled $A^{k} x$ for $c_{1} v_{1}$ is dominated by the term $c_{2} (λ_{2} / λ_{1})^{k} v_{2}$ , so the convergence is geometric with ratio $∣ λ_{2} / λ_{1} ∣$ .

In the worked run $λ_{2} = 1$ and $λ_{1} = 7$ , so the ratio is $1/7$ and the estimates gained roughly one digit per step.

Error against iteration for two ratios; a small

∣ λ_{2} / λ_{1} ∣ = 1/7

(steep line) converges far faster than a ratio near

0.95

(shallow line).

When $∣ λ_{2} / λ_{1} ∣$ is near $1$ , the plain power method is too slow and other approaches are preferred.

The Rayleigh quotient

The largest-entry estimate $μ_{k}$ is crude. A better eigenvalue estimate from an approximate eigenvector uses the Rayleigh quotient. If $A x = λ x$ , then $x^{⊤} A x = λ (x^{⊤} x)$ , so

R (x) = \frac{x ^{⊤} A x}{x ^{⊤} x}

equals $λ$ exactly for an eigenvector, and is close to $λ$ when $x$ is close to an eigenvector.

From the same iterate the Rayleigh quotient carries about twice the correct digits of the largest-entry estimate for a symmetric matrix.

This doubling is why symmetric problems — the setting of the Spectral Theorem — are especially amenable to iteration.

The inverse power method

The power method finds only the dominant eigenvalue. To target any eigenvalue $λ$ , given a rough estimate $α$ near it, apply the power method to

B = (A - α I)^{- 1} .

If the eigenvalues of $A$ are $λ_{1}, \dots, λ_{n}$ , the eigenvalues of $B$ are $1/ (λ_{i} - α)$ , with the same eigenvectors. When $α$ is close to a particular $λ_{i}$ , the value $1/ (λ_{i} - α)$ is huge — strictly dominant in $B$ — so the power method on $B$ homes in on that eigenvector rapidly.

Shifting by

α

near an eigenvalue makes

1/ (λ_{i} - α)

dominate; the reciprocal blows up the eigenvalue nearest the shift.

In practice $(A - α I)^{- 1}$ is never formed. Instead each step solves a linear system, and the eigenvalue of $A$ is recovered from the estimate for $B$ .

Algorithm:

\textsc{InversePowerMethod}(A, \alpha, x_0)

— estimate an eigenvalue near

\alpha

1
require x_0 has largest entry equal to 1
2
for k = 0, 1, 2, ... do
3
solve (A - alpha I) y = x_k for y
reuse one LU factorization each step
4
mu <- an entry of y with the largest absolute value
5
lambda <- alpha + 1 / mu
eigenvalue estimate for A
6
x_{k+1} <- (1 / mu) y
7
until lambda and x_k stop changing
8
return (lambda, x_k)

Because the same matrix $A - α I$ is solved against every step, one LU factorization is computed once and reused. If no estimate of the smallest eigenvalue is available, taking $α = 0$ (so $B = A^{- 1}$ ) targets the eigenvalue nearest zero.

The methods in context

Method	Finds	Needs	Convergence rate
Power method	dominant eigenvalue $λ_{1}$	a strictly dominant eigenvalue	$∣ λ_{2} / λ_{1} ∣$
Rayleigh quotient	sharper $λ$ from an approximate eigenvector	(best for symmetric $A$ )	doubles the digits
Inverse power method	eigenvalue nearest $α$	a good initial guess $α$	$∣ (λ_{i} - α) / (λ_{j} - α) ∣$

These algorithms handle simple situations well. Production software instead uses the QR algorithm, which generates a sequence of matrices all similar to $A$ that converge to triangular form, reading off every eigenvalue at once. It underlies the standard eigenvalue routines and is developed with the rest of the numerical eigenvalue problem.¹

Lay, Linear Algebra and Its Applications, 5th ed., §5.8 — Iterative Estimates for Eigenvalues: the power method for a strictly dominant eigenvalue and its $∣ λ_{2} / λ_{1} ∣$ convergence rate, the Rayleigh quotient $R (x) = x^{⊤} A x / x^{⊤} x$ , and the inverse power method using $(A - α I)^{- 1}$ to target an eigenvalue near $α$ . ↩