Eigenvalue Computation: From Math Derivation to Native C++

The Algebraic Illusion

Every engineer learns eigenvalue computation the exact same way in introductory linear algebra: write down the matrix $A$, construct the characteristic polynomial by computing $\det (A-\lambda I)=0$, factor the polynomial to find its roots ${\lambda }_i$, and solve the resulting homogeneous system $(A-{\lambda }_{i}I)v=0$ for the null space basis vectors.

On paper, this method feels satisfying and complete. In production software, it is an absolute disaster.

The breakdown happens for two distinct reasons: algebraic impossibility and catastrophic numerical instability. First, by the Abel-Ruffini theorem, there is no general algebraic closed-form solution for polynomial roots of degree $n\ge 5$. Second, even for modest $4\times 4$ matrices, computing eigenvalues by finding the roots of an expanded characteristic polynomial is notoriously ill-conditioned. A classic example is Wilkinson's polynomial: shifting a single coefficient in a 20th-degree polynomial by as little as $2^{-23}$ can shift its roots from distinct integers into widely scattered complex numbers.

Numerical linear algebra solves eigensystems by discarding root-finding entirely. Instead of expanding polynomials, production solvers rely on iterative orthogonal similarity transformations that preserve the spectrum while progressively driving off-diagonal elements toward zero.

Understanding how we transition from pencil-and-paper derivations to production-grade LAPACK routines and SIMD-vectorized C++ is essential for building fast, stable systems in robotics, graphics, computer vision, and machine learning.

Spectral Mechanics and Mathematical Invariants

Let $A\in {\mathbb{R}}^{n\times n}$ denote a real square matrix representing a linear transformation on the vector space ${\mathbb{R}}^n$. A non-zero vector $v\in {\mathbb{C}}^n$ is an eigenvector of $A$ associated with eigenvalue $\lambda \in \mathbb{C}$ if and only if:

$Av=\lambda v⟺(A-\lambda I)v=0$

A non-trivial solution $v\ne 0$ exists if and only if the linear operator $(A-\lambda I)$ is singular, which requires its determinant to vanish:

$p(\lambda )=\det (A-\lambda I)=0$

The characteristic polynomial $p(\lambda )$ is an $n$-th degree monic polynomial in $\lambda$. By the Fundamental Theorem of Algebra, $p(\lambda )$ has exactly $n$ complex roots (counting algebraic multiplicity $\mu (\lambda )$, which is the multiplicity of $\lambda$ as a root of $p(\lambda )$).

Two fundamental algebraic invariants tie the operator directly to its spectrum:

1. Trace Invariant: The sum of all diagonal elements equals the sum of all eigenvalues: $\text{Tr}(A)={\sum }_{i=1}^n{A}_{ii}={\sum }_{i=1}^n{\lambda }_i$

2. Determinant Invariant: The determinant equals the product of all eigenvalues: $\det (A)={\prod }_{i=1}^n{\lambda }_i$

Structural Matrix Taxonomy

The matrix structure determines the spectral nature and whether the transformation can be decoupled into orthogonal principal components:

• Real Symmetric Matrices ($A=A^T$): The spectral theorem guarantees that all eigenvalues ${\lambda }_i\in \mathbb{R}$ are strictly real. Eigenvectors corresponding to distinct eigenvalues are mutually orthogonal. Consequently, $A$ admits an orthogonal diagonalization $A=Q\Lambda Q^T$, where $Q$ is an orthogonal matrix ($Q^{T}Q=I$) whose columns are normalized eigenvectors, and $\Lambda =\text{diag}({\lambda }_1,\dots ,{\lambda }_n)$.
• Real Non-Symmetric Matrices ($A\ne A^T$): Eigenvalues can appear as complex conjugate pairs $\lambda =\alpha \pm i\beta$, accompanied by complex conjugate eigenvector pairs $v=u\pm iw$.
• Defective Matrices: A matrix is defective if it possesses an eigenvalue whose geometric multiplicity $\gamma (\lambda )=\dim (\ker (A-\lambda I))$ is strictly less than its algebraic multiplicity $\mu (\lambda )$: $\gamma (\lambda )<\mu (\lambda )$ Defective matrices lack a complete basis of linearly independent eigenvectors, making them non-diagonalizable. Their canonical reduction requires the Jordan normal form with non-trivial Jordan blocks.

