Eigenvalue Problems — Math Reference

Eigenvalue Problems

1. At a glance

The eigen-decomposition is one of the most consequential ideas in applied mathematics. Given a square matrix A, the question “for which vectors v does A act as scalar multiplication?” — i.e. Av = λv — exposes the intrinsic axes along which the linear map stretches space. Once you know those axes, every question about powers A^k, exponentials exp(tA), inverses, resolvents, and long-time behavior becomes algebraic instead of analytic.

Eigenvalue problems are everywhere:

• Vibration analysis — natural frequencies and mode shapes of structures.
• PCA — directions of maximum variance in data.
• PageRank — the dominant eigenvector of a stochastic web-link matrix.
• Quantum mechanics — eigenvalues of the Hamiltonian are energy levels.
• Stability of dynamical systems — sign of Re(λ) controls growth/decay.
• Network analysis — eigenvector centrality, spectral graph theory.
• Markov chains — stationary distribution is the left eigenvector with λ = 1.
• Graph Laplacians — smallest eigenvalues reveal cluster structure.
• Spectral methods (PDEs) — separation of variables uses eigenfunctions.
• Control theory — pole placement, Lyapunov equations, controllability gramians.
• Machine learning — kernel PCA, diffusion maps, Laplacian eigenmaps, spectral clustering, normalized cuts.

This note covers the formulations, the algorithms used to solve them in practice, the software stack, and the numerical pitfalls.

2. Standard eigenvalue problem

For an n × n matrix A (real or complex), a scalar λ and nonzero vector v form an eigenpair iff:

    A v = λ v,   v ≠ 0

Equivalently (A − λI) v = 0, which has a nonzero solution iff A − λI is singular, iff its determinant vanishes:

    p_A(λ) = det(A − λI) = 0          (characteristic polynomial)

p_A is a degree-n polynomial in λ with leading coefficient (−1)^n. By the fundamental theorem of algebra it has exactly n roots in C counted with multiplicity.

Algebraic vs geometric multiplicity

• Algebraic multiplicity m_a(λ) = multiplicity of λ as a root of p_A.
• Geometric multiplicity m_g(λ) = dim ker(A − λI) = number of linearly independent eigenvectors associated with λ.

Always 1 ≤ m_g(λ) ≤ m_a(λ). When m_g < m_a the matrix is defective: it cannot be diagonalized over C. The canonical example is the Jordan block

    J = [[λ, 1], [0, λ]]

with m_a = 2, m_g = 1.

Diagonalization

If A has n linearly independent eigenvectors (i.e. m_g = m_a for every λ), stack them as columns of V and the eigenvalues on the diagonal of Λ = diag(λ_1, ..., λ_n). Then

    A V = V Λ        →        A = V Λ V^{-1}

A is diagonalizable. This is the cleanest decomposition; almost all of practical linear algebra rests on the cases where it exists or has a close substitute (Schur form, SVD).

Trace and determinant

    tr(A) = Σ λ_i ,    det(A) = Π λ_i

— useful sanity checks after a numerical computation.

3. Generalized eigenvalue problem

Many physical formulations land on:

    A v = λ B v

where B is also n × n. The pair (A, B) is called a matrix pencil and its spectrum is the set of λ for which A − λB is singular.

When B is invertible the problem reduces to standard via B^{-1} A, but doing the inversion explicitly destroys symmetry and conditioning, so dedicated algorithms (QZ for nonsymmetric, dsygv for symmetric-definite) are preferred.

When B is singular the pencil is singular; some eigenvalues become infinite (degeneracy in B) and others may be undefined. The complete canonical form is the Kronecker canonical form.

Where it appears

