skip to content

How does power iteration find the dominant eigenvector of a large matrix?

level: seniorimportance: nice to knowfreq 30%

answer

  1. just multiply, over and over
  2. the biggest eigenvalue drowns out the rest
  3. renormalize so nothing overflows
  4. error falls by |lambda2 / lambda1| each step
  5. reads the value off the Rayleigh quotient

basics

~20 s

Start from a random nonzero vector, repeatedly multiply by the matrix, and renormalize each step. Components along smaller eigenvalues shrink relative to the largest one, so the vector converges to the dominant eigenvector at rate |lambda2 / lambda1| per iteration.

solid answer

~50 s

Write the start vector in the eigenbasis as `x = c1*v1 + c2*v2 + ...`. Applying `A` repeatedly gives `A^k x = c1*lambda1^k*v1 + c2*lambda2^k*v2 + ...`. Dividing through by `lambda1^k`, every term other than the first carries a factor `(lambda_i / lambda1)^k`, which decays to zero as long as `|lambda1|` is strictly largest. So after normalizing at each step to stop overflow, the iterate lines up with `v1`. Convergence is linear at rate `|lambda2 / lambda1|`: with `lambda1 = 5` and `lambda2 = 2`, the error shrinks by a factor `2/5` per iteration, so about ten steps buy four decimal digits. The eigenvalue is then read off with the Rayleigh quotient `v^T A v / (v^T v)`. Its appeal is that it needs only matrix-vector products, so it scales to huge sparse matrices that could never be factorized.

code

python · 19 lines
python
A = [[4, 1], [2, 3]]              # eigenvalues 5 and 2

def matvec(M, v):
    return [sum(M[i][j] * v[j] for j in range(len(v))) for i in range(len(M))]

def normalize(v):
    norm = sum(x * x for x in v) ** 0.5
    return [x / norm for x in v]

v = normalize([1.0, 0.0])         # arbitrary nonzero start
for k in range(1, 13):
    v = normalize(matvec(A, v))
    rayleigh = sum(a * b for a, b in zip(v, matvec(A, v)))
    print(k, [round(x, 6) for x in v], round(rayleigh, 6))

# k=1  [0.894427, 0.447214]  5.0
# k=4  [0.720636, 0.693313]  5.018197
# k=12 [0.707116, 0.707098]  5.000013
# the direction error shrinks by a factor of 2/5 per step

go deeper

for a junior

Know the loop: multiply by the matrix, renormalize, repeat, and it settles on the direction belonging to the largest eigenvalue in magnitude.

for a middle

Be able to derive the convergence from the eigenbasis expansion, state the rate as the ratio of the two largest eigenvalue magnitudes, and explain why normalization is needed in floating point.

for a senior

Show the operating knowledge: a small spectral gap, not matrix size, drives cost; equal-magnitude top eigenvalues break it; and only matrix-vector products are needed, which is why it survives at scale.

for a principal

Own the build-versus-approximate call. Decide when the top one or two directions genuinely suffice, what accuracy the downstream use actually needs, and how iteration budgets and restarts are set so a slow-converging job cannot stall a pipeline.