Step-by-Step Analytical Derivations

Working through three classical matrix topologies demonstrates how real distinct, defective, and complex spectra behave when evaluated analytically.

1. Distinct Real Spectrum ($2\times 2$)

Consider the non-symmetric matrix:

$A=\left(\begin{array}{cc}4& 1\\ 2& 3\end{array}\right)$

Formulating the characteristic equation:

$\det (A-\lambda I)=\det \left(\begin{array}{cc}4-\lambda & 1\\ 2& 3-\lambda \end{array}\right)=(4-\lambda )(3-\lambda )-(1)(2)={\lambda }^2-7\lambda +10=0$

Factoring gives $(\lambda -5)(\lambda -2)=0$, yielding two distinct real eigenvalues: ${\lambda }_1=5$ and ${\lambda }_2=2$.

To compute the eigenspace for ${\lambda }_1=5$:

$(A-5I)v=\left(\begin{array}{cc}-1& 1\\ 2& -2\end{array}\right)\left(\begin{array}{c}{v}_1\\ v_2\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)⟹-v_1+v_2=0⟹v_1=v_2$

Choosing $v_2=1$ gives the basis vector $v=[1,1{]}^T$. Normalizing under the Euclidean norm $\Vert v{\Vert }_2=\sqrt{1^2+1^2}=\sqrt{2}$:

$u_1=\frac{1}{\sqrt{2}}\left(\begin{array}{c}1\\ 1\end{array}\right)\approx \left(\begin{array}{c}0.7071\\ 0.7071\end{array}\right)$

To compute the eigenspace for ${\lambda }_2=2$:

$(A-2I)v=\left(\begin{array}{cc}2& 1\\ 2& 1\end{array}\right)\left(\begin{array}{c}{v}_1\\ v_2\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)⟹2v_1+v_2=0⟹v_2=-2v_1$

Choosing $v_1=1$ gives $v=[1,-2{]}^T$. Normalizing under $\Vert v{\Vert }_2=\sqrt{1^2+(-2{)}^2}=\sqrt{5}$:

$u_2=\frac{1}{\sqrt{5}}\left(\begin{array}{c}1\\ -2\end{array}\right)\approx \left(\begin{array}{c}0.4472\\ -0.8944\end{array}\right)$

Checking trace and determinant invariants:

• $\text{Tr}(A)=4+3=7={\lambda }_1+{\lambda }_2$
• $\det (A)=(4)(3)-(1)(2)=10={\lambda }_1{\lambda }_2$

2. Defective Operator with Repeated Roots ($3\times 3$)

Consider the matrix:

$A=\left(\begin{array}{ccc}1& 2& 2\\ 0& 2& 1\\ -1& 2& 2\end{array}\right)$

Evaluating the characteristic equation $\det (A-\lambda I)=0$:

$p(\lambda )=-(\lambda -1)(\lambda -2{)}^2=0$

The spectrum consists of ${\lambda }_1=1$ with algebraic multiplicity $\mu (1)=1$, and ${\lambda }_2=2$ with algebraic multiplicity $\mu (2)=2$.

For ${\lambda }_1=1$:

$(A-I)v=\left(\begin{array}{ccc}0& 2& 2\\ 0& 1& 1\\ -1& 2& 1\end{array}\right)\left(\begin{array}{c}{v}_1\\ v_2\\ v_3\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\end{array}\right)$

From row 2, $v_2+v_3=0⟹v_2=-v_3$. Substituting into row 3 gives $-v_1+2(-v_3)+v_3=0⟹v_1=-v_3$. Setting $v_3=1$ yields the eigenvector basis vector $v_{\lambda =1}=[0,-1,1{]}^T$.

