Power Iteration: Compute Dominant Eigenvalues & Eigenvectors

In many real-world applications—from Google's PageRank algorithm to Principal Component Analysis (PCA) in machine learning—we need to find the largest eigenvalue and its corresponding eigenvector of a large matrix. Computing all eigenvalues via full eigendecomposition is often too expensive for huge matrices (e.g., millions of dimensions). The power iteration method offers a simple, iterative approach to compute the dominant eigenvalue and eigenvector. It is the foundation of more advanced algorithms like the QR algorithm and Lanczos method. In this tutorial, you'll learn how power iteration works, implement it in Python, understand its convergence properties, and see how to use it in production systems. We'll also cover common pitfalls, debugging strategies, and interview questions to solidify your understanding.

What is Power Iteration?

Power iteration is an iterative algorithm to find the eigenvalue with the largest absolute value (dominant eigenvalue) and its corresponding eigenvector of a square matrix. Given a matrix A (n×n) with eigenvalues λ₁, λ₂, ..., λₙ such that |λ₁| > |λ₂| ≥ ... ≥ |λₙ|, the algorithm starts with a random vector b₀ and repeatedly applies the update:

b_{k+1} = A b_k / ||A b_k||

Under mild conditions, b_k converges to the eigenvector associated with λ₁, and the Rayleigh quotient λ_k = (b_k^T A b_k) / (b_k^T b_k) converges to λ₁.

The algorithm is simple, memory-efficient (only matrix-vector products), and works well for large sparse matrices. However, it only finds the dominant eigenvalue and requires a good eigenvalue gap.

import numpy as np

def power_iteration(A, num_simulations=100, tolerance=1e-6):
    n = A.shape[0]
    b = np.random.rand(n)
    b = b / np.linalg.norm(b)
    
    for i in range(num_simulations):
        b_new = A @ b
        b_new_norm = np.linalg.norm(b_new)
        b_new = b_new / b_new_norm
        
        # Rayleigh quotient
        eigenvalue = b_new @ (A @ b_new)
        
        # Check convergence
        if np.linalg.norm(b_new - b) < tolerance:
            break
        b = b_new
    
    return eigenvalue, b

# Example usage
A = np.array([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])
eigval, eigvec = power_iteration(A)
print("Dominant eigenvalue:", eigval)
print("Corresponding eigenvector:", eigvec)

Convergence Analysis

The convergence rate of power iteration is linear, with the error decreasing by a factor of |λ₂/λ₁| each iteration. If the eigenvalues are close, convergence is slow. The algorithm converges if: - The matrix has a strictly dominant eigenvalue (|λ₁| > |λ₂|). - The initial vector has a component in the direction of the dominant eigenvector (almost always true for random vectors).

The error after k iterations is O(|λ₂/λ₁|^k). For example, if |λ₂/λ₁| = 0.9, after 100 iterations the error is about 0.9^100 ≈ 0.000026, which is good. But if |λ₂/λ₁| = 0.999, you need ~2300 iterations for similar accuracy.

import numpy as np
import matplotlib.pyplot as plt

def power_iteration_track(A, num_simulations=100):
    n = A.shape[0]
    b = np.random.rand(n)
    b = b / np.linalg.norm(b)
    errors = []
    true_eigvec = np.linalg.eig(A)[1][:,0]  # dominant eigenvector
    for i in range(num_simulations):
        b_new = A @ b
        b_new = b_new / np.linalg.norm(b_new)
        # error: angle between current and true eigenvector
        cos_angle = np.abs(b_new @ true_eigvec)
        errors.append(np.arccos(cos_angle))
        b = b_new
    return errors

A = np.diag([1, 0.9, 0.5])  # eigenvalues: 1, 0.9, 0.5
errors = power_iteration_track(A, 50)
print("First few errors:", errors[:5])

Implementing Power Iteration in Python

Here's a robust implementation that handles sparse matrices and includes convergence checks. We'll use SciPy's sparse matrix format for efficiency.

Key steps: 1. Initialize a random vector. 2. Normalize it. 3. Multiply by the matrix (using @ for dense, or dot for sparse). 4. Normalize the result. 5. Compute Rayleigh quotient to estimate eigenvalue. 6. Check convergence: if the change in eigenvector norm is below tolerance, stop. 7. Return eigenvalue and eigenvector.

We also include a maximum iteration limit to avoid infinite loops.

import numpy as np
from scipy.sparse import issparse, csr_matrix

def power_iteration_robust(A, max_iter=1000, tol=1e-8):
    if issparse(A):
        n = A.shape[0]
    else:
        n = A.shape[0]
    b = np.random.randn(n)
    b = b / np.linalg.norm(b)
    
    for i in range(max_iter):
        if issparse(A):
            b_new = A.dot(b)
        else:
            b_new = A @ b
        b_new_norm = np.linalg.norm(b_new)
        if b_new_norm == 0:
            raise ValueError("Matrix is singular or zero vector produced.")
        b_new = b_new / b_new_norm
        
        # Rayleigh quotient
        if issparse(A):
            eigenvalue = b_new.dot(A.dot(b_new))
        else:
            eigenvalue = b_new @ (A @ b_new)
        
        # Check convergence
        if np.linalg.norm(b_new - b) < tol:
            break
        b = b_new
    
    return eigenvalue, b