• Vibration modal analysis — undamped free vibration M ẍ + K x = 0 yields (K − ω² M) v = 0, i.e. K v = ω² M v with mass M and stiffness K. Eigenvalues are squared natural frequencies; eigenvectors are mode shapes.
• Structural eigenvalue buckling — (K − λ K_g) v = 0; the smallest positive λ is the critical load factor.
• Generalized PCA / canonical correlation analysis — Σ_xy Σ_yy^{-1} Σ_yx v = λ Σ_xx v.
• Fisher linear discriminant analysis — between-class scatter S_B and within-class scatter S_W give S_B v = λ S_W v.
• Quantum chemistry (generalized eigenvalue Hartree-Fock / DFT) — F C = S C ε with overlap matrix S.

4. Symmetric / Hermitian case

A is symmetric if A = A^T (real entries) or Hermitian if A = A* (complex). The spectral theorem is the cornerstone:

If A is Hermitian then all eigenvalues are real, eigenvectors corresponding to distinct eigenvalues are orthogonal, and there exists a unitary Q such that

    A = Q Λ Q*           (Λ real diagonal, Q* Q = I)

For real symmetric A, Q is orthogonal: A = Q Λ Q^T.

Consequences worth memorizing:

• A is always diagonalizable (no defective Hermitian matrix exists).
• The eigenvalues are the stationary values of the Rayleigh quotient (§9).
• A is positive (semi)definite iff all eigenvalues are > 0 (≥ 0).
• ‖A‖_2 = max |λ_i|, cond_2(A) = |λ_max|/|λ_min| for SPD.

This is the most common case in physics, engineering, statistics, and ML — and it admits the most efficient and numerically stable algorithms.

5. Symmetric eigensolvers

The standard dense pipeline for n × n symmetric A:

1. Tridiagonalize via a sequence of Householder reflectors: A → T = Q^T A Q with T symmetric tridiagonal. Cost: (4/3) n³ flops.
2. Diagonalize T — three competing algorithms:
   • Symmetric QR (with Wilkinson shifts) — O(n²) per sweep, robust.
   • Divide and conquer (Cuppen 1981) — O(n²·³) empirical, with deflation often much faster; parallelizes well.
   • MRRR — Multiple Relatively Robust Representations (Dhillon and Parlett 2004) — O(n²) to compute all eigenpairs without re-orthogonalization; basis of LAPACK’s dsyevr.

LAPACK routines

• dsyev — basic QR-based driver.
• dsyevd — divide-and-conquer driver (good when you want all eigenpairs).
• dsyevr — MRRR driver (most efficient; supports subset by index or by interval).
• dsyevx — bisection + inverse iteration (legacy; superseded by dsyevr).

Cost

• All n eigenvalues only: O(n³) with small constant.
• All eigenpairs (vectors too): O(n³).
• k extremal eigenpairs via Lanczos: O(n²k) if A is dense, or O(nnz · k · m) for sparse A where m is the Krylov dimension.

6. Non-symmetric eigenproblem

For general A ∈ C^{n×n}:

• Eigenvalues may be complex even for real A (occurring in conjugate pairs).
• Eigenvectors may be linearly dependent — A may be defective.
• The right place to land numerically is the Schur decomposition:

    A = Q T Q*

with Q unitary and T upper triangular (block triangular with 2×2 blocks for real Schur). The Schur form always exists, eigenvalues appear on the diagonal, and Schur is numerically far better behaved than attempting to compute the Jordan form.

LAPACK routines

• dgehrd — Hessenberg reduction A → H = Q^T A Q.
• dhseqr — QR algorithm on Hessenberg H to get Schur T.
• dgees — full Schur decomposition driver.
• dgeev — driver for eigenvalues and (left/right) eigenvectors.
• dggev — generalized nonsymmetric eigenproblem via QZ.

Eigenvector reconstruction

Once T is computed, eigenvectors are obtained by back-substitution on triangular systems (T − λI) x = 0. This step can be ill-conditioned when eigenvalues cluster.

7. Power method

The simplest iterative scheme:

    v_{k+1} = A v_k / ‖A v_k‖

Under mild conditions (dominant eigenvalue |λ_1| > |λ_2|, starting vector with nonzero λ_1-component) the iterate converges to the dominant eigenvector. The convergence rate is linear with ratio |λ_2 / λ_1|.