For the repeated root ${\lambda }_2=2$:

$(A-2I)v=\left(\begin{array}{ccc}-1& 2& 2\\ 0& 0& 1\\ -1& 2& 0\end{array}\right)\left(\begin{array}{c}{v}_1\\ v_2\\ v_3\end{array}\right)=\left(\begin{array}{c}0\\ 0\\ 0\end{array}\right)$

Row 2 forces $v_3=0$. Substituting $v_3=0$ into row 1 and row 3 leaves $-v_1+2v_2=0⟹v_1=2v_2$.

Setting $v_2=1$ gives $v_{\lambda =2}=[2,1,0{]}^T$.

Because there is only one free variable, the geometric multiplicity is $\gamma (2)=1$. Since $\gamma (2)=1<\mu (2)=2$, matrix $A$ is defective. No invertible matrix $P$ exists such that $P^{-1}AP$ is purely diagonal; $A$ can only be reduced to a Jordan form containing a $2\times 2$ Jordan block.

3. Planar Rotation Matrix with Complex Conjugate Spectrum

Consider the standard $90^{\circ }$ planar rotation matrix:

$R=\left(\begin{array}{cc}0& -1\\ 1& 0\end{array}\right)$

The characteristic equation:

$\det (R-\lambda I)={\lambda }^2+1=0⟹\lambda =\pm i$

For ${\lambda }_1=i$:

$(R-iI)v=\left(\begin{array}{cc}-i& -1\\ 1& -i\end{array}\right)\left(\begin{array}{c}{v}_1\\ v_2\end{array}\right)=\left(\begin{array}{c}0\\ 0\end{array}\right)⟹v_1-i{v}_2=0⟹v_1=i{v}_2$

Setting $v_2=1$ yields the unnormalized vector $[i,1{]}^T$. Normalizing with the Hermitian inner product $\Vert v{\Vert }_2=\sqrt{|i{|}^2+|1{|}^2}=\sqrt{2}$:

$u_1=\frac{1}{\sqrt{2}}\left(\begin{array}{c}i\\ 1\end{array}\right)$

For ${\lambda }_2=-i$, taking the complex conjugate gives:

$u_2=\frac{1}{\sqrt{2}}\left(\begin{array}{c}-i\\ 1\end{array}\right)$

Because rotation preserves vector lengths while continuously changing their directions, no real vector can maintain its direction under $R$. The eigenvectors exist strictly in the complexified vector space ${\mathbb{C}}^2$.

Numerical Algorithms for Dense Eigensystems

Because characteristic polynomial root-finding fails for general matrices, practical algorithms rely on iterative transformations.

Power Iteration: Isolating the Dominant Mode

Power iteration finds the largest eigenvalue in absolute value ($|{\lambda }_1|>|{\lambda }_2|\ge \cdots \ge |{\lambda }_n|$) and its corresponding eigenvector.

Starting with an initial vector $b_0$ with non-zero component along the dominant eigenvector $v_1$, the recurrence relation applies $A$ repeatedly:

$b_{k+1}=\frac{A{b}_k}{\Vert A{b}_k{\Vert }_2}$

As $k\to \infty$, the sequence $b_k$ aligns with the dominant eigenvector $v_1$. The dominant eigenvalue is then extracted using the Rayleigh quotient:

${\lambda }_1=\frac{b_k^{T}A{b}_k}{b_k^T{b}_k}$

The asymptotic error at iteration $k$ scales as:

$\mathcal{O}\left({\left|\frac{{\lambda }_2}{{\lambda }_1}\right|}^k\right)$

When the spectral gap is narrow ($|{\lambda }_2/{\lambda }_1|\approx 1$), Power Iteration converges very slowly. In such scenarios, acceleration methods like Rayleigh Quotient Iteration or spectral shifts $(A-\sigma I{)}^{-1}$ are required.

