Mathematics.

numerical linear algebra

Numerical Eigenvalue Algorithms

Numerical Analysis35 minDifficulty7 out of 10

You should know: eigenvalues and eigenvectors, gaussian elimination

Overview

Finding eigenvalues exactly means finding roots of the characteristic polynomial det(A − λI) = 0, but for n larger than 4 there is no general algebraic formula for polynomial roots (a consequence of the Abel-Ruffini theorem), so every practical eigenvalue algorithm is necessarily iterative and approximate. Power iteration is the simplest such method: repeatedly multiply a vector by A and rescale, and the vector's direction converges to the eigenvector of the largest-magnitude (dominant) eigenvalue while the scaling factor converges to that eigenvalue itself. Production numerical libraries (LAPACK, and hence NumPy, MATLAB, etc.) instead use the QR algorithm, which repeatedly factors a matrix as A = QR and reassembles it as A' = RQ; this sequence of similar matrices converges to a triangular (or block-triangular) form whose diagonal entries are all the eigenvalues at once, typically after first reducing A to a cheaper-to-iterate Hessenberg form. Both families underlie everything from Google's original PageRank (power iteration on the web graph) to vibration-mode analysis in structural engineering (QR-based solvers on stiffness/mass matrices).

Intuition

Think of power iteration as repeatedly stretching a rubber vector using the matrix A as a stretching machine. Any starting vector is some mix of all the eigenvector directions, but each pass through A stretches the dominant eigenvector's component the most (that's what 'dominant eigenvalue' means), so after repeated applications and rescaling, the other directions shrink into irrelevance and only the dominant direction survives — the vector 'snaps' into alignment with the top eigenvector, and the amount it grew each step reveals the eigenvalue. The QR algorithm is subtler: it doesn't chase a single eigenvector but instead repeatedly reshuffles the whole matrix through a rotate-and-recombine step that provably drives all the off-diagonal 'coupling' between eigen-directions toward zero simultaneously, revealing every eigenvalue on the diagonal at once — like slowly untangling several ropes at the same time rather than pulling out one at a time.

Formal Definition

Definition

Power iteration starts from an initial vector v⁽⁰⁾ and repeatedly applies A, normalizing after each step; the Rayleigh quotient extracts the current eigenvalue estimate:

w(k)=Av(k1),v(k)=w(k)w(k)w^{(k)} = A v^{(k-1)}, \qquad v^{(k)} = \dfrac{w^{(k)}}{\|w^{(k)}\|_\infty}
Power iteration update and normalization
λ(k)=(v(k))TAv(k)(v(k))Tv(k)\lambda^{(k)} = \dfrac{(v^{(k)})^{T} A\, v^{(k)}}{(v^{(k)})^{T} v^{(k)}}
Rayleigh quotient — the eigenvalue estimate at step k
A(0)=A;A(k1)=Q(k)R(k),A(k)=R(k)Q(k)    triangularA^{(0)} = A; \quad A^{(k-1)} = Q^{(k)} R^{(k)}, \quad A^{(k)} = R^{(k)} Q^{(k)} \;\to\; \text{triangular}
QR algorithm — similar matrices converging to (block) triangular form

Worked Examples

  1. Multiply A by v⁽⁰⁾.

    w(1)=Av(0)=(21+1011+20)=(21)w^{(1)} = A v^{(0)} = \begin{pmatrix}2\cdot 1 + 1\cdot 0\\ 1\cdot 1 + 2\cdot 0\end{pmatrix} = \begin{pmatrix}2\\1\end{pmatrix}
  2. Normalize by the largest entry (2) to get v⁽¹⁾, and read off the Rayleigh quotient estimate.

    v(1)=(1, 0.5),λ(1)=v(1)TAv(1)v(1)Tv(1)=2.81.25=2.8v^{(1)} = (1,\ 0.5), \qquad \lambda^{(1)} = \dfrac{v^{(1)T}Av^{(1)}}{v^{(1)T}v^{(1)}} = \dfrac{2.8}{1.25} = 2.8
  3. Second application: multiply A by v⁽¹⁾.

    w(2)=Av(1)=(2(1)+1(0.5)1(1)+2(0.5))=(2.52)w^{(2)} = A v^{(1)} = \begin{pmatrix}2(1)+1(0.5)\\ 1(1)+2(0.5)\end{pmatrix} = \begin{pmatrix}2.5\\2\end{pmatrix}
  4. Normalize by 2.5 to get v⁽²⁾, and recompute the Rayleigh quotient.

    v(2)=(1, 0.8),λ(2)=v(2)TAv(2)v(2)Tv(2)=7.242.642.9758v^{(2)} = (1,\ 0.8), \qquad \lambda^{(2)} = \dfrac{v^{(2)T}Av^{(2)}}{v^{(2)T}v^{(2)}} = \dfrac{7.24}{2.64} \approx 2.9758

Answer: λ⁽¹⁾ = 2.8, λ⁽²⁾ ≈ 2.976, rapidly converging to the true dominant eigenvalue λ = 3 (A's eigenvalues are 3 and 1, with dominant eigenvector direction (1,1)).

Practice Problems

Difficulty 5/10

For A = [[2,1],[1,2]] and v⁽⁰⁾=(1,0), compute w⁽¹⁾ = Av⁽⁰⁾ (before normalization).

Difficulty 6/10

Using v⁽¹⁾ = (1, 0.5) from the worked example, compute the Rayleigh quotient λ⁽¹⁾ = (v⁽¹⁾ᵀAv⁽¹⁾)/(v⁽¹⁾ᵀv⁽¹⁾) for A = [[2,1],[1,2]].

Difficulty 7/10

Google's original PageRank algorithm models web importance as the dominant eigenvector of a huge (billions of rows) sparse link matrix. Why is power iteration, rather than computing the full characteristic polynomial, the only feasible approach here?

Quiz

Power iteration converges to the eigenvector associated with:
Why must eigenvalue algorithms be iterative rather than using a closed-form root formula for n ≥ 5?
The QR algorithm produces all eigenvalues of A by:

Summary

  • Since no algebraic root formula exists for polynomials of degree ≥ 5 (Abel-Ruffini), all practical eigenvalue algorithms for n ≥ 5 are iterative.
  • Power iteration repeatedly applies A and rescales, converging to the dominant eigenvector; the Rayleigh quotient gives the eigenvalue estimate at each step, e.g. 2.8 then ≈2.976 toward the true value 3 for A=[[2,1],[1,2]].
  • The QR algorithm (used in LAPACK/NumPy/MATLAB) repeatedly factors A=QR and reassembles RQ, producing a sequence of similar matrices that converges to triangular form revealing all eigenvalues at once.
  • These methods underlie PageRank (power iteration on the web graph) and vibration/stability analysis in engineering (QR-based solvers on stiffness/mass matrices).

References

  1. BookGolub, G. H. & Van Loan, C. F. Matrix Computations, 4th ed. Johns Hopkins University Press, 2013.