Estimate the eigenvalue via the Rayleigh quotient λ ≈ v_k^T A v_k.

Why it matters

• PageRank — the dominant left eigenvector of the column-stochastic web matrix. One iteration is a sparse matrix-vector product, perfect for the web.
• Markov chain stationary distribution — largest eigenvalue is 1, the eigenvector gives the stationary distribution.
• Dominant vibration mode — quick estimate of the lowest frequency.
• Spectral radius estimation — useful in numerical analysis of iterative methods.

Variants

• Subspace iteration — power method on a block of vectors with periodic orthogonalization; gives the top-k eigenpairs.
• Shifted power iteration — apply (A − σI) to suppress chosen eigenvalues.

8. Inverse iteration

For an approximate eigenvalue μ ≈ λ, iterate:

    Solve  (A − μI) w = v_k
    v_{k+1} = w / ‖w‖

This is the power method applied to (A − μI)^{-1}, whose dominant eigenvalue is 1 / (λ − μ) — huge when μ is close to λ. Convergence is therefore very fast (one or two steps once you’re near).

The linear solve (A − μI) w = v_k is dangerous-looking because A − μI is nearly singular, but this is precisely what makes inverse iteration work; the near-null direction is exactly the eigenvector we want.

Used to refine an eigenpair after a global solver has located μ, and as the workhorse inside many subspace methods.

9. Rayleigh quotient

For Hermitian A:

    R(v) = (v* A v) / (v* v)

Properties:

• R(v) = λ when v is an eigenvector with eigenvalue λ.
• min_v R(v) = λ_min, max_v R(v) = λ_max.
• Courant-Fischer min-max theorem — the k-th eigenvalue (sorted) equals

    λ_k = min_{S, dim S = k}  max_{0≠v∈S}  R(v)
        = max_{S, dim S = n−k+1}  min_{0≠v∈S}  R(v)

This characterization is the key tool for proving interlacing theorems (Cauchy interlacing, Weyl inequalities) and underlies many iterative algorithms.

Rayleigh-quotient iteration

Combine §8 and §9: at each step, set μ_k = R(v_k) and solve (A − μ_k I) w = v_k. For Hermitian matrices the convergence is cubic in the neighborhood of an eigenvalue — extraordinarily fast.

10. QR algorithm

Discovered independently by Francis (1961) and Kublanovskaya (1961). The unshifted version is shockingly simple:

    A_0 = A
    Repeat:
        Q_k R_k = A_k        (QR factorization)
        A_{k+1} = R_k Q_k    (multiply back in reverse)

A_{k+1} = Q_k^T A_k Q_k, so the sequence is a similarity transformation — eigenvalues are preserved. Under mild conditions the sequence converges to upper triangular (Schur) form, with eigenvalues on the diagonal in order of decreasing magnitude.

Practical implementations add three crucial ingredients:

1. Hessenberg reduction — first reduce A to upper Hessenberg H (symmetric → tridiagonal) via Householder reflectors. The QR algorithm preserves Hessenberg structure, cutting the per-iteration cost from O(n³) to O(n²).
2. Shifts — replace A_k with A_k − μ_k I before the factorization, then undo. Single-shift (Rayleigh shift) or double-shift (Wilkinson, Francis) give cubic convergence for symmetric matrices and ensure complex conjugate pairs are found in real arithmetic.
3. Deflation — once a subdiagonal entry becomes negligibly small, split the problem.

The QR algorithm is the standard dense eigensolver for both symmetric and nonsymmetric problems. Total cost O(n³).

11. Lanczos algorithm

Lanczos (1950) attacked the large sparse symmetric eigenproblem with a Krylov-subspace approach. Given a starting vector v_1, the m-step Lanczos process builds an orthonormal basis V_m = [v_1, ..., v_m] of the Krylov subspace

    K_m(A, v_1) = span{ v_1, A v_1, A² v_1, ..., A^{m−1} v_1 }

