Eigenvalue Decomposition

Eigenvalue Decomposition (cross-cutting/numerics/eigenvalue-decomposition)

Use when the task is computing eigenvalues and eigenvectors of a square matrix with the Python standard library only (no numpy). The leaf covers the power iteration for the dominant eigenvalue- eigenvector pair of a general square matrix, deflation for the next pairs, the Jacobi eigenvalue algorithm for the full spectrum of a real symmetric matrix, the convergence criteria of both methods, and the residual A v - lambda v that verifies every returned pair.

Domain quick reference

• Power iteration finds the dominant pair. Start from a unit vector v, multiply by A, normalize to unit 2-norm, repeat; the Rayleigh quotient lam = v^T A v estimates the eigenvalue at each iterate. Converges when the pair residual ||A v - lam v||_inf falls at or below the tolerance. Worked anchor: A = [[2, 0], [0, 1]] converges to lam = 2.0 with v = [1, 0]; the quotient sequence 1.5, 1.8, 1.95, ... settles on 2.0.
• Convergence criteria: power iteration stops on the pair residual (the Rayleigh quotient change between iterates is the monitor, but the quotient converges quadratically while the eigenvector converges only linearly, so a quotient-change-only stop returns a vector accurate to only about sqrt(tol)). Jacobi stops when the off-norm, the square root of the sum of squared off-diagonal entries, falls at or below the sweep tolerance.
• Deflation exposes the next eigenvalue. Subtract the found pair with the rank-one update A - lam * v v^T (v unit norm); the remaining eigenvalues and eigenvectors are unchanged, so power iteration on the deflated matrix yields the next dominant pair. Worked anchor: deflating [[2, 0], [0, 1]] by (2.0, [1, 0]) gives [[0, 0], [0, 1]], whose dominant eigenvalue 1.0 is the second eigenvalue of A. Call power_spectrum(A, count) to get the count largest pairs with the deflation built in.
• Jacobi diagonalizes a real symmetric matrix. Each sweep visits every off-diagonal pair (p, q) and applies the plane rotation that zeroes A[p][q]; rotations accumulate in the eigenvector matrix. Worked anchor: [[2, 1], [1, 2]] gives eigenvalues 3.0 and 1.0 with unit-norm eigenvectors [0.707, 0.707] and [0.707, -0.707], and the residual A v - lam v is zero to machine precision. A = [[4, 1, 1], [1, 4, 1], [1, 1, 4]] gives 6.0, 3.0, 3.0 (6.0 belongs to the all-ones vector).
• Normalization conventions: eigenvectors are returned with unit 2-norm, eigenvalues sorted in descending order (largest first). The sign of an eigenvector is arbitrary: v and -v are the same eigenvector, so compare results with the residual, not with a memorized sign.
• Residual verification: residual_norm(A, v, lam) = max_i |sum_j A[i][j] * v[j] - lam * v[i]|, the natural check that (lam, v) is a true pair. For a unit-norm v it is the absolute defect; a value at machine-epsilon scale confirms the pair. Worked anchor: [[2, 0], [0, 1]] with (2.0, [1, 0]) gives 0.0; with (1.0, [1, 0]) it gives 1.0.
• Domain fit: symmetric matrices arise as covariance, stiffness, and graph Laplacian matrices; the Jacobi path fits them. General square matrices (including non-symmetric ones) use the power iteration path, which returns only the count largest pairs, not the full spectrum. Cost: power iteration is O(n^2) per iterate; Jacobi is O(n^3) per sweep; both are dense-matrix tools.

Workflow

1. Confirm the matrix is square, numeric, and non-empty; the logic validates this and raises ValueError.
2. Pick the method by the matrix and the need: a symmetric matrix needing its full spectrum goes to jacobi_eigen(A); a general square matrix needing the dominant pair or a few top pairs goes to power_iteration(A) or power_spectrum(A, count).
3. For the single dominant pair, call power_iteration(A) -> (lam, v); pass v0 when the default all-ones start is orthogonal to the dominant eigenspace.
4. For the count largest pairs of a general matrix, call power_spectrum(A, count) -> [(lam, v), ...]; the module deflates between pairs with a deterministic mixed start vector.
5. For the full spectrum of a symmetric matrix, call jacobi_eigen(A) -> (eigenvalues, eigenvectors) with eigenvalues descending and unit-norm eigenvector columns.
6. Verify every returned pair with residual_norm(A, v, lam); a residual at machine-epsilon scale confirms the pair, a larger residual means the iteration did not converge, the wrong method was used, or the tolerance was too loose.
7. Report the eigenvalues, the eigenvectors, the convergence criterion used, and the residuals so the result is checkable.