# Example with sparse matrix
A_sparse = csr_matrix([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])
eigval, eigvec = power_iteration_robust(A_sparse)
print("Eigenvalue:", eigval)
print("Eigenvector:", eigvec)

Applications of Power Iteration

Power iteration is used in many fields:

1. PageRank: Google's algorithm uses power iteration on the Google matrix to compute the principal eigenvector (PageRank scores). The matrix is huge but sparse, making power iteration ideal.
2. Principal Component Analysis (PCA): For high-dimensional data, power iteration can find the first principal component (dominant eigenvector of the covariance matrix). This is often done via the power method on the data matrix directly (without forming the covariance matrix).
3. Network Analysis: Eigenvector centrality measures the influence of nodes in a network. It is the dominant eigenvector of the adjacency matrix.
4. Markov Chains: The stationary distribution of a Markov chain is the dominant eigenvector of the transition matrix.
5. Spectral Clustering: Uses eigenvectors of the Laplacian matrix; power iteration can find the second eigenvector (Fiedler vector) by using a shift.

In each case, the matrix is often large and sparse, and only the dominant eigenvector is needed.

import numpy as np

def pagerank_power_iteration(adj_matrix, damping=0.85, max_iter=100, tol=1e-6):
    n = adj_matrix.shape[0]
    # Build Google matrix: G = d*M + (1-d)/n * ones
    # M is column-stochastic: each column sums to 1
    out_degree = adj_matrix.sum(axis=0)
    M = adj_matrix / out_degree  # handle zero columns separately
    # Power iteration on G
    r = np.ones(n) / n
    for i in range(max_iter):
        r_new = damping * M @ r + (1-damping) * np.ones(n) / n
        if np.linalg.norm(r_new - r) < tol:
            break
        r = r_new
    return r

# Example: small web graph
adj = np.array([[0,1,1,0], [1,0,0,0], [1,0,0,1], [0,0,1,0]])
pr = pagerank_power_iteration(adj)
print("PageRank scores:", pr)

Advanced Variants: Shifted and Inverse Iteration

To find eigenvalues other than the dominant one, we can use shifted power iteration or inverse iteration.

Shifted Power Iteration: Apply power iteration to (A - μI)⁻¹. This converges to the eigenvalue closest to μ. The shift μ can be chosen based on prior knowledge.

Inverse Iteration: Same as shifted with μ=0, finds the smallest eigenvalue.

Rayleigh Quotient Iteration: Uses the current eigenvalue estimate as shift, updating each iteration. This achieves cubic convergence but requires solving a linear system each step.

These variants are more powerful but also more computationally expensive per iteration.

import numpy as np
from scipy.linalg import solve

def inverse_iteration(A, mu=0, max_iter=100, tol=1e-8):
    n = A.shape[0]
    b = np.random.randn(n)
    b = b / np.linalg.norm(b)
    I = np.eye(n)
    for i in range(max_iter):
        # Solve (A - mu I) x = b
        try:
            x = solve(A - mu * I, b)
        except np.linalg.LinAlgError:
            raise ValueError("Matrix is singular; shift may be an eigenvalue.")
        b_new = x / np.linalg.norm(x)
        # Rayleigh quotient
        eigenvalue = b_new @ (A @ b_new)
        if np.linalg.norm(b_new - b) < tol:
            break
        b = b_new
    return eigenvalue, b

A = np.array([[2, -1, 0], [-1, 2, -1], [0, -1, 2]])
# Find eigenvalue near 0.5
eigval, eigvec = inverse_iteration(A, mu=0.5)
print("Eigenvalue near 0.5:", eigval)

Common Pitfalls and How to Avoid Them

1. Initial vector orthogonal to dominant eigenvector: Very unlikely with random vectors, but if the matrix has symmetries, it can happen. Use a random vector with a small random perturbation.
2. Slow convergence due to close eigenvalues: Monitor the Rayleigh quotient; if it oscillates, consider using a shift or Aitken acceleration.
3. Numerical instability: Normalize each iteration to avoid overflow. Use double precision.
4. Non-diagonalizable matrices: Power iteration still works for defective matrices, but convergence may be slower.
5. Complex eigenvalues: Power iteration works for real matrices with complex eigenvalues, but the vector may not converge (oscillate). Use the real part or consider using the QR algorithm.
6. Zero eigenvalue: If the dominant eigenvalue is zero, the matrix is singular and power iteration will converge to zero vector. Check for singularity beforehand.

import numpy as np

# Example: matrix with equal magnitude eigenvalues
A = np.array([[0, 1], [1, 0]])  # eigenvalues: 1 and -1
# Power iteration may oscillate between eigenvectors
b = np.array([1.0, 0.0])
for i in range(10):
    b = A @ b
    b = b / np.linalg.norm(b)
    print(f"Iter {i}: {b}")