via a three-term recurrence:

    β_{j+1} v_{j+1} = A v_j − α_j v_j − β_j v_{j−1}
    α_j = v_j^T A v_j
    β_{j+1} = ‖A v_j − α_j v_j − β_j v_{j−1}‖

The projection T_m = V_m^T A V_m is symmetric tridiagonal, and its eigenvalues (Ritz values) approximate eigenvalues of A. Extreme eigenvalues converge fastest, often after m ≪ n iterations.

Why it is sparse-friendly

The only operation involving A is the matrix-vector product A v_j. If A has nnz nonzeros, that costs O(nnz). The rest of the iteration is O(n) per step. Total O(nnz · m + n m²) for m Lanczos steps.

Loss of orthogonality

In floating-point arithmetic, V_m rapidly loses orthogonality once Ritz values converge. Remedies: full reorthogonalization, selective reorthogonalization (Parlett-Scott 1979), partial reorthogonalization (Simon 1984). Modern implementations restart implicitly when memory is bounded.

12. Arnoldi algorithm

For nonsymmetric A, the analog is Arnoldi (1951): build an orthonormal basis of K_m(A, v_1) with a (full) Gram-Schmidt instead of a three-term recurrence because V_m^T A V_m is upper Hessenberg, not tridiagonal:

    h_{i,j} = v_i^T A v_j  for i ≤ j+1
    H_m = V_m^T A V_m   (upper Hessenberg, (m+1)×m for the residual form)

Eigenvalues of H_m are the Ritz approximations.

ARPACK

Lehoucq, Sorensen, and Yang (1998) packaged the implicitly restarted Arnoldi (IRAM) method as the ARPACK library. Symmetric problems use the equivalent implicitly restarted Lanczos. ARPACK is the engine behind:

• MATLAB eigs
• SciPy scipy.sparse.linalg.eigs (nonsymmetric) and eigsh (symmetric)
• Octave eigs
• Many other front-ends

IRAM keeps Krylov dimension bounded (e.g. 20-50) while extracting k desired eigenpairs by repeatedly compressing the basis.

Shift-invert

For interior or smallest-magnitude eigenvalues, replace the operator with (A − σI)^{-1}. Largest eigenvalues of this operator are eigenvalues of A nearest σ. ARPACK’s “shift-invert mode” requires a sparse factorization of A − σI, then each “matvec” is a triangular solve.

13. Jacobi-Davidson

Sleijpen and van der Vorst (1996) introduced an iterative method targeted at interior eigenvalues, where Krylov methods converge slowly. At each step it:

1. Picks a Ritz value θ and Ritz vector u from the current subspace.

2. Solves the correction equation approximately:

       (I − u u*) (A − θ I) (I − u u*) t = −r,    where r = A u − θ u
   

3. Expands the subspace by t and repeats.

The correction equation is solved inexactly (e.g. a few GMRES steps), making the method inner-outer iterative. JD handles generalized eigenproblems and interior eigenvalues naturally.

The PRIMME library (Stathopoulos and McCombs 2010) implements Jacobi-Davidson plus several near-related methods and is one of the strongest options for hard eigenvalue problems.

14. Sparse eigenvalue libraries

Library	Strength	Language / interface
ARPACK / ARPACK-NG	IRAM, mature, ubiquitous	Fortran 77, bindings in MATLAB / SciPy / Julia
SLEPc	PETSc-based, parallel, generalized eigenvalues, polynomial / nonlinear eigenvalue problems	C, MPI
PRIMME	Jacobi-Davidson, GD+k, interior eigenvalues, preconditioned	C, with MATLAB / Python wrappers
FEAST (Polizzi 2009)	Contour-integration eigensolver — projects onto eigenspace of an interval via numerical contour integral; embarrassingly parallel	Fortran / C; in Intel MKL Extended Eigensolvers
Anasazi (Trilinos)	Block Krylov-Schur, block Davidson, LOBPCG	C++, parallel
LOBPCG (Knyazev 2001)	Locally Optimal Block Preconditioned CG; symmetric only, very effective with a preconditioner	In SLEPc, BLOPEX, SciPy
MAGMA	LAPACK-equivalent for GPUs	CUDA / HIP
cuSOLVER	NVIDIA’s dense and sparse eigensolvers	CUDA

