Chapter 19: Power Method and Inverse Iteration

1. Introduction and Motivation

There are two broad approaches to numerically solving the eigenvalue problem $A\mathbf{v} = \lambda \mathbf{v}$ for an $N \times N$ matrix $A$.

1. Methods that compute all eigenvalues simultaneously: the QR algorithm (Chapter 18)
2. Iterative methods that selectively compute specific eigenvalues: the power method and its variants

This chapter covers approach (2) in detail. The power method is the most fundamental technique for finding the eigenvalue of largest absolute value, and inverse iteration is a variant that can efficiently find the eigenvalue closest to any target value.

2. Power Method

2.1 Prerequisites

2.2 Algorithm

3. Convergence Proof

3.1 Basic Property: $A^k \mathbf{v}_n = \lambda_n^k \mathbf{v}_n$

From the eigenvalue equation $A\mathbf{v}_n = \lambda_n \mathbf{v}_n$, the following holds by induction.

Proof: By induction on $k$. The case $k = 1$ holds by definition. Assuming the result holds for $k-1$,

3.2 Convergence via Eigenvector Expansion

Proof: Expand the initial vector in the eigenvector basis.

Applying $A^k$ and using equation $\eqref{eq:Akvn}$,

Since $|\lambda_1| > |\lambda_i|$ for $i \geq 2$, we have $|\lambda_i / \lambda_1| < 1$, so

Therefore, as $k \to \infty$ the bracketed expression converges to $\mathbf{v}_1$, and after normalization,

$\square$

3.3 Eigenvalue Extraction

With the normalization $\|\mathbf{x}^{(k)}\| = 1$, the eigenvalue can be approximated via the Rayleigh quotient.

Since $\mathbf{x}^{(k)} \to \mathbf{v}_1 / \|\mathbf{v}_1\|$, we have $A\mathbf{x}^{(k)} \to \lambda_1 \mathbf{x}^{(k)}$, so

3.4 Necessity of Normalization

If $|\lambda_1| \neq 1$, iterating without normalization causes overflow ($|\lambda_1| > 1$) or underflow ($|\lambda_1| < 1$) due to the factor $\lambda_1^k$. Normalization at each step prevents this.

4. Convergence Rate Analysis

The convergence ratio $r = |\lambda_2/\lambda_1|$ takes values between 0 and 1; the smaller $r$, the faster the convergence. The number of iterations required to reduce the relative eigenvalue error below $\varepsilon$ is

5. Inverse Power Method (Inverse Iteration)

5.1 Principle

If $A$ has eigenvalues $\lambda_1, \lambda_2, \ldots, \lambda_N$, then $A^{-1}$ has eigenvalues $1/\lambda_1, 1/\lambda_2, \ldots, 1/\lambda_N$ with the same eigenvectors. Applying the power method to $A^{-1}$ therefore converges to the eigenvector corresponding to the largest $|1/\lambda_i|$, i.e., the smallest $|\lambda_i|$.

5.2 Algorithm

6. Shifted Inverse Power Method

6.1 Introducing the Shift

Introduce a target shift value $\sigma$ and apply the power method to $(A - \sigma I)^{-1}$. Since the eigenvalues of $(A - \sigma I)$ are $\lambda_i - \sigma$, the eigenvalues of $(A - \sigma I)^{-1}$ are

For the eigenvalue $\lambda_j$ closest to $\sigma$, $|1/(\lambda_j - \sigma)|$ is the largest, so the power method converges to the eigenvector of $\lambda_j$.

6.2 Algorithm

6.3 Convergence Rate

7. Rayleigh Quotient Iteration

7.1 Properties of the Rayleigh Quotient

7.2 Algorithm

Update the shift $\sigma$ in shifted inverse iteration at every step using the Rayleigh quotient.

8. Numerical Examples

8.1 Power Method

Python implementation

import numpy as np

