Power Iteration

Power Iteration (also known as the power method or von Mises iteration) is an iterative eigenvalue algorithm in numerical linear algebra and computation and algorithms that computes the dominant eigenvalue λ₁ (the eigenvalue of largest modulus) and its corresponding dominant eigenvector v₁ of a square matrix A ∈ ℂ^n×n. Formulated by Austrian mathematician Richard von Mises and Hilda Pollaczek-Geiringer in 1929, the method constructs a sequence of normalized vectors inside the Krylov subspace span{b₀, Ab₀, A²b₀, ..., A^kb₀} through repeated matrix-vector multiplications and Rayleigh quotient evaluations.

Because power iteration requires only a matrix-vector product y = Ab_k at each step rather than explicit matrix storage or O(n³) similarity transformations, it operates in O(nnz(A)) time per iteration for sparse matrices with nnz(A) nonzero entries. While full-spectrum solvers such as the QR algorithm and singular value decomposition compute all n eigenvalues of dense matrices, power iteration and its shifted-inverse variants remain the primary computational engines for large-scale network centrality rankings—most notably Google's PageRank algorithm on web graphs with billions of vertices—as well as principal component initialization, Markov chain stationary distributions, and dominant Lyapunov or buckling mode extraction.

The Normalized Power Iteration Recurrence and Rayleigh Quotient

Let A ∈ ℂ^n×n be a diagonalizable matrix whose eigenvalues are ordered by decreasing magnitude with a strict spectral gap behind the dominant eigenvalue:

Starting from an initial unit vector b₀ ∈ ℂⁿ (‖b₀‖₂ = 1) that has a nonzero projection onto the dominant eigenspace, power iteration computes for k = 1, 2, 3, ...:

Periodic normalization by the Euclidean norm ‖w_k‖₂ prevents floating-point overflow (when |λ₁| > 1) or underflow (when |λ₁| < 1) under IEEE 754 arithmetic. At step k, the dominant eigenvalue λ₁ is estimated via the Rayleigh quotient r(b_k):

Spectral Expansion and Asymptotic Rate of Convergence

To analyze convergence, expand the starting vector b₀ in the basis of normalized eigenvectors {v₁, v₂, ..., *v_n*} satisfying Av_j = λ_jv_j:

Multiplying b₀ by A^k and factoring out c₁λ₁^k reveals how the subdominant eigencomponents attenuate geometrically with k:

Because |λ_j / λ₁| ≤ |λ₂ / λ₁| < 1 for all j ≥ 2, every subdominant term decays at the geometric damping ratio q = |λ₂ / λ₁|. Consequently, the angle θ(b_k, v₁) between the iterate b_k and the true dominant eigenvector v₁ satisfies:

For a Hermitian or real symmetric matrix (A = A^*), the eigenvectors {*v_j*} form an orthonormal basis, and the Rayleigh quotient r(b_k) is stationary at v₁ (∇r(v₁) = 0). As a result, while the eigenvector error ‖b_k − ±v₁‖₂ decays as O(|λ₂/λ₁|^k), the Rayleigh quotient eigenvalue error converges quadratically faster with respect to the spectral ratio:

Inverse Iteration and Cubic Rayleigh Quotient Iteration (RQI)

Direct power iteration can find only the largest-magnitude eigenvalue λ₁ and converges slowly when the spectral ratio |λ₂/λ₁| is close to 1 (for example, if |λ₂/λ₁| = 0.99, gaining one decimal digit of accuracy requires ln(10) / |ln(0.99)| ≈ 230 iterations). Both limitations are overcome by shifted inverse iteration (introduced by Helmut Wielandt in 1944): for any scalar shift σ ∈ ℂ not in the spectrum of A, the shifted-and-inverted matrix (A − σI)⁻¹ shares the identical eigenvectors v_j as A, but its eigenvalues are (λ_j − σ)⁻¹.

If σ is chosen close to a target interior eigenvalue λ_m so that |λ_m − σ| ≪ |λ_j − σ| for all j ≠ m, then (λ_m − σ)⁻¹ becomes the overwhelmingly dominant eigenvalue of (A − σI)⁻¹. Instead of computing the explicit matrix inverse, each step solves the linear system (A − σI)w_k = b_k−1 using a precomputed PA = LU factorization via Gaussian elimination or LU decomposition, converging at rate:

When the shift σ is updated dynamically at every step to equal the current Rayleigh quotient σ_k = r(b_k), the procedure becomes Rayleigh Quotient Iteration (RQI). For Hermitian matrices, RQI achieves local cubic convergence:

typically reaching full 16-digit IEEE double-precision accuracy in 3 to 5 iterations, at the cost of a fresh O(n³) (or O(n) tridiagonal) linear solve per step.

By the Perron–Frobenius theorem, every irreducible non-negative matrix A ≥ 0 possesses a positive real dominant eigenvalue r = ρ(A) (the spectral radius) with a strictly positive eigenvector v₁ > 0, and if A is primitive (aperiodic), power iteration converges to v₁ from any non-negative starting vector. In Google's PageRank algorithm (Brin & Page, 1998), the web hyperlink transition matrix P is modified with a teleportation damping factor α ∈ (0, 1) (typically α = 0.85) to form the primitive stochastic Google matrix G = αP + (1 − α)*ve*^T. Because the second eigenvalue of G satisfies |λ₂(G)| ≤ α = 0.85 while λ₁(G) = 1, power iteration reduces the residual by a factor of 0.85⁵⁰ ≈ 2.96 × 10⁻⁴ in just 50 sparse matrix-vector sweeps regardless of web graph size.

In data science and machine learning, applying power iteration to A^TA (via sequential products z = Ab_k, w = A^Tz without forming A^TA explicitly) extracts the leading right singular vector v₁ and largest singular value σ₁ = ‖A‖₂ = √λ₁(A^TA), providing fast spectral norm regularization and randomized singular value decomposition initialization.

What happens if the starting vector b₀ is orthogonal to the dominant eigenvector v₁ (c₁ = 0)?

In exact real arithmetic, if c₁ = 0, power iteration would converge to the second eigenvector v₂. In finite-precision IEEE 754 arithmetic, however, rounding errors in Ab_k introduce a nonzero perturbation of magnitude ~ε_mach ≈ 10⁻¹⁶ along v₁, which is then amplified by (λ₁/λ₂)^k until the iterate eventually rotates toward v₁. Choosing a random Gaussian starting vector b₀ ensures c₁ ≠ 0 with probability 1.

How can power iteration find the second dominant eigenvalue λ₂?

Once (λ₁, v₁) of a symmetric matrix A is computed, Hotelling's deflation forms the deflated matrix A₁ = A − λ₁v₁v₁^T, which replaces λ₁ with 0 while preserving all remaining eigenpairs (λ₂, v₂), ..., (λ_n, v_n). Alternatively, orthogonal subspace iteration iterates on an n × p block of vectors B_k and re-orthonormalizes its columns at each step via QR factorization—the direct precursor to the full QR algorithm.

Why does power iteration fail when the dominant eigenvalues form a complex conjugate pair?

For a real non-symmetric matrix A whose dominant eigenvalues are complex conjugates λ_1,2 = r e^±iθ (with |λ₁| = |λ₂| = r > |λ₃|), the ratio |λ₂/λ₁| = 1 violates the strict spectral gap condition. A single real vector b_k oscillates perpetually within the two-dimensional invariant subspace spanned by Re(v₁) and Im(v₁) without converging to a single direction, requiring a 2-column block power iteration or double-shift QR step to resolve the complex pair.

1. von Mises, R., & Pollaczek-Geiringer, H. (1929). "Praktische Verfahren der Gleichungsauflösung." Zeitschrift für Angewandte Mathematik und Mechanik (ZAMM), 9(1), 58–77. https://doi.org/10.1002/zamm.19290090105
2. Parlett, B. N. (1998). The Symmetric Eigenvalue Problem (Classics in Applied Mathematics, Vol. 20, Chapter 4: Power Iterations). SIAM. https://doi.org/10.1137/1.9781611971163
3. Page, L., Brin, S., Motwani, R., & Winograd, T. (1999). The PageRank Citation Ranking: Bringing Order to the Web. Stanford InfoLab Technical Report 1999-66. http://ilpubs.stanford.edu:8090/422/
4. Golub, G. H., & Van Loan, C. F. (2013). Matrix Computations (4th ed., Chapter 7.3: Power Iterations and Chapter 8.2: Symmetric Power Methods). Johns Hopkins University Press. https://doi.org/10.56021/9781421407944
5. Saad, Y. (2011). Numerical Methods for Large Eigenvalue Problems (Revised ed., Chapter 4: Subspace Iteration). SIAM. https://doi.org/10.1137/1.9781611970739