BrevFeed
Sign in

41. Iterative methods (intro): power iteration for dominant eigenvalue

📖 Reading · 10 min
💡 Every code box below is live — edit it and hit Run.

Every eigenvalue method so far — the characteristic polynomial (§19.0), Schur's QR algorithm (§22.1) — computes all of a matrix's eigenvalues at once, and needs genuine linear-algebra machinery to do it. This closing numerical-methods lesson asks a narrower question with a startlingly simple answer: what if only the largest eigenvalue is needed? Repeated matrix-vector multiplication, nothing else, finds it.

The algorithm

Given A with a strictly dominant eigenvalue \lambda_1 (meaning |\lambda_1|>|\lambda_2|\ge\cdots\ge|\lambda_n| — a genuine gap, not a tie), start from almost any vector \vec x_0 and repeat:

\vec x_{k+1}=\frac{A\vec x_k}{\|A\vec x_k\|}

(normalizing at each step, purely to prevent the vector's length from exploding or vanishing — it doesn't change the direction the iteration converges to). As k\to\infty, \vec x_k converges to the eigenvector for \lambda_1, and the Rayleigh quotient \vec x_k^TA\vec x_k (with \vec x_k unit length) converges to \lambda_1 itself.

Why it works

Assume A is diagonalizable with eigenbasis \vec v_1,\dots,\vec v_n (§19.1) and write \vec x_0=c_1\vec v_1+\cdots+c_n\vec v_n (possible since it's a basis). Then, before normalizing,

A^k\vec x_0=c_1\lambda_1^k\vec v_1+c_2\lambda_2^k\vec v_2+\cdots+c_n\lambda_n^k\vec v_n=\lambda_1^k\left[c_1\vec v_1+c_2\left(\frac{\lambda_2}{\lambda_1}\right)^k\vec v_2+\cdots\right]

(factoring \lambda_1^k out of every term). Since \lambda_1 is strictly dominant, |\lambda_i/\lambda_1|<1 for every i\ge2, so every term but the first shrinks to zero as k grows — normalizing away the runaway \lambda_1^k scale factor at each step leaves the vector converging to (a multiple of) \vec v_1 alone. This is the exact mechanism behind §19.1's power formula, A^k=PD^kP^{-1}, run in reverse: instead of using known eigenvalues to compute a power, the power itself is used to reveal the largest eigenvalue.

Convergence rate, and when it's slow

The rate is governed by |\lambda_2/\lambda_1| — the second-largest eigenvalue ratio. If it's close to 1 (a near-tie for dominance), convergence crawls; if \lambda_1 is well-separated from the rest, it's fast. If |\lambda_1|=|\lambda_2| exactly (no strict dominance — e.g. a rotation matrix's complex-conjugate pair, §18.2/§21.5), power iteration doesn't converge to a single direction at all; the method requires the gap assumption stated at the outset.

Doing it in Python

import numpy as np

A = np.array([[4., 1.], [2., 3.]])   # eigenvalues 5, 2 -- a real gap

x = np.array([1., 0.])
for i in range(15):
    x = A @ x
    x = x / np.linalg.norm(x)

rayleigh = x @ A @ x
print("converged eigenvector:", [round(v, 4) for v in x.tolist()])
print("Rayleigh quotient (eigenvalue estimate):", round(float(rayleigh), 4))

true_eigvals = np.linalg.eigvals(A)
print("true eigenvalues:", sorted([round(v, 4) for v in true_eigvals], reverse=True))

converged eigenvector: [0.7071, 0.7071]
Rayleigh quotient (eigenvalue estimate): 5.0
true eigenvalues: [5.0, 2.0]

Watching convergence speed depend on the eigenvalue gap directly:

import numpy as np

def power_iteration(A, iters):
    x = np.array([1., 1.])   # NOT already an eigenvector -- both directions present
    estimates = []
    for _ in range(iters):
        x = A @ x
        x = x / np.linalg.norm(x)
        estimates.append(float(x @ A @ x))
    return estimates

A_fast = np.array([[10., 0.], [0., 1.]])   # ratio 1/10 -- fast
A_slow = np.array([[1.1, 0.], [0., 1.]])   # ratio 1/1.1 -- slow

for name, A in [("A_fast", A_fast), ("A_slow", A_slow)]:
    est = power_iteration(A, 5)
    print(f"{name}: estimates after 5 steps = {[round(e, 4) for e in est]}")

A_fast: estimates after 5 steps = [9.9109, 9.9991, 10.0, 10.0, 10.0]
A_slow: estimates after 5 steps = [1.0548, 1.0594, 1.0639, 1.0682, 1.0722]

Worked example

Estimate the dominant eigenvalue of A=\begin{pmatrix}2&1\\1&2\end{pmatrix} using two steps of power iteration from \vec x_0=(1,0).

\vec x_1=A\vec x_0=(2,1), \|\vec x_1\|=\sqrt5, normalized: (2/\sqrt5,1/\sqrt5)\approx(0.8944,0.4472).

\vec x_2=A\vec x_1'\approx(2(0.8944)+0.4472,\,0.8944+2(0.4472))=(2.2361,1.7889), \|\vec x_2\|\approx2.8636, normalized \approx(0.7810,0.6247).

Rayleigh quotient: \vec x_2'^TA\vec x_2'\approx0.7810(2(0.7810)+0.6247)+0.6247(0.7810+2(0.6247)) \approx0.7810(2.1867)+0.6247(2.0303)\approx1.7078+1.2684\approx2.976.

\boxed{\lambda_1\approx2.976\text{ after two steps}}

Sanity check. The true eigenvalues (§19.0's trace/determinant shortcut: \operatorname{tr}=4,\det=3, \lambda^2-4\lambda+3=0) are \lambda=1,3 — the estimate 2.976 after just two steps is already close to the true dominant eigenvalue 3, consistent with a reasonably fast convergence rate (|\lambda_2/\lambda_1|=1/3, not a near-tie).

Your turn

1. Why does power iteration fail to find the dominant eigenvalue if started from \vec x_0=\vec v_2 exactly (an eigenvector for a non-dominant eigenvalue)?

2. For A=\operatorname{diag}(5,5,1), does power iteration converge to a single eigenvector? Why or why not?

3. True or false: power iteration requires computing a characteristic polynomial at any point.

Solutions

1. If \vec x_0=\vec v_2 exactly, then c_1=0 in the eigenbasis expansion — there's no \vec v_1 component to survive at all, and A^k\vec x_0=\lambda_2^k\vec v_2 forever, converging (trivially) to \vec v_2 instead. In practice this is rarely an issue: floating-point rounding almost always introduces a tiny nonzero c_1 component, which the same amplification mechanism from "Why it works" above eventually grows dominant anyway — but it is the one genuine failure mode of the method.

2. No, not to a single eigenvector. |\lambda_1|=|\lambda_2|=5 — not strictly dominant (a tie), violating the convergence assumption stated at the top. The iteration keeps a mixture of the two \lambda=5 eigenvectors' directions rather than settling on one particular vector (any vector in that eigenspace is equally "dominant," so there is no single direction being selected for).

3. False. This is the entire point of the method: it uses only repeated matrix-vector multiplication and normalization — no determinant, no characteristic polynomial, no elimination at all. This is exactly why it scales to enormous matrices (search-engine-scale graphs, §24.2's spectral graph theory) where forming a characteristic polynomial would be computationally hopeless.

Check yourself in code

Run 10 steps of power iteration on A=\begin{pmatrix}6&2\\2&3\end{pmatrix} from \vec x_0=(1,0), and report the final Rayleigh quotient (eigenvalue estimate).

Print exactly this:

eigenvalue estimate: 7.0

import numpy as np

A = np.array([[6., 2.], [2., 3.]])
x = np.array([1., 0.])
for _ in range(10):
    x = A @ x
    x = x / np.linalg.norm(x)

# print "eigenvalue estimate: <rayleigh quotient, rounded to 4 dp>"

import numpy as np

A = np.array([[6., 2.], [2., 3.]])
x = np.array([1., 0.])
for _ in range(10):
    x = A @ x
    x = x / np.linalg.norm(x)

rayleigh = x @ A @ x
print("eigenvalue estimate:", round(float(rayleigh), 4))

Power iteration finds a matrix's dominant eigenvalue and eigenvector using nothing but repeated matrix-vector multiplication and normalization, converging because the dominant eigenvalue's direction grows fastest under repeated application — the exact mechanism behind §19.1's power formula, run to reveal an unknown eigenvalue instead of use a known one. It needs a strict dominance gap, and converges faster when that gap is larger.

Next, closing Module 22 with the theorem that guarantees power iteration's setup — a strictly dominant real, positive eigenvalue — actually holds for an important class of matrices: Perron-Frobenius.

← Previous Next →