Mathematical foundations of Inverse Iteration & Shift.
Rayleigh Quotient Iteration (RQI)—conceived in its modern form by American numerical analyst Richard H. Ostrowski in 1958 and rigorously analyzed by James H. Wilkinson in 1965—is widely regarded as one of the most astonishingly fast algorithms in numerical linear algebra. While standard Power Iteration is restricted to finding the single dominant exterior eigenvalue with slow linear convergence, applied scientists frequently need to locate specific interior resonant frequencies—such as the acoustic modes of a concert hall, the flutter frequencies of an aircraft wing, or the quantum energy levels of a molecule. Inverse iteration isolates any targeted interior eigenvalue by applying power iteration to the inverted shifted matrix (A − σI)⁻¹. The eigenvalue closest to the shift σ is mapped to the largest eigenvalue of the inverted operator, isolating the desired mode. The stroke of mathematical genius in Rayleigh Quotient Iteration is to make the shift dynamic: instead of keeping σ constant, the algorithm updates the shift at every single step using the current vector's own Rayleigh quotient: σ_k = (x_kᵀAx_k) / (x_kᵀx_k). When applied to symmetric matrices, this creates an extraordinary synergy: the Rayleigh quotient provides quadratic accuracy in the eigenvalue (|σ_k − λ*| = O(‖e_k‖²)), while the shifted inverse solve provides quadratic accuracy in the vector, resulting in an exact cubic convergence order α = 3.000. In cubic convergence, the number of correct significant decimal digits triples at every iteration: an initial estimate with 2 correct digits advances to 6 digits in one step, 18 digits in two steps, and 54 digits in three steps—hitting 64-bit machine precision in a mere trio of iterations. Furthermore, in one of the most celebrated paradoxes in numerical analysis, J. H. Wilkinson proved that although the matrix A − σ_k I becomes nearly singular as σ_k approaches λ* (with condition number exploding toward infinity), the resulting numerical error aligns almost perfectly with the desired eigenvector, meaning ill-conditioning actually accelerates convergence! In this computational laboratory, students evaluate Rayleigh Quotient Iteration applied to an interior eigenvalue λ* = 3.000000 of a symmetric matrix, track dynamic shift updates, and extract the theoretical cubic convergence order α = 3.000 via logarithmic linear regression.
Viva voce questions & analytical theory
What is the mathematical mechanism behind Inverse Power Iteration with a constant shift?
Let A have eigenvalues λ_i and eigenvectors v_i. The shifted matrix A − σI has eigenvalues λ_i − σ with the exact same eigenvectors. Its inverse (A − σI)⁻¹ has eigenvalues μ_i = 1 / (λ_i − σ). When power iteration is applied to (A − σI)⁻¹, the dominant eigenvalue corresponds to whichever λ_i is closest to the shift σ: |λ_target − σ| = min_i |λ_i − σ|. Because 1/|λ_target − σ| >> 1/|λ_other − σ|, the algorithm converges to the target interior eigenvector v_target with linear rate |(λ_target − σ) / (λ_closest − σ)|.
Why does Rayleigh Quotient Iteration (RQI) achieve cubic convergence (α = 3.0) on symmetric matrices?
In Rayleigh Quotient Iteration, the shift σ_k is not held constant: it is updated at every iteration using the current vector's Rayleigh quotient: σ_k = (x_kᵀAx_k) / (x_kᵀx_k). For symmetric matrices, the Rayleigh quotient has quadratic accuracy: |σ_k − λ*| = O(‖x_k − v*‖²). When this quadratically accurate shift is substituted into the inverse iteration step, the error in the next vector scales as ‖x_{k+1} − v*‖ = O(|σ_k − λ*| · ‖x_k − v*‖) = O(‖x_k − v*‖³). D. K. Faddeev and J. H. Wilkinson proved that this mutual feedback produces an extraordinary cubic convergence order α = 3.000, tripling the number of correct digits per iteration.
Why doesn't the near-singularity of (A − σ_k I) corrupt the computed eigenvector as σ_k → λ*?
As σ_k approaches the exact eigenvalue λ*, the matrix A − σ_k I becomes nearly singular, with condition number κ → ∞. In standard linear systems, an ill-conditioned matrix produces disastrous errors. However, Wilkinson proved a celebrated result: the round-off error incurred when solving (A − σ_k I)w = x_k lies almost entirely in the direction of the desired eigenvector v*! What would be fatal numerical noise in a general linear system is, in inverse iteration, precisely the component that aligns w with the true eigenvector. Consequently, ill-conditioning actually accelerates convergence rather than ruining it.
How does RQI compare to the QR algorithm with Wilkinson shifts?
The QR algorithm with Wilkinson shifts computes all n eigenvalues simultaneously with cubic convergence in O(n³) operations per sweep. In contrast, RQI computes a single specific eigenpair (λ*, v*) from a targeted initial guess in only 2 to 4 iterations. RQI is therefore vastly faster when only a few interior vibrational modes or buckling modes are needed in engineering dynamics.
What is the computational cost per iteration of RQI compared to fixed-shift inverse iteration?
In fixed-shift inverse iteration, the matrix A − σI is factored into LU or Cholesky only once at the beginning (costing (2/3)n³ FLOPs), and subsequent steps require only O(n²) back-substitutions. In RQI, because the shift σ_k changes at every iteration, the matrix A − σ_k I must be refactored at every single step, costing (2/3)n³ FLOPs per iteration. However, because RQI converges in only 3 steps (cubic convergence) compared to dozens of steps for fixed shift, RQI is often cheaper in total operations.
What happens if RQI is initialized with a vector equidistant between two eigenvectors?
If x^{(0)} is perfectly symmetric between two eigenvectors (e.g. x^{(0)} = (v₁ + v₂)/√2), the Rayleigh quotient lands midway between the two eigenvalues, and the iteration can theoretically stall or oscillate. However, in finite-precision arithmetic, round-off noise immediately breaks the symmetry, causing RQI to plunge rapidly toward one of the two eigenvalues with cubic speed.
Can RQI converge to a different eigenvalue than the one closest to the initial guess?
Yes. While RQI converges locally with cubic speed, its global basin of attraction has a fractal boundary in the complex plane. A poorly chosen initial guess x^{(0)} can occasionally jump to an unexpected eigenvalue far from the initial shift.
Algorithm precautions and stability rules
- Factor A − σ_k I using partial pivoting (or LDLᵀ with Bunch-Kaufman pivoting) to prevent catastrophic division by zero.
- Do not abort computation when the condition number κ(A − σ_k I) becomes huge; this signals that the shift is extremely close to the true eigenvalue.
- Normalize the eigenvector iterate w^{(k+1)} at every step: x^{(k+1)} = w^{(k+1)} / ‖w^{(k+1)}‖₂.
- Verify that matrix A is symmetric; unsymmetric matrices yield quadratic convergence (α = 2) rather than cubic (α = 3).
- Select an initial shift σ₀ within the basin of attraction of the desired interior eigenvalue.
- Maintain 64-bit IEEE 754 precision throughout; 3 iterations are typically sufficient to hit machine epsilon.
Theoretical limitations
- Requires refactoring the matrix A − σ_k I at every iteration, costing O(n³) operations per step.
- Basins of attraction can have complex fractal boundaries, occasionally converging to an unintended eigenvalue if initialized poorly.
- For non-symmetric matrices, convergence drops from cubic (α = 3) to quadratic (α = 2).