## The algorithm Power iteration is three lines: 1. Pick a random nonzero start vector `x0`. 2. Repeat: `x <- A x`, then `x <- x / ||x||`. 3. Estimate the eigenvalue as the Rayleigh quotient `lambda ~= (x^T A x) / (x^T x)`, which is just `x^T A x` when `x` is a unit vector. The normalization step does no mathematical work — it does not change the direction — but it is essential in floating point, because without it the entries grow like `lambda1^k` and overflow (or, when `|lambda1| < 1`, underflow to zero). ## Why it converges Assume `A` has an eigenbasis with eigenvalues ordered `|lambda1| > |lambda2| >= ... >= |lambda_n|`. Expand the start vector in that basis: ``` x0 = c1*v1 + c2*v2 + ... + cn*vn ``` Applying `A` scales each component by its own eigenvalue, so after `k` steps: ``` A^k x0 = c1*lambda1^k*v1 + c2*lambda2^k*v2 + ... = lambda1^k * [ c1*v1 + c2*(lambda2/lambda1)^k*v2 + ... ] ``` Every ratio `|lambda_i / lambda1|` is strictly less than 1, so every term but the first decays geometrically. The bracket converges to `c1*v1`, and after normalization the iterate converges to the unit dominant eigenvector (up to sign, which may flip each step when `lambda1` is negative). The overall factor `lambda1^k` is exactly what the normalization strips away. ## Convergence rate The error decays like `|lambda2 / lambda1|^k`, which is **linear** (geometric) convergence: a constant factor per iteration, not a doubling of accurate digits. With `lambda1 = 5` and `lambda2 = 2` the factor is `0.4`, so each iteration removes 60% of the remaining error and roughly every 2.5 iterations buys one decimal digit. The practical consequence is stark: when the top two eigenvalues are close in magnitude — a ratio of `0.99`, say — the method crawls, needing hundreds of iterations for a few digits. The spectral gap, not the matrix size, sets the cost. A useful asymmetry: the eigenvector converges at rate `r = |lambda2/lambda1|`, while for a symmetric matrix the Rayleigh-quotient eigenvalue estimate converges at rate `r^2`, because the quotient is stationary at an eigenvector. The eigenvalue is usually accurate well before the direction is. ## When it fails or misbehaves - **No strictly dominant eigenvalue.** If `|lambda1| = |lambda2|` — for example eigenvalues `+3` and `-3`, or a complex conjugate pair — the two components never decay relative to each other and the iterate oscillates instead of settling. - **A start vector with no dominant component.** If `c1 = 0`, the theory converges to the second eigenvector instead. With a random start this has probability zero, but a structured start (all ones, a basis vector) can hit it in a symmetric problem; rounding error usually reintroduces a small `c1` and rescues it slowly. - **A tiny spectral gap.** Not a failure, just intolerably slow; this is the common real-world complaint. - **Sign flapping.** With `lambda1 < 0` the iterate alternates direction each step. The axis is converging even though the arrow keeps flipping. ## Why anyone uses it The method touches the matrix only through the product `A x`. That means it never needs `A` stored densely or factorized: a sparse matrix, or even a rule that computes `A x` implicitly, is enough. Memory is one vector. For a matrix with millions of rows where a full decomposition is impossible and only the leading direction is wanted, this is decisive. ## Extensions worth naming - **Deflation.** After converging to `v1`, subtract its contribution (for a symmetric matrix, work with `A - lambda1*v1*v1^T`) and iterate again to get the next eigenpair. - **Block or subspace iteration.** Iterate several orthonormal vectors at once, re-orthogonalizing each step, to get the top `k` directions together. - **Shifting.** The eigenvalues of `A - s*I` are `lambda_i - s` with the same eigenvectors, so a well-chosen shift can enlarge the effective gap and speed convergence dramatically. - **Inverse iteration.** Applying power iteration to `(A - s*I)^-1` converges to the eigenvector whose eigenvalue is *closest* to `s`, which is how a specific interior eigenvector is targeted. All of these keep the same core insight: repeated multiplication amplifies whichever eigen-direction has the largest magnitude, and everything else is arranging for that direction to be the one you want.

  • When does power iteration fail to converge to a single eigenvector?
    When there is no strictly dominant eigenvalue in magnitude. Eigenvalues `+3` and `-3`, or a complex conjugate pair, leave two components that never shrink relative to each other, so the iterate oscillates forever. It also stalls, without failing, when the top two magnitudes are merely close, since the error factor is their ratio.
  • Why normalize the vector on every iteration rather than at the end?
    Because the entries scale by `lambda1` each step, so in floating point they overflow when `|lambda1| > 1` and underflow to zero when it is below 1. Normalizing costs nothing mathematically, since scaling does not change the direction, and it keeps every iterate in a numerically safe range.
  • How would you get the second eigenvector after finding the first?
    Deflate: for a symmetric matrix, iterate again on `A - lambda1*v1*v1^T`, which zeroes out the dominant eigenvalue while leaving all other eigenpairs untouched, so the same procedure now converges to the second. Alternatively iterate a block of orthonormal vectors together, re-orthogonalizing each step, to obtain several leading directions at once.
  • Why is power iteration attractive for very large matrices?
    Because it touches the matrix only through the product `A x`. Nothing is factorized and nothing dense is stored, so a sparse matrix or even an implicit rule for computing `A x` suffices, with memory for a single vector. When only the leading direction of a matrix with millions of rows is needed, a full decomposition is not an option and this is.

Play a chord through an amplifier that boosts one frequency far more than the others, feed the output back in, and after enough passes you hear only that one note. Power iteration amplifies the dominant eigen-direction the same way.

saying these in an interview costs you the question

  • Claims convergence is quadratic rather than linear
  • Says the rate depends on matrix size rather than the eigenvalue gap
  • Omits normalization and ignores overflow
  • Thinks it returns all eigenvectors at once
  • Believes it works when the top two eigenvalues have equal magnitude
  • Confuses the Rayleigh quotient with the vector norm

context

resolving links…