The QR Algorithm: Extracting the Full Spectrum

The QR algorithm extracts the entire spectrum of a matrix $A$ by executing iterative orthogonal similarity transformations.

At step $k$:

1. Decompose the current matrix $A_k$ into an orthogonal matrix $Q_k$ and an upper triangular matrix $R_k$: $A_k=Q_k{R}_k$
2. Recombine them in reverse order to form $A_{k+1}$: $A_{k+1}=R_k{Q}_k$

Substituting $R_k=Q_k^T{A}_k$ into the second equation:

$A_{k+1}=Q_k^T{A}_k{Q}_k$

Because $A_{k+1}$ is an orthogonal similarity transformation of $A_k$, both matrices have identical eigenvalues across all iterations. As $k\to \infty$, $A_k$ converges to an upper triangular Schur matrix (or quasi-triangular real Schur form with $2\times 2$ blocks for complex conjugate pairs), where the eigenvalues sit directly on the main diagonal.

In production numerical libraries (LAPACK, Eigen3), this raw QR loop is structured into an optimized pipeline:

Code Implementations Across Four Software Paradigms

Let us implement and compare eigensystem routines across four different software layers: symbolic CAS, pure NumPy algorithms, production SciPy/LAPACK, and native SIMD-accelerated C++ with Eigen3.

1. Exact Symbolic Computation via SymPy

SymPy provides exact algebraic derivations without floating-point rounding errors:

import sympy as sp

# Define the 3x3 matrix from our analytical derivation
A = sp.Matrix([
    [ 1, 2, 2],
    [ 0, 2, 1],
    [-1, 2, 2]
])

lam = sp.Symbol('lambda')

# 1. Compute characteristic polynomial symbolically
char_poly = A.charpoly(lam)
print("Characteristic Polynomial:", char_poly.as_expr())

# 2. Extract exact eigenvalues and algebraic multiplicities
evals = A.eigenvals()
print("Eigenvalues (Value: Multiplicity):", evals)

# 3. Extract exact eigenspaces (value, multiplicity, [basis vectors])
evecs = A.eigenvects()
for val, mult, basis in evecs:
    print(f"Eigenvalue {val} (Alg Multiplicity: {mult}, Geo Multiplicity: {len(basis)}):")
    for vec in basis:
        print("  Basis Vector:", vec.T)

Running this script outputs:

Characteristic Polynomial: lambda**3 - 5*lambda**2 + 8*lambda - 4
Eigenvalues (Value: Multiplicity): {1: 1, 2: 2}
Eigenvalue 1 (Alg Multiplicity: 1, Geo Multiplicity: 1):
  Basis Vector: Matrix([[0, -1, 1]])
Eigenvalue 2 (Alg Multiplicity: 2, Geo Multiplicity: 1):
  Basis Vector: Matrix([[2, 1, 0]])

The output confirms our manual derivation: $\lambda =2$ has algebraic multiplicity 2 but only a 1-dimensional eigenspace, proving defectivity.

2. Custom Iterative Solvers in NumPy

The following script implements basic Power Iteration and the QR Algorithm directly from scratch:

import numpy as np

def power_iteration(A: np.ndarray, num_iterations: int = 100) -> tuple[float, np.ndarray]:
    np.random.seed(42)
    b_k = np.random.rand(A.shape[1])
    b_k = b_k / np.linalg.norm(b_k)
    
    for _ in range(num_iterations):
        b_k1 = np.dot(A, b_k)
        norm = np.linalg.norm(b_k1)
        if norm == 0:
            break
        b_k = b_k1 / norm
        
    dominant_eigenvalue = float(np.dot(b_k.T, np.dot(A, b_k)) / np.dot(b_k.T, b_k))
    return dominant_eigenvalue, b_k