def power_method(A, max_iter=200, tol=1e-12):
    """Return the dominant eigenvalue and eigenvector of A.

    A: (N, N) ndarray, diagonalizable with |λ_1| > |λ_2|.
    Returns: (eigenvalue, eigenvector, iters)
    """
    n = A.shape[0]
    x = np.random.default_rng(0).standard_normal(n)
    x /= np.linalg.norm(x)
    mu_prev = np.inf
    for k in range(1, max_iter + 1):
        y = A @ x
        x = y / np.linalg.norm(y)
        mu = x @ A @ x                      # Rayleigh quotient
        if abs(mu - mu_prev) < tol * max(1.0, abs(mu)):
            return mu, x, k
        mu_prev = mu
    raise RuntimeError(f"did not converge in {max_iter} iterations")

A = np.array([[2.0, 1.0], [1.0, 3.0]])
mu, v, k = power_method(A)
print(mu, v, k)
# 3.6180339887498945 [0.5257311  0.85065081] ~60 iters

To extend this to shifted inverse iteration, compute lu = scipy.linalg.lu_factor(A - sigma*np.eye(n)) once before the loop and replace the matrix-vector product with y = scipy.linalg.lu_solve(lu, x).

8.2 Shifted Inverse Iteration

8.3 Rayleigh Quotient Iteration

9. Convergence Visualization

9.1 Power Method Convergence Graph

9.2 Comparison of Three Methods

10. Implementation Notes

10.1 Choice of Initial Vector

Using a random vector as the initial value virtually guarantees $c_1 \neq 0$. In practice, a random initial vector is sufficient, but if prior information is available, it can be exploited to reduce the number of iterations.

10.2 Convergence Criteria

Either of the following can be used:

• Eigenvalue change: $|\mu^{(k)} - \mu^{(k-1)}| < \varepsilon$
• Residual: $\|A\mathbf{x}^{(k)} - \mu^{(k)}\mathbf{x}^{(k)}\| < \varepsilon$ (residual norm)

Residual-based criteria are more reliable. When the eigenvalue is very small, absolute error is insufficient, so implementations often use the relative error

10.3 When the Power Method Fails

10.4 Computational Cost Comparison

10.5 Relationship to Other Methods

• QR Algorithm: Computes all eigenvalues simultaneously. Can be viewed as an implicit parallel execution of the power method.
• Lanczos Method: Eigenvalue computation for large sparse matrices. An extension of the power method to Krylov subspaces.
• LAPACK Strategy: Compute eigenvalues with the QR algorithm, then compute eigenvectors using shifted inverse iteration.

! Fortran: all eigenvalues + eigenvectors
CALL DSYEV('V', 'U', N, A, LDA, W, WORK, LWORK, INFO)
! 'V' = also compute eigenvectors, 'U' = reference upper triangle

11. Summary

• The power method is the most fundamental technique for finding the eigenvalue of largest absolute value. Its convergence rate is $O(|\lambda_2/\lambda_1|^k)$, and convergence is slow when $|\lambda_2/\lambda_1|$ is close to 1.
• Inverse iteration applies the power method to $A^{-1}$ to find the eigenvalue of smallest absolute value.
• Shifted inverse iteration applies the power method to $(A - \sigma I)^{-1}$ to find the eigenvalue closest to $\sigma$. This is the standard LAPACK strategy for computing eigenvectors.
• Rayleigh quotient iteration dynamically updates the shift using the Rayleigh quotient, achieving cubic convergence for symmetric matrices.

References

• G. H. Golub, C. F. Van Loan, Matrix Computations, 4th ed., Johns Hopkins University Press, 2013.
• L. N. Trefethen, D. Bau III, Numerical Linear Algebra, SIAM, 1997.
• B. N. Parlett, The Symmetric Eigenvalue Problem, SIAM, 1998.
• Y. Saad, Numerical Methods for Large Eigenvalue Problems, Revised ed., SIAM, 2011.
• Wikipedia: Power iteration