Mathematical foundations of Power Iteration.
Power Iteration—developed in 1929 by the Austrian mathematician Richard von Mises and Hilda Geiringer—is one of the oldest, simplest, and most widely deployed algorithms in numerical linear algebra. In modern data science, web search, quantum mechanics, and structural dynamics, matrices can span millions or billions of dimensions. Calculating all eigenvalues using full factorizations like the QR algorithm requires O(n³) operations and O(n²) memory, which is completely impossible for massive graphs. However, in countless real-world applications, one does not need all n eigenvalues: one needs only the single dominant eigenvalue λ₁ (having the largest absolute magnitude |λ₁| > |λ₂|) and its corresponding principal eigenvector v₁. Power iteration solves this with breathtaking simplicity: starting from an arbitrary unit vector x^{(0)}, one repeatedly multiplies by the matrix A and normalizes: x^{(k+1)} = Ax^{(k)} / ‖Ax^{(k)}‖₂. Expanding the starting vector in the eigenbasis reveals that multiplying by A scales each eigenvector component v_i by λ_i^k. Because |λ₁| is strictly larger than all other eigenvalues, the component along v₁ dominates all others by an exponential factor of (λ₁/λ_i)^k. Consequently, the vector x^{(k)} rotates monotonically toward alignment with the dominant eigenvector v₁. The convergence speed is governed entirely by the spectral gap ratio |λ₂/λ₁|: the error attenuates as O(|λ₂/λ₁|^k). Once the eigenvector direction is established, the Rayleigh quotient μ_k = (x^{(k)})ᵀAx^{(k)} extracts the eigenvalue with doubled precision: the error in μ_k scales quadratically as O(|λ₂/λ₁|^{2k}). The most famous modern deployment of Power Iteration is Google's PageRank algorithm, developed in 1998 by Larry Page and Sergey Brin to rank web pages by computing the principal eigenvector of the global web transition matrix. In this computational laboratory, students evaluate Power Iteration applied to a symmetric matrix with eigenvalues λ = [6.0, 4.0, 1.0] (spectral ratio |λ₂/λ₁| = 4/6 = 0.667, dominant ratio |λ₁/λ₂| = 1.500), observe eigenvector alignment on the unit sphere, compute Rayleigh quotients, and extract the dominant ratio |λ₁/λ₂| = 1.500 via linear regression.
Viva voce questions & analytical theory
What is the mathematical proof of convergence for Power Iteration?
Let A have a dominant eigenvalue |λ₁| > |λ₂| ≥ ... ≥ |λ_n| with corresponding orthonormal eigenvectors v₁, ..., v_n. Expanding the initial vector x^{(0)} in the eigenbasis: x^{(0)} = c₁v₁ + c₂v₂ + ... + c_nv_n with c₁ ≠ 0. Multiplying by A^k yields: A^k x^{(0)} = c₁λ₁^k v₁ + c₂λ₂^k v₂ + ... + c_nλ_n^k v_n = c₁λ₁^k [v₁ + ∑_{i=2}^n (c_i/c₁)(λ_i/λ₁)^k v_i]. Because |λ₁| > |λ_i|, the ratios (λ_i/λ₁)^k vanish exponentially as k → ∞. Thus, normalized iterates x^{(k)} converge directly to the dominant eigenvector v₁, with error dominated by the largest subdominant ratio: ‖x^{(k)} − v₁‖ = O(|λ₂/λ₁|^k).
What is the Rayleigh Quotient and why does it converge twice as fast as the eigenvector?
For any vector x, the Rayleigh quotient is defined as r(x) = (xᵀAx) / (xᵀx). If x is an approximation to an eigenvector v₁ with error ‖x − v₁‖ = O(ε), then Taylor expansion shows that r(x) − λ₁ = O(ε²). For symmetric matrices, the Rayleigh quotient doubles the number of correct digits: while the eigenvector x^{(k)} converges linearly with rate |λ₂/λ₁|, the Rayleigh quotient eigenvalue estimate μ_k converges with rate |λ₂/λ₁|², achieving quadratic accuracy relative to the vector error.
What causes Power Iteration to fail or converge slowly?
Power iteration fails if: (1) The initial vector has zero overlap with the dominant eigenvector (c₁ = 0), although round-off error typically introduces a tiny component that eventually grows; (2) The dominant eigenvalues are complex conjugates or opposite in sign (|λ₁| = |λ₂|), causing iterates to oscillate without settling; or (3) The spectral gap is narrow (|λ₂/λ₁| ≈ 1.0, e.g. 0.999), causing excruciatingly slow convergence requiring thousands of iterations.
How does Google's original PageRank algorithm utilize Power Iteration?
In 1998, Larry Page and Sergey Brin formulated the web link graph as a massive Markov transition matrix M. The PageRank score of web pages is the stationary distribution vector π satisfying Mπ = π, which is precisely the principal eigenvector of M corresponding to eigenvalue λ₁ = 1. To guarantee a unique dominant eigenvalue and prevent dead ends, they introduced a damping factor d = 0.85: G = dM + (1 − d)E. Google computed PageRank across billions of web pages using distributed Power Iteration, where each step requires only one sparse matrix-vector multiplication.
What is Wielandt Deflation and how does it extract subsequent eigenvalues?
Once the dominant eigenpair (λ₁, v₁) is computed, Wielandt deflation constructs the deflated matrix A₁ = A − λ₁ v₁ v₁ᵀ (for symmetric A). The eigenvalues of A₁ are {0, λ₂, λ₃, ..., λ_n}, with eigenvectors identical to A. Running Power Iteration on A₁ then converges to λ₂, the second largest eigenvalue of A. This process can be repeated to extract several leading eigenvalues.
How does Inverse Power Iteration overcome the slow convergence of standard Power Iteration?
Inverse iteration applies power iteration to (A − σI)⁻¹ rather than A, where σ is a scalar shift close to a target eigenvalue λ_j. The eigenvalues of (A − σI)⁻¹ are 1/(λ_i − σ). If σ is very close to λ_j, the eigenvalue 1/(λ_j − σ) becomes enormous compared to all other eigenvalues, creating a vast spectral gap that produces lightning-fast convergence in 2 to 4 steps.
Why is vector normalization required at every single iteration?
If |λ₁| > 1, the unnormalized vector A^k x^{(0)} grows exponentially, causing IEEE 754 floating-point overflow (+Inf) within a few dozen iterations. If |λ₁| < 1, the vector decays exponentially to zero (underflow). Normalizing x^{(k+1)} = Ax^{(k)} / ‖Ax^{(k)}‖ at every step keeps the vector on the unit sphere S^{n-1}, maintaining numerical stability indefinitely.
Algorithm precautions and stability rules
- Always normalize the iterate vector at every step (x ← y / ‖y‖₂) to prevent floating-point overflow or underflow.
- Verify that the initial vector x^{(0)} is not orthogonal to the dominant eigenvector v₁; use random entries.
- Check for eigenvalue clustering (|λ₂| ≈ |λ₁|); if the ratio exceeds 0.99, use Chebyshev acceleration or Inverse Iteration.
- Guard against complex conjugate dominant pairs; if iterates oscillate without converging, switch to Arnoldi iteration.
- Use the Rayleigh quotient to compute the eigenvalue estimate μ = xᵀAx rather than simple coordinate ratios.
- Maintain 64-bit precision to detect convergence down to machine tolerance.
Theoretical limitations
- Only computes the single dominant eigenvalue λ₁ and its corresponding eigenvector; cannot find interior eigenvalues directly.
- Convergence can be excruciatingly slow if the spectral gap |λ₁| − |λ₂| is narrow.
- Fails completely if the matrix has multiple dominant eigenvalues of equal magnitude but different phase.