For nonlinear eigenproblems T(λ) v = 0 (e.g. quadratic (λ² M + λ C + K) v = 0 for damped vibration), SLEPc’s NEP module and the NLEIGS algorithm (Güttel et al. 2014) are state of the art.

15. Generalized eigenvalue solvers

For the pair (A, B):

Case	LAPACK routine	Method
A symmetric, B SPD	dsygv, dsygvd, dsygvx	Cholesky B = L L^T, solve (L^{-1} A L^{-T}) y = λ y
A, B general	dggev	QZ algorithm (generalized Schur)
A, B Hermitian, B HPD	zhegv	Cholesky reduction

The Cholesky reduction trick: when B is SPD, factor B = L L^T, define C = L^{-1} A L^{-T} (still symmetric), and solve the standard problem C y = λ y. Recover original eigenvectors via v = L^{-T} y. This preserves symmetry — crucial for numerical accuracy.

For sparse generalized problems use SLEPc, ARPACK with shift-invert mode, or PRIMME — all support (A, B) natively.

16. Singular (defective) eigenproblem

A defective matrix has fewer linearly independent eigenvectors than its dimension. Examples appear:

• Optimization — saddle Hessians with repeated eigenvalues.
• Control — uncontrollable / unobservable modes.
• Mechanics — degenerate vibration modes (square plate, symmetric drum).

The exact decomposition is the Jordan canonical form A = P J P^{-1} with J block diagonal, each block of the form

    J_k(λ) = [[λ, 1, 0, ...],
              [0, λ, 1, ...],
              ...
              [0, 0, ..., λ]]

The Jordan form is infinitely numerically unstable: an arbitrarily small perturbation of A makes it diagonalizable. Therefore:

Do not compute the Jordan form numerically. Use Schur instead, and treat near-defective behavior with care (rank-revealing factorizations, deflation by tolerance).

17. Condition number

For a general matrix, the 2-norm condition number is

    κ_2(A) = σ_max(A) / σ_min(A)        (singular values, from the SVD)

For SPD A, singular values equal eigenvalues, so

    κ_2(A) = λ_max / λ_min

κ_2 predicts numerical sensitivity: a perturbation in b of relative size ε produces a perturbation in the solution of A x = b of relative size up to κ_2 · ε. Eigenvalue computations themselves are well-conditioned for symmetric matrices (perturbations stay bounded by Weyl’s inequality |λ_i(A + E) − λ_i(A)| ≤ ‖E‖_2) but ill-conditioned for non-normal matrices — this is the realm of pseudospectra (Trefethen and Embree 2005).

18. Stability of dynamical systems

Continuous-time linear system:

    ẋ = A x

Solutions are linear combinations of exp(λ_i t) v_i. Asymptotic stability (x(t) → 0 for every initial condition) iff every eigenvalue of A has Re(λ) < 0.

Discrete-time system x_{k+1} = A x_k is stable iff every eigenvalue has |λ| < 1 (spectral radius < 1).

For nonlinear systems near a fixed point, the same condition applied to the Jacobian gives local asymptotic stability (Hartman-Grobman). This is the foundation of:

• ODE numerical stability analysis ([[Math/ode-numerical-methods]]).
• Structural dynamics (vibration damping requires Re(λ) < 0).
• Control synthesis (place poles in left half-plane).
• Reinforcement learning value iteration convergence.
• Lyapunov methods — A^T P + P A = −Q solvable with P ≻ 0, Q ≻ 0 iff A is stable.

19. Applications