Pitfalls

• Confusion with matrix-operations: solving Ax = b, the determinant, and the inverse are different objects from the spectrum. An eigenvalue pair satisfies A v = lam v; a solve returns x for one right-hand side b. A zero eigenvalue means A is singular (det 0, no inverse), which is the only overlap: read it off the spectrum, but route "solve the linear system, determinant, inverse" tasks to matrix-operations and "eigenvalues, eigenvectors, spectrum" tasks here.
• Confusion with convergence-verification: that leaf checks whether a CFD or structural mesh sequence has converged with Richardson extrapolation and the grid convergence index. This leaf's "convergence" is the internal stopping criterion of the iterative algorithms (pair residual, off-diagonal sum). Route "grid convergence index, refinement study, discretization error" tasks there.
• Confusion with finite-difference-derivatives: differentiating a function at a point estimates a derivative, not an eigenvalue. A Jacobian matrix of partial derivatives is a matrix you might later decompose, but "compute the Jacobian" tasks route to finite-difference-derivatives; "eigenvalues of the linearized system" tasks route here.
• Confusion with ode-solvers: time marching dy/dt = f(t, y) steps the state forward; it does not decompose a matrix. Eigenvalue analysis of the linearized system matrix is a stability and analysis task, not an integration. Route integration tasks to ode-solvers; route spectrum-of-the-state-matrix tasks here.
• Jacobi is only for symmetric matrices: feeding a non-symmetric matrix to jacobi_eigen raises ValueError. A non-symmetric matrix (or one that just needs a few top pairs) takes the power iteration path.
• Power iteration finds only the dominant pair: the sequence converges to the largest-magnitude eigenvalue. Additional pairs need deflation (power_spectrum) or Jacobi; equal-magnitude eigenvalues (for example 3.0 and -3.0) converge slowly or mix.
• Start vector orthogonality: a v0 lying exactly in another eigenspace returns that pair instead, because the iteration cannot leave the eigenspace. The all-ones default can hit this on deflated matrices; power_spectrum avoids it with a mixed start, and a direct power_iteration call should pass a custom v0 and re-verify with the residual.
• Rayleigh quotient vs eigenvector accuracy: the quotient converges quadratically while the eigenvector converges linearly, so a stop based only on quotient change can return an eigenvector accurate to only about sqrt(tol). The residual-based stop in the logic module avoids this; keep the residual check in your report.
• Loose tolerance: a large tol stops the iteration early and the pair has a small-but-nonzero residual. Tighten tol until the residual is at machine-epsilon scale for the analysis.
• Sign and scale: eigenvectors are defined up to sign and scale; this leaf returns unit 2-norm vectors. Compare pairs via the residual, never via a memorized sign.
• Non-square or ragged input: both routines validate and raise ValueError; validate before calling rather than relying on downstream errors.
• Stdlib-only contract: the logic module must stay free of numpy and other third-party imports; the gate 3 contract test enforces stdlib-only imports.

Behavior contract (gate 3)

The Jacobi, power iteration, deflation, and residual logic is exercised by the gate 3 contract test: scripts/test_eigenvalue_logic.py against scripts/eigenvalue_logic.py (stdlib unittest, offline). Run:

python3 scripts/test_eigenvalue_logic.py

Anchors: jacobi_eigen([[2, 1], [1, 2]]) returns eigenvalues 3.0 and 1.0 with unit-norm eigenvectors satisfying A v = lam v to 1e-8; [[4, 1, 1], [1, 4, 1], [1, 1, 4]] returns 6.0, 3.0, 3.0; power_iteration([[2, 0], [0, 1]]) returns 2.0 with v = [1, 0]; after deflation the second eigenvalue is 1.0; power_spectrum matches jacobi_eigen on the 2x2 anchor. Trend properties: Jacobi eigenvectors are orthonormal, eigenvalues sort descending, and residual_norm is 0.0 for an exact pair. ValueError on non-square, ragged, non-numeric, or bool input, and on non-symmetric input to the Jacobi routine.

Compliance

• NACA Report 824 is US government work (public domain); the pack anchor per standards-map.yaml. Jacobi's method (1846) and the power iteration are classical numerical linear algebra, generic textbook methodology, not RTCA, SAE, or IAQG content; summary and formulas only.
• compliance: STANDARDS-REF, gated: false.