Eigenvalue Algorithms

Numerical eigenvalue methods range from power iteration for a single dominant value to the QR algorithm for the full spectrum of a matrix.

Why iteration is required

Eigenvalues are the roots of the characteristic polynomial, but computing that polynomial and finding its roots is numerically disastrous for all but tiny matrices. Instead, eigenvalue algorithms are iterative, applying orthogonal transformations that reveal the eigenvalues while keeping the problem well conditioned.

Power iteration

Power iteration repeatedly multiplies a vector by A and normalizes it; the vector aligns with the eigenvector of the largest-magnitude eigenvalue, and the Rayleigh quotient gives that eigenvalue. It is simple and needs only matrix-vector products, but it finds only the dominant eigenpair and converges at a rate set by the ratio of the two largest eigenvalues.

import numpy as np
def power_iteration(A, iters=1000):
    x = np.random.rand(A.shape[0])
    for _ in range(iters):
        x = A @ x; x = x/np.linalg.norm(x)
    return (x @ A @ x)/(x @ x), x   # eigenvalue, eigenvector

Shifts and inverse iteration

Inverse iteration applies power iteration to (A - sigma I)^{-1}, converging to the eigenvalue nearest the shift sigma. With a good shift it converges very fast and finds interior eigenvalues, at the cost of a linear solve each step. Rayleigh-quotient iteration updates the shift adaptively for cubic convergence.

The full spectrum

To find all eigenvalues, the QR algorithm is the standard: it repeatedly factors the matrix into Q times R and reforms R times Q, driving the matrix toward triangular form whose diagonal holds the eigenvalues. For large sparse matrices, Krylov methods (Lanczos for symmetric, Arnoldi for general) find a subset of the spectrum efficiently.

Eigenvalue computation underlies stability and mode analysis in physics, including magnetohydrodynamic stability spectra examined in breeder Hyperion design studies.