def qr_algorithm(A: np.ndarray, num_iterations: int = 50) -> tuple[np.ndarray, np.ndarray]:
    A_k = np.copy(A).astype(float)
    n = A.shape[0]
    Q_accum = np.eye(n)
    
    for _ in range(num_iterations):
        Q, R = np.linalg.qr(A_k)
        A_k = R @ Q
        Q_accum = Q_accum @ Q
        
    eigenvalues = np.diag(A_k)
    return eigenvalues, Q_accum

# Real symmetric test matrix
A_sym = np.array([
    [ 6.0, -2.0,  1.0],
    [-2.0,  3.0, -2.0],
    [ 1.0, -2.0,  6.0]
], dtype=float)

dom_val, dom_vec = power_iteration(A_sym, num_iterations=100)
qr_vals, qr_evecs = qr_algorithm(A_sym, num_iterations=50)

print(f"Power Iteration Dominant Eigenvalue: {dom_val:.6f}")
print(f"Power Iteration Dominant Eigenvector: {dom_vec}")
print(f"QR Algorithm Extracted Spectrum: {np.sort(qr_vals)[::-1]}")

Running this code produces:

Power Iteration Dominant Eigenvalue: 8.000000
Power Iteration Dominant Eigenvector: [ 0.70710678 -0.         -0.70710678]
QR Algorithm Extracted Spectrum: [8. 5. 2.]

Power Iteration isolates the dominant mode $\lambda =8$, while the QR algorithm extracts the entire spectrum $\{8,5,2\}$.

3. Production LAPACK Drivers via NumPy and SciPy

High-level Python libraries wrap Fortran LAPACK routines (dgeev for general non-symmetric matrices, dsyev for symmetric matrices, and dggev for generalized eigensystems):

import numpy as np
import scipy.linalg as la

# 1. Real Symmetric Matrix -> dsyev (faster, guaranteed orthogonal real vectors)
A_sym = np.array([[6.0, -2.0, 1.0], [-2.0, 3.0, -2.0], [1.0, -2.0, 6.0]])
evals_sym, evecs_sym = la.eigh(A_sym)
print("Symmetric Eigenvalues (la.eigh -> dsyev):", evals_sym)

# 2. General Non-Symmetric Matrix -> dgeev
A_gen = np.array([[4.0, 1.0], [2.0, 3.0]])
evals_gen, evecs_gen = la.eig(A_gen)
print("General Eigenvalues (la.eig -> dgeev):", evals_gen.real)

# 3. Planar Rotation Matrix (Complex Spectrum)
D_rot = np.array([[0.0, -1.0], [1.0, 0.0]])
evals_rot, evecs_rot = la.eig(D_rot)
print("Rotation Matrix Complex Eigenvalues:", evals_rot)

# 4. Generalized Eigensystem: A * v = lambda * B * v
B_metric = np.array([[2.0, 1.0], [1.0, 2.0]])
gen_evals, gen_evecs = la.eig(A_gen, B_metric)
print("Generalized Eigenvalues:", gen_evals.real)

Running this script yields:

Symmetric Eigenvalues (la.eigh -> dsyev): [2. 5. 8.]
General Eigenvalues (la.eig -> dgeev): [5. 2.]
Rotation Matrix Complex Eigenvalues: [0.+1.j 0.-1.j]
Generalized Eigenvalues: [1. 3.]

For defective matrices, numerical eigensolvers return duplicate eigenvectors that are linearly dependent within machine precision tolerances ($\approx 10^{-16}$). When complex conjugate roots appear, NumPy automatically promotes real floating-point buffers to complex128.

4. Native C++ with Eigen3 and SIMD Vectorization

In latency-critical applications (such as robotics state estimation, physics engines, and high-frequency trading), Python's FFI overhead is prohibitive. Compiled C++ using Eigen3 executes matrix decomposition directly on CPU registers:

#include <iostream>
#include <complex>
#include <Eigen/Dense>
#include <Eigen/Eigenvalues>

int main() {
    // Stack-allocated 3x3 double-precision matrix
    Eigen::Matrix3d A;
    A <<  6.0, -2.0,  1.0,
         -2.0,  3.0, -2.0,
          1.0, -2.0,  6.0;

    // SelfAdjointEigenSolver exploits symmetry for maximum performance
    Eigen::SelfAdjointEigenSolver<Eigen::Matrix3d> solver(A);

    if (solver.info() != Eigen::Success) {
        std::cerr << "Eigensolver computation failed to converge!" << std::endl;
        return 1;
    }

    Eigen::Vector3d eigenvalues = solver.eigenvalues();
    Eigen::Matrix3d eigenvectors = solver.eigenvectors();

    std::cout << "Computed Real Eigenvalues:\n" << eigenvalues << "\n\n";
    std::cout << "Computed Orthogonal Eigenvectors (Columns):\n" << eigenvectors << "\n\n";

    // Verify residual norm: ||A * v - lambda * v||_2 for dominant mode
    double lambda_dom = eigenvalues[2]; // Sorted in increasing order
    Eigen::Vector3d v_dom = eigenvectors.col(2);

    Eigen::Vector3d residual = A * v_dom - lambda_dom * v_dom;
    std::cout << "Residual Norm ||A*v - lambda*v||_2: " << residual.norm() << std::endl;

    return 0;
}

Compiling with aggressive compiler flags:

g++ -O3 -march=native -DNDEBUG -I/usr/include/eigen3 eigen_demo.cpp -o eigen_demo
./eigen_demo

Output:

Computed Real Eigenvalues:
2
5
8

Computed Orthogonal Eigenvectors (Columns):
 0.408248 -0.57735   0.707107
 0.816497  0.57735  -0.000000
-0.408248  0.57735   0.707107

Residual Norm ||A*v - lambda*v||_2: 0

Architectural Comparison and Hardware Performance

The choice of execution paradigm introduces fundamental trade-offs across mathematical precision, memory allocation, and instruction-level parallelism.

The Foreign Function Interface (FFI) Tax

When computing eigensystems for massive matrices ($N>2000$), NumPy and SciPy perform almost identically to native C++ because the execution time is entirely dominated by compute-intensive Level-3 BLAS operations inside the LAPACK routine.

However, for low-dimensional systems ($N\in [2,6]$) evaluated millions of times inside tight real-time loops (e.g., orientation kinematics, extended Kalman filtering, surface normal estimation):

• Python FFI Overhead: Every call to np.linalg.eig or scipy.linalg.eigh pays for Python function call overhead, argument type checking, C-FFI pointer marshaling, temporary buffer allocations, and GIL acquisition. This overhead often requires 1.5 to 5 microseconds per call.
• Native Eigen3 Stack Execution: A fixed-size matrix (Eigen::Matrix3d) allocates its 9 double values directly on the stack ($72$ bytes). The compiler inlines the solver, unrolls all transformation loops, and generates direct vectorized SIMD instructions (vfmadd213pd), executing the entire eigensystem extraction in under 40 nanoseconds.

When to Use Which Tool

Selecting the right eigensystem strategy depends directly on your system's scale, sparsity, and latency budget:

• Use SymPy when deriving theoretical proofs, checking algebraic vs. geometric multiplicities, or verifying exact rational eigenspaces of small symbolic operators.
• Use Custom Iterative Solvers when studying convergence dynamics, spectral shifts, or testing specialized numerical preconditions.
• Use NumPy / SciPy (scipy.linalg.eigh / eig) for data pipelines, offline ML research, PCA, clustering, and general scientific computing where development velocity outweighs sub-microsecond latency.
• Use Native C++ (Eigen3) for real-time control, robotics kinematics, physical simulations, and embedded hardware where deterministic execution and stack-allocated SIMD performance are non-negotiable.

When matrices become sparse and high-dimensional ($N>10^6$, such as web graphs or graph neural networks), dense solvers become completely intractable. In those regimes, the focus shifts to Krylov subspace methods: Arnoldi iteration for non-symmetric matrices and Lanczos factorization for symmetric systems, computing only the extremal eigenvalues without ever materializing the full dense operator.