19.1 Vibration modal analysis

Free undamped vibration of an n-DoF system:

    M ẍ + K x = 0

Ansatz x(t) = v · cos(ω t) gives:

    K v = ω² M v

— a symmetric generalized eigenvalue problem. Eigenvalues ω_i² are squared natural frequencies, eigenvectors v_i are mode shapes. With M-orthogonality v_i^T M v_j = δ_{ij}, the modal coordinates decouple the equations of motion. See [[Engineering/vibration-dynamics]] and [[Engineering/structural-dynamics]].

For damped vibration M ẍ + C ẋ + K x = 0 the eigenproblem becomes quadratic (λ² M + λ C + K) v = 0 — handled by linearization to a 2n × 2n standard generalized eigenproblem, or by direct nonlinear-eigenvalue solvers.

19.2 Principal Component Analysis

Given centered data X ∈ R^{n × p}, the sample covariance is Σ = X^T X / (n − 1). Its eigen-decomposition

    Σ = V Λ V^T

gives principal directions (columns of V) sorted by variance (diagonal of Λ). Equivalent and numerically preferred: SVD of X directly. See [[Math/svd-pca-spectral]].

19.3 PageRank

Define the column-stochastic matrix M = α A + (1 − α) E where A is the normalized link matrix and E = (1/n) 1 1^T (teleportation). PageRank is the dominant eigenvector π with M π = π. The power method converges geometrically with rate α (typically α = 0.85).

19.4 Graph Laplacian + spectral clustering

For weighted graph with weight matrix W and degree matrix D = diag(W 1):

    L = D − W                (unnormalized Laplacian)
    L_sym = I − D^{−1/2} W D^{−1/2}    (symmetric normalized)
    L_rw  = I − D^{−1} W              (random-walk Laplacian)

The number of zero eigenvalues equals the number of connected components. The Fiedler vector (eigenvector for the second-smallest eigenvalue) gives a good 2-way graph cut (Fiedler 1973). For k clusters take the k smallest eigenvectors and run k-means on the rows — this is spectral clustering (Shi and Malik 2000; Ng, Jordan, Weiss 2001).

19.5 Quantum mechanics

The time-independent Schrödinger equation is the eigenvalue equation

    Ĥ ψ = E ψ

with Hermitian Hamiltonian Ĥ. Eigenvalues E_k are energy levels; eigenfunctions ψ_k are stationary states. Discretizing Ĥ on a grid or in a basis turns it into a finite-dimensional symmetric eigenproblem of dimension that can reach 10^9+ — driving the development of sparse eigensolvers like Lanczos, FEAST, and SLEPc.

19.6 Structural buckling

Linearized buckling: (K − λ K_g) v = 0, generalized eigenvalue with geometric stiffness K_g. The smallest positive λ is the critical buckling load factor; v is the buckling mode shape. Standard in FEA ([[Engineering/fem-fea]]).

19.7 MEMS resonators and RF filters

Designed by placing natural-frequency modes — the eigenvalue spacing of the mechanical or electromagnetic eigenproblem controls the filter’s pass-band and stop-band.

19.8 Lyapunov stability

Linear stability of ẋ = A x is decided by Re(λ_i(A)). For nonlinear systems, the Jacobian at a fixed point gives local stability via the same test.

19.9 Kalman filtering and covariance evolution

Algebraic Riccati equation solutions and the steady-state error covariance involve eigenvalues of the augmented Hamiltonian matrix. See [[Robotics/bayesian-estimation]].

19.10 Manifold learning

• Laplacian eigenmaps (Belkin and Niyogi 2003) — embed graph by the smallest d nontrivial eigenvectors of L.
• Diffusion maps (Coifman and Lafon 2006) — eigenvectors of a normalized Markov transition matrix on the data.
• Isomap, LLE — eigen-decompositions of derived matrices.

19.11 Transformers and attention

The covariance / Gram matrix of token embeddings, and the spectrum of the attention matrix, are studied to understand learning dynamics, capacity, and the emergence of in-context learning. See [[Compute/transformer-architecture]].

20. Numerical caveats

1. Defective matrices — the matrix of eigenvectors becomes ill-conditioned or singular. Fall back to Schur form; never attempt the Jordan form numerically.
2. Close eigenvalues — individual eigenvectors lose accuracy even when the invariant subspace is fine. Report the subspace (block) rather than individual eigenvectors when eigenvalues are within ε · ‖A‖.
3. Loss of orthogonality in Lanczos/Arnoldi — re-orthogonalize: full (expensive but simple), selective, or partial. Implicit restart in ARPACK bounds the basis size and handles this naturally.
4. Do not compute the characteristic polynomial then root-find — the coefficients of p_A(λ) are catastrophically ill-conditioned (the Wilkinson polynomial example), even though the eigenvalues themselves may be well-conditioned. Always use a similarity-transformation-based algorithm (QR, divide-and-conquer, Lanczos).
5. Iterative methods converge fast at extremes, slow in the middle. For interior eigenvalues use shift-invert ((A − σI)^{-1}), Jacobi-Davidson, or FEAST (contour integration over a specified region of the complex plane).
6. Non-normal matrices — ‖A‖_2 can be far from max|λ|. Pseudospectra reveal transient growth that spectrum alone hides. Critical in fluid stability, neural network dynamics, and resolvent analysis.
7. Scaling and balancing — for non-symmetric problems, scale rows and columns first (LAPACK dgebal) to reduce sensitivity.
8. Convergence tests — relative residual ‖A v − λ v‖ / (‖A‖ · ‖v‖), not absolute, especially when eigenvalues vary by many orders of magnitude.

21. Software stack

Dense

• LAPACK — reference Fortran implementation of every routine mentioned above (dsyev*, dgeev, dgees, dsygv*, dggev, …).
• Intel MKL / OpenBLAS / BLIS — high-performance LAPACK back-ends.
• cuSOLVER — NVIDIA’s GPU LAPACK; covers dense syev and geev.
• MAGMA — heterogeneous CPU + GPU LAPACK, often the fastest large dense solver.
• Eigen3 (C++) — header-only template library; SelfAdjointEigenSolver, EigenSolver, GeneralizedEigenSolver.
• Armadillo (C++) — wraps LAPACK with a MATLAB-like API.

Sparse / large-scale

• ARPACK / ARPACK-NG — implicit restart Arnoldi/Lanczos, the de facto standard.
• SLEPc — PETSc-based, parallel, supports polynomial and nonlinear eigenproblems.
• PRIMME — Jacobi-Davidson, preconditioned, especially good for interior eigenvalues.
• FEAST — contour integration; Intel MKL Extended Eigensolvers.
• Anasazi (Trilinos) — Krylov-Schur, block Davidson, LOBPCG.

High-level languages

• Python — NumPy numpy.linalg.eig, eigh, eigvals, eigvalsh, eigvalsh_tridiagonal. SciPy scipy.sparse.linalg.eigsh (symmetric) and eigs (general), lobpcg. JAX jax.numpy.linalg.eigh (differentiable).
• MATLAB — eig, eigs, polyeig, qz, schur.
• Julia — LinearAlgebra.eigen, eigvals, eigvecs, KrylovKit.jl (eigsolve), Arpack.jl, ArnoldiMethod.jl.
• R — eigen, RSpectra package for sparse.

22. Cross-references

• [[Math/linear-algebra-essentials]] — bases, rank, projections, norms.
• [[Math/numerical-linear-algebra]] — Householder, QR factorization, conditioning, BLAS/LAPACK structure.
• [[Math/svd-pca-spectral]] — singular values, SVD, PCA in detail.
• [[Math/ode-numerical-methods]] — stability of integrators uses eigenvalues of the Jacobian.
• [[Math/convex-optimization]] — semidefinite cones, eigenvalue functionals, log-det.
• [[Engineering/vibration-dynamics]] — modal analysis applications.
• [[Engineering/structural-dynamics]] — full structural eigenanalysis.
• [[Engineering/fem-fea]] — discretization producing the sparse matrices fed to ARPACK/SLEPc.
• [[Robotics/bayesian-estimation]] — Kalman and information matrices.
• [[Compute/transformer-architecture]] — spectral properties of attention and embedding matrices.

23. Citations

• Trefethen, L. N. and Bau, D. (1997). Numerical Linear Algebra. SIAM. Lectures 24-29 on the eigenvalue problem.
• Golub, G. H. and Van Loan, C. F. (2013). Matrix Computations. 4th ed., Johns Hopkins. Chapters 7-8.
• Wilkinson, J. H. (1965). The Algebraic Eigenvalue Problem. Oxford University Press. The classical reference for error analysis.
• Parlett, B. N. (1998). The Symmetric Eigenvalue Problem. SIAM Classics.
• Saad, Y. (2011). Numerical Methods for Large Eigenvalue Problems. 2nd ed., SIAM. The reference for iterative methods.
• Francis, J. G. F. (1961). “The QR transformation, parts I and II.” Computer Journal.
• Kublanovskaya, V. N. (1961). “On some algorithms for the solution of the complete eigenvalue problem.” USSR Computational Mathematics.
• Lanczos, C. (1950). “An iteration method for the solution of the eigenvalue problem of linear differential and integral operators.” J. Res. Natl. Bur. Stand.
• Arnoldi, W. E. (1951). “The principle of minimized iterations in the solution of the matrix eigenvalue problem.” Quart. Appl. Math.
• Cuppen, J. J. M. (1981). “A divide and conquer method for the symmetric tridiagonal eigenproblem.” Numer. Math.
• Dhillon, I. S. and Parlett, B. N. (2004). “Multiple representations to compute orthogonal eigenvectors of symmetric tridiagonal matrices.” Linear Algebra Appl.
• Sleijpen, G. L. G. and van der Vorst, H. A. (1996). “A Jacobi-Davidson iteration method for linear eigenvalue problems.” SIAM J. Matrix Anal. Appl.
• Lehoucq, R. B., Sorensen, D. C., and Yang, C. (1998). ARPACK Users’ Guide: Solution of Large-Scale Eigenvalue Problems with Implicitly Restarted Arnoldi Methods. SIAM.
• Polizzi, E. (2009). “Density-matrix-based algorithm for solving eigenvalue problems.” Phys. Rev. B. (FEAST.)
• Knyazev, A. V. (2001). “Toward the optimal preconditioned eigensolver: Locally Optimal Block Preconditioned Conjugate Gradient method.” SIAM J. Sci. Comput. (LOBPCG.)
• Stathopoulos, A. and McCombs, J. R. (2010). “PRIMME: PReconditioned Iterative MultiMethod Eigensolver.” ACM Trans. Math. Softw.
• Fiedler, M. (1973). “Algebraic connectivity of graphs.” Czechoslovak Math. J.
• Shi, J. and Malik, J. (2000). “Normalized cuts and image segmentation.” IEEE Trans. Pattern Anal. Mach. Intell.
• Ng, A. Y., Jordan, M. I., and Weiss, Y. (2001). “On spectral clustering: Analysis and an algorithm.” NeurIPS.
• Belkin, M. and Niyogi, P. (2003). “Laplacian eigenmaps for dimensionality reduction and data representation.” Neural Comput.
• Coifman, R. R. and Lafon, S. (2006). “Diffusion maps.” Appl. Comput. Harmon. Anal.
• Trefethen, L. N. and Embree, M. (2005). Spectra and Pseudospectra. Princeton University Press.
• Güttel, S., Van Beeumen, R., Meerbergen, K., and Michiels, W. (2014). “NLEIGS: A class of fully rational Krylov methods for nonlinear eigenvalue problems.” SIAM J. Sci. Comput.