Why Eigenvalue Methods Matters
A matrix describes how a system transforms a vector. Most vectors change direction, but some special directions are only stretched or reversed. These are eigenvectors, and their stretch factors are eigenvalues. They reveal dominant directions in data, vibration modes in mechanical models, and growth or decay in repeated linear transformations.
Different questions call for different methods. Power iteration seeks a dominant direction by repeated multiplication. Inverse iteration uses shifted linear solves to favor an eigenvalue near a chosen target. The QR algorithm changes coordinates repeatedly to expose a whole spectrum. Understanding their convergence conditions is as important as implementing their loops.
This lesson explains the broader ideas, then scopes its complete programs to real symmetric matrices. Such matrices have real eigenvalues and an orthonormal eigenbasis, making the geometry and residual checks especially clear. The implementations include iteration limits, pivoted solves for inverse iteration, and an unshifted Householder QR algorithm. They explicitly report nonconvergence rather than assuming every input will work.
Eigenvalue methods search for directions satisfying $Av=\lambda v$. Multiplication emphasizes large-magnitude eigenvalues, a shifted inverse emphasizes nearby eigenvalues, and orthogonal similarity transformations preserve the spectrum while making its structure easier to read.
Core Concepts and Terminology
An eigenpair consists of a scalar $\lambda$ and a nonzero vector $v$ with $Av=\lambda v$. Multiplying $v$ by any nonzero scalar gives the same eigendirection, so algorithms normalize their iterates. Negative eigenvalues can reverse a vector's sign on every multiplication without preventing convergence of its direction.
For a nonzero real vector $x$, the Rayleigh quotient is $\theta(x)=x^TAx/(x^Tx)$. It estimates an eigenvalue, but the estimate alone is insufficient. Always inspect the eigenpair residual $r=Ax-\theta(x)x$. Our nonzero-matrix stopping measure is
$$ \rho(x)=\frac{\lVert Ax-\theta(x)x\rVert_2} {\lVert A\rVert_F\lVert x\rVert_2}, \qquad \lVert A\rVert_F=\sqrt{\sum_{i,j}A_{ij}^2}. $$
The code evaluates the equivalent expression after uniformly scaling $A$, using hypot for norms. The zero matrix is treated separately through a reference value of $1$; every nonzero vector is then an exact eigenvector with eigenvalue zero.
A dominant eigenvalue means largest magnitude, not necessarily largest algebraic value. If $\lvert\lambda_1\rvert\gt\lvert\lambda_2\rvert$ and the start has a nonzero component in the first eigendirection, power iteration tends toward that direction. If this component is absent, another eigenpair can have zero residual while the intended dominant target is missed.
For symmetric $A$, write $A=V\Lambda V^T$, with orthogonal $V$. A general real nonsymmetric matrix can instead have complex eigenvalues, nonorthogonal eigenvectors, or no full eigenbasis. The code rejects nonsymmetric inputs; it is not a general real-Schur solver.
| Term | Meaning | Role or example |
|---|---|---|
| Eigenpair $(\lambda,v)$ | $Av=\lambda v$ with $v\ne0$ | For rows $(2,1),(1,2)$, direction $(1,1)^T$ has eigenvalue $3$ |
| Spectral radius | $\max_i\lvert\lambda_i\rvert$ | Equals $3$ in the sample |
| Rayleigh quotient $\theta(x)$ | $x^TAx/(x^Tx)$ | First power update gives $14/5$ |
| Residual $r$ | $Ax-\theta x$ | Tests an approximate eigenpair against original $A$ |
| Relative residual $\rho$ | Residual normalized by $\lVert A\rVert_F\lVert x\rVert_2$ | The program stops when $\rho\le10^{-10}$ |
| Magnitude gap | Separation between the largest eigenvalue magnitudes | Power's sample direction error decreases at a factor near $1/3$ per step |
| Shift $\mu$ | Target used in $A-\mu I$ | $\mu=3/2$ favors eigenvalue $1$ over $3$ |
| Inverse iteration | Repeated solves with $A-\mu I$, followed by normalization | Factor the shifted matrix once for a fixed shift |
| Orthogonal similarity | $T_{k+1}=Q_k^TT_kQ_k$ | Preserves all eigenvalues in exact arithmetic |
| Accumulated basis $Z_k$ | Product of QR orthogonal factors | Its columns approach eigenvectors when the symmetric iterate becomes diagonal |
| Iteration budget | Maximum number of update steps | Reaching it without the residual test passing is reported as failure |
How Eigenvalue Methods Works, Step by Step
Power iteration. Start with a nonzero vector and normalize it. At every check, compute its Rayleigh quotient and residual. If the residual passes, return the approximate pair. Otherwise calculate $y=Ax$, normalize $x\leftarrow y/\lVert y\rVert_2$, and repeat within the budget. An exactly zero product needs care: if the current vector lies in the nullspace, its zero eigenpair should already pass the residual test.
For symmetric $A$, expanding $x_0=\sum_i c_iv_i$ gives $A^kx_0=\sum_i c_i\lambda_i^kv_i$. Relative to a unique largest-magnitude eigenvalue, other components decay like $(\lambda_i/\lambda_1)^k$, provided $c_1\ne0$. This is a convergence condition, not a promise for every start or every spectrum.
Inverse iteration with a fixed shift. Choose finite $\mu$ such that $A-\mu I$ is nonsingular. Factor this matrix once with partial row pivoting. At each update, solve
$$ (A-\mu I)y=x, \qquad x\leftarrow y/\lVert y\rVert_2. $$
The transformed eigenvalues are $1/(\lambda_i-\mu)$. A uniquely nearest eigenvalue is therefore dominant in the transformed problem, assuming the starting component is present. The Netlib inverse-iteration reference explains this shift-and-solve construction and why an explicit inverse is unnecessary. Our code estimates the final eigenvalue with the Rayleigh quotient of original $A$, not by reporting a transformed inverse eigenvalue.
The packed LU helper swaps entire rows, including earlier stored multipliers, and records their final order. A shifted pivot at most $10^{-14}$ times the shifted matrix's initial maximum entry magnitude is rejected. A shift too close to an eigenvalue can trigger this threshold; changing the shift requires a new factorization.
Unshifted QR iteration. Initialize $T_0=A$ and $Z_0=I$. Repeatedly compute an unpivoted factorization and reverse its factors:
$$ T_k=Q_kR_k, \qquad T_{k+1}=R_kQ_k, \qquad Z_{k+1}=Z_kQ_k. $$
Since $R_k=Q_k^TT_k$, each update is an orthogonal similarity. For symmetric inputs, a diagonal limit exposes the eigenvalues, and accumulated columns of $Z_k$ supply their eigenvectors. Our code checks each accumulated column against original $A$ using its Rayleigh quotient, rather than relying on unchanged diagonal entries alone. On success, it sorts the eigenvalues and moves the corresponding eigenvector columns with them.
The QR helper uses Householder reflectors with an opposite-sign target to avoid cancellation. It skips an exactly zero tail, allowing a factorization even when the matrix is singular. Do not substitute column-pivoted least-squares QR and blindly use $RQ$: the additional permutation changes the similarity argument.
The Cornell numerical-analysis notes on QR iteration connect this process to orthogonal iteration and discuss the need for shifts and structural reduction in practical solvers. The simple unshifted version here can stagnate or cycle.
flowchart TD
A["Validate symmetric A, tolerance, and update budget"] --> B{"Method?"}
B -->|Power| C["Normalize starting vector"]
B -->|Inverse| D["Factor shifted matrix with partial row pivoting"]
D --> C
C --> E{"Original-matrix eigenpair residual passes?"}
E -->|Yes| F["Return approximate eigenpair"]
E -->|No| G{"Budget remains?"}
G -->|No| H["Report nonconvergence"]
G -->|Yes| I["Multiply or shifted-solve; normalize"]
I --> E
B -->|QR| J["Initialize T and accumulated orthogonal basis"]
J --> K{"All basis-column residuals pass?"}
K -->|Yes| L["Sort eigenpairs and verify orthogonality"]
K -->|No| M{"Budget remains?"}
M -->|No| H
M -->|Yes| N["Unpivoted QR; T becomes R Q; accumulate Q"]
N --> K
Worked Example with a Full Trace
Use
$$ A=\begin{pmatrix}2&1\\1&2\end{pmatrix}, \qquad x_0=(1,0)^T, \qquad \mu=3/2. $$
Direct multiplication verifies eigenpairs $\lambda_+=3$, $v_+=(1,1)^T/\sqrt2$ and $\lambda_-=1$, $v_-=(1,-1)^T/\sqrt2$. The Frobenius norm is $\sqrt{10}$.
For a complete short hand run, set relative tolerance $\tau=0.025$ and an update budget of $3$. The programs later use the tighter $10^{-10}$ tolerance. We show all three updates and all intermediate iterates for this hand run. Signs of eigenvectors are immaterial.
Power iteration. The first raw product is $(2,1)^T$, with norm $\sqrt5$, so $x_1=(2,1)^T/\sqrt5$. Then $Ax_1=(5,4)^T/\sqrt5$, giving $\theta_1=14/5$ and residual norm $3/5$. The next raw product normalizes to $x_2=(5,4)^T/\sqrt{41}$, and the third to $x_3=(14,13)^T/\sqrt{365}$.
Inverse iteration setup. The shifted matrix is
$$ B=A-(3/2)I=\begin{pmatrix}1/2&1\\1&1/2\end{pmatrix}. $$
Partial pivoting swaps its rows. Its packed factors represent lower rows $(1,0),(1/2,1)$ and upper rows $(1,1/2),(0,3/4)$. Solving $By=x_0$ gives $y_1=(-2/3,4/3)^T$, so the normalized vector is $u_1=(-1,2)^T/\sqrt5$. Reusing the same factors gives
$$ y_2=\frac{(10,-8)^T}{3\sqrt5}, \quad u_2=\frac{(5,-4)^T}{\sqrt{41}}, \qquad y_3=\frac{(-26,28)^T}{3\sqrt{41}}, \quad u_3=\frac{(-13,14)^T}{\sqrt{365}}. $$
The desired inverse eigenvalue of $B^{-1}$ is negative, so these vector signs alternate. Their directions converge toward the eigenspace for original eigenvalue $1$.
The following table records every vector state and check for the hand run. The absolute residual norms are the same for both methods; divide by $\sqrt{10}$ to obtain the last column.
| Update $k$ | Power vector $x_k$ | Power $\theta_k$ | Inverse vector $u_k$ | Inverse $\theta_k$ | Absolute residual norm | Relative residual |
|---|---|---|---|---|---|---|
| $0$ | $(1,0)^T$ | $2$ | $(1,0)^T$ | $2$ | $1$ | $0.316228$ |
| $1$ | $(2,1)^T/\sqrt5$ | $14/5$ | $(-1,2)^T/\sqrt5$ | $6/5$ | $3/5$ | $0.189737$ |
| $2$ | $(5,4)^T/\sqrt{41}$ | $122/41$ | $(5,-4)^T/\sqrt{41}$ | $42/41$ | $9/41$ | $0.069416$ |
| $3$ | $(14,13)^T/\sqrt{365}$ | $1094/365$ | $(-13,14)^T/\sqrt{365}$ | $366/365$ | $27/365$ | $0.023392$ |
Both vector methods first pass $\tau=0.025$ at update $3$. Their returned eigenvalues are approximations near $3$ and $1$, not exact answers merely because the check passed.
QR iteration. An unpivoted Householder factorization of $T_0=A$ gives
$$ Q_0=\frac1{\sqrt5}\begin{pmatrix}-2&1\\-1&-2\end{pmatrix}, \qquad R_0=\begin{pmatrix}-\sqrt5&-4/\sqrt5\\0&-3/\sqrt5\end{pmatrix}. $$
Reversing the factors and repeating produces all hand-run matrix states:
$$ T_1=\begin{pmatrix}14/5&3/5\\3/5&6/5\end{pmatrix}, \qquad T_2=\begin{pmatrix}122/41&9/41\\9/41&42/41\end{pmatrix}, \qquad T_3=\begin{pmatrix}1094/365&27/365\\27/365&366/365\end{pmatrix}. $$
Up to column signs, the accumulated basis is $Z_k=\begin{pmatrix}a_k&b_k\\b_k&-a_k\end{pmatrix}/\sqrt{a_k^2+b_k^2}$, with $(a_k,b_k)=(1,0),(2,1),(5,4),(14,13)$ for $k=0,1,2,3$. Thus all accumulated eigenvector states are specified too. Each column has the same relative residual listed in the table, so QR also passes the loose check after three updates. Every $T_k$ has trace $4$ and determinant $3$, preserving the spectrum.
Tighter program run. Let $q=3^k$. For these three runs in exact arithmetic, the residual norm is $2q/(q^2+1)$. At update $20$, the relative residual is about $1.81\times10^{-10}$; at update $21$, it is about $6.05\times10^{-11}$. Therefore the programs first pass $10^{-10}$ after $21$ updates. They scale $A$ by $2$ internally and scale the shift by the same amount; normalization and the relative residual are unchanged.
Trace It by Hand
Starting from $x_2=(5,4)^T/\sqrt{41}$, predict the next power direction before calculating its Rayleigh quotient. Reveal: the raw product is $(14,13)^T/\sqrt{41}$, so normalization gives $(14,13)^T/\sqrt{365}$ and the quotient is $1094/365$.
Now start power iteration with $(1,-1)^T$ instead. Check: the vector is already an eigenvector for eigenvalue $1$, with residual zero. The routine can return immediately even though $3$ is the dominant eigenvalue. Explain which convergence assumption the start violates.
Finally, change the inverse shift to $5/2$. Prediction: with a nonorthogonal start, inverse iteration now favors eigenvalue $3$, whose distance to the shift is $1/2$, rather than eigenvalue $1$, whose distance is $3/2$.
Implementation in Python, C++, and Java
Each panel is a complete program implementing power iteration, fixed-shift inverse iteration, and unshifted QR. Run Python 3 directly, compile C++14 or newer, or save Java as Main.java. Inputs are finite, nonempty, exactly symmetric square arrays. A nonzero start must have the matching length. Each call preserves its input arrays, validates the tolerance and iteration limit, and either returns a residual-checked approximation or raises an error.
prepare scales $A$ by its largest entry magnitude and scales the inverse shift consistently. The inverse routine factors once and reuses a packed, partially pivoted LU solve. Its pivot threshold is a separate fixed $10^{-14}$ policy; it is not the eigenpair stopping tolerance. The QR helper is deliberately unpivoted and supports exactly zero column tails. The QR eigenvalue routine returns ascending values and matching eigenvector columns. None of these choices guarantees convergence for every symmetric input.
The vector routines return an eigenvalue, normalized eigenvector, update count, and relative residual. QR returns all values, their vector matrix, its update count, and the largest relative residual. A zero update limit allows only an initial-state check. To make signs reproducible within a run, canonical makes the largest-magnitude component nonnegative, retaining the first on an exact tie. Different roundoff or sign conventions can still produce equivalent opposite eigenvectors on other inputs.
The sample prints a scientific residual rather than hiding it behind six decimal places. QR orthogonality is checked separately using the largest entry magnitude of $Z^TZ-I$; its printed zero is rounded, not a claim of exact floating-point orthogonality. Scaling reduces overflow risk but does not recover information lost through extreme dynamic range. These simple dense routines emphasize the algorithms and failure reporting rather than production performance.
All three programs print the same nine lines:
power lambda = 3.000000power vector = [0.707107, 0.707107]inverse lambda = 1.000000inverse vector = [-0.707107, 0.707107]QR lambdas = [1.000000, 3.000000]iterations = [21, 21, 21]max relative residual = 6.05e-11QR orthogonality error = 0.000000check = PASS
Complexity, Correctness, and Invariants
The number of iterations is part of the cost. There is no matrix-size-only worst-case convergence bound for these unrestricted elementary methods: a gap may be tiny, a target component may be missing, or an unshifted QR sequence may cycle. With a fixed budget $K$, the code terminates by returning a checked result or reporting failure.
Dense power iteration costs $\Theta(n^2)$ per update and residual check, with $O(n)$ vector workspace beyond its matrix copy. A fixed-shift inverse run has $\Theta(n^3)$ LU setup and $\Theta(n^2)$ work per solve/check. The simple QR implementation performs dense Householder factorization, dense products, and checks all $n$ candidate eigenvectors, giving $\Theta(n^3)$ per sweep and $\Theta(n^2)$ storage. For a fixed number of updates, existing zeros do not remove the dense matrix-product costs.
An initial exact eigenpair can return before any vector update; matrix validation and copying still cost $\Theta(n^2)$. Inverse iteration also prepares its shifted factorization before that initial check. QR checks an entire initial basis and therefore does cubic checking work even when it returns at update zero. For a successful run using $k$ updates, the bounds below include setup and checks.
| Method | Dense total time for $k$ updates | Stored data and workspace | Target and convergence conditions |
|---|---|---|---|
| Power iteration | $\Theta((k+1)n^2)$ | $\Theta(n^2)$ matrix copy plus $O(n)$ vectors | One eigenpair; a unique largest magnitude and nonzero target component give the standard convergence argument |
| Fixed-shift inverse iteration | $\Theta(n^3+(k+1)n^2)$ | $\Theta(n^2)$ matrix and LU factors | One nearby eigenpair; shifted solve must pass pivots and the nearest eigenspace must be represented |
| Unshifted dense QR here | $\Theta((k+1)n^3)$ | $\Theta(n^2)$ matrices | All symmetric eigenpairs if convergence occurs; equal magnitudes can obstruct it |
| Structured QR with shifts and deflation | Typically cubic dense reduction plus cheaper structured sweeps | $O(n^2)$ if eigenvectors are retained | Practical eigensolvers exploit Hessenberg or symmetric tridiagonal form |
| Sparse Lanczos for symmetric matrices | Depends on matrix-vector products, subspace size, and reorthogonalization | Sparse matrix plus a retained basis | Selected eigenpairs or invariant subspaces of large sparse problems |
For power iteration, the invariant is $x_k$ proportional to $A^kx_0$, as long as normalization is defined. Orthogonal eigenbasis expansion then proves component suppression under the stated gap and starting-vector assumptions. The Rayleigh quotient of a normalized vector gives a scalar to test, and the residual verifies the resulting pair.
For inverse iteration, the corresponding invariant is $u_k$ proportional to $(A-\mu I)^{-k}u_0$. Its eigenvectors are unchanged, and the transformed eigenvalues are reciprocal distances with sign. Under a unique nearest eigenvalue and a nonzero component in that direction, unwanted components decay at a rate governed by $\lvert(\lambda_{\mathrm{target}}-\mu)/(\lambda_{\mathrm{next}}-\mu)\rvert$.
For QR, induction gives
$$ T_k=Z_k^TAZ_k, \qquad Z_k^TZ_k=I. $$
The update preserves this identity because $T_{k+1}=Q_k^TT_kQ_k$ and $Z_{k+1}=Z_kQ_k$. If $T_k$ becomes diagonal, its diagonal entries and columns of $Z_k$ form exact eigenpairs of original $A$. The scaled implementation obeys the corresponding identity with $A/s$, then restores eigenvalue units.
For a real symmetric matrix, any approximate pair satisfies
$$ \min_i\lvert\theta-\lambda_i\rvert \le\frac{\lVert Ax-\theta x\rVert_2}{\lVert x\rVert_2}. $$
This follows by expanding the residual in the orthonormal eigenbasis. It locates some nearby eigenvalue, but does not certify which eigenvalue was found or a unique eigenvector when eigenvalues cluster. Floating-point orthogonality and similarity drift are additional reasons to check residuals against the original input.
Power and inverse iteration preserve their respective repeated-operator directions; QR preserves an orthogonal similarity and accumulates its change of basis. These invariants explain what a converged state means, while the spectral gap and starting components determine whether the intended target is reached.
Variants, Pitfalls, and Optimizations
A valid eigenpair may be the wrong target. Starting the sample power routine in direction $(1,-1)^T$ returns eigenvalue $1$ with zero residual. That is a correct eigenpair but not the dominant one. Multiple starts or subspace methods help when a single starting direction is unreliable. Residuals certify the equation, not the requested ordering of the spectrum.
Sign changes and equal magnitudes. A negative dominant eigenvalue can make successive vectors alternate sign, so $\lVert x_{k+1}-x_k\rVert$ is a poor sole stopping test. With $A=\begin{pmatrix}0&1\\1&0\end{pmatrix}$ and start $(1,0)^T$, power alternates coordinate axes. The Rayleigh quotient remains $0$ while the residual norm remains $1$. Unshifted QR also fails to diagonalize this example; an unchanged diagonal is not convergence. The budget must report failure.
A shift changes the problem being solved internally. Do not choose an exact eigenvalue as a fixed inverse shift: the shifted matrix is singular. Near-eigenvalue shifts can speed direction convergence but make the solves delicate or trigger pivot rejection. Rayleigh quotient iteration updates the shift and can converge rapidly near a simple symmetric eigenpair, but generally needs new shifted factorizations. An inverse estimate is not itself an eigenvalue of original $A$.
QR needs similarity, shifts, and deflation. A shifted sweep factors $T_k-\mu_kI=Q_kR_k$, then forms $T_{k+1}=R_kQ_k+\mu_kI$. Practical methods choose shifts and split converged blocks using deflation. Symmetric solvers usually reduce to tridiagonal form first; general real solvers work toward a real Schur form that may contain $2\times2$ blocks for complex conjugate pairs. Columns of an accumulated nonsymmetric Schur basis are not generally individual eigenvectors. The code's symmetric restriction avoids that false interpretation.
Repeated and clustered eigenvalues. Eigenvectors within a repeated eigenspace are not unique. Even a very small residual may leave the individual direction poorly determined when nearby eigenvalues are tightly clustered. Compare invariant subspaces, residuals, and orthogonality rather than requiring one particular signed vector. Arnoldi, Lanczos, block iteration, and restart strategies extend the basic repeated-multiplication idea for larger problems.
Exploit structure without changing the target. Sparse power iteration can use a sparse matrix-vector product rather than dense arrays. Sparse shift-invert methods reuse sparse factors or suitable solves, while fill-in affects memory. For covariance directions, applying an SVD to centered data can avoid explicitly forming a covariance matrix. For graph Laplacians, the smallest eigenvalue may be a known zero; the next-smallest direction requires excluding that nullspace or using an appropriate shifted/subspace method.
Always check original-matrix residuals and enforce an update limit; unchanged Rayleigh quotients or diagonals can hide cycling. Track eigenvector signs and eigenvalue targets separately, and use unpivoted QR in the similarity iteration rather than blindly reusing column-pivoted least-squares QR.
Practice Problems and Real-World Use
Validate the methods on matrices whose spectrum is known independently. For example, construct $A=V\Lambda V^T$ from an orthogonal matrix and a chosen diagonal spectrum, then verify returned eigenvalues, residuals, and QR orthogonality. For a simple symmetric eigenvalue, compare vectors up to sign. For repeated eigenvalues, compare their spanned subspace instead.
Include negative eigenvalues, a zero matrix, a singular but valid QR input, a start orthogonal to the dominant direction, an exact singular shift, a shift needing an LU row swap, and the equal-magnitude cycling example. Repeat tests after uniform scaling and preserve the original input for checks. A nonconvergence test is a required outcome for a limited method, not automatically an implementation defect.
| Exercise or application | Concrete task or checkpoint | Idea tested | Difficulty |
|---|---|---|---|
| Hand-trace comparison | Reproduce all three sample methods through update $3$ at $\tau=0.025$ | Different operators and the same residual criterion | Beginner |
| Negative dominant eigenvalue | Use $\operatorname{diag}(-4,1)$ with both starting components nonzero | Magnitude dominance and sign alternation | Intermediate |
| Wrong-target regression | Start the sample with $(1,-1)^T$ and obtain eigenvalue $1$ immediately | A zero residual does not certify dominance | Intermediate |
| Cycling regression | Use rows $(0,1),(1,0)$ with start $(1,0)^T$ and expect a budget failure | Equal magnitudes and misleading unchanged quotients | Intermediate |
| Shift targeting | Change the sample shift from $3/2$ to $5/2$; switch from target $1$ to $3$ | Closest-eigenvalue selection and factor reuse | Intermediate |
| Principal component direction | Find a largest covariance eigendirection, or compare with an SVD of centered data | Dominant variance direction | Applied |
| Graph spectral partitioning | Find an appropriate small Laplacian eigenspace after handling the constant nullspace | Targeting and structural zero modes | Advanced |
| Vibration and stability models | Analyze real symmetric mode matrices or the appropriate generalized problem | Eigenmodes, units, and modeling assumptions | Applied |