Linear Algebra II: Eigenvalues & SVD

The determinant of a square matrix A, written det(A) or |A|, is a single number that encodes a fundamental geometric fact: how does A scale volumes? If you take the unit cube (or unit square in 2D) and apply A, the resulting parallelepiped has volume |det(A)| times the original. If det(A) = 0, the transformation collapses space into a lower dimension — it is singular and non-invertible. If det(A) < 0, the transformation includes a reflection (orientation is flipped, like a mirror image).

Computing determinants:

• For 2×2: det([[a,b],[c,d]]) = ad − bc. Geometrically, this is the signed area of the parallelogram spanned by the column vectors.
• For 3×3 and beyond: Laplace expansion along any row or column. For row i: det(A) = Σⱼ (−1)^(i+j) · aᵢⱼ · det(Mᵢⱼ), where Mᵢⱼ is the (n−1)×(n−1) minor obtained by deleting row i and column j.
• In practice, LU decomposition computes det(A) in O(n³) as the product of diagonal entries of U (times a sign for row swaps).

Key algebraic properties:

• det(AB) = det(A) · det(B): the volume-scaling of two transformations in sequence multiplies.
• det(Aᵀ) = det(A): transpose doesn't change the determinant.
• det(A⁻¹) = 1/det(A): the inverse undoes the scaling.
• det(cA) = cⁿ det(A) for n×n matrix: scaling every entry by c scales each of n dimensions.

Why determinants matter for eigenvalues: The key bridge is that A is invertible if and only if det(A) ≠ 0. We will use this when defining eigenvalues via det(A − λI) = 0.

The Leibniz formula gives a direct but computationally expensive expression: det(A) = Σ_{σ ∈ Sₙ} sgn(σ) · Π_i a_{i,σ(i)}, summing over all n! permutations of {1,...,n}.

import numpy as np

# 2x2 determinant — geometric interpretation
A = np.array([[3., 1.],
              [1., 2.]])
print(f"det(A) = {np.linalg.det(A):.4f}")  # 5.0
# Area of parallelogram spanned by columns [3,1] and [1,2]

# Volume scaling: det(AB) = det(A)*det(B)
B = np.array([[2., 0.], [0., 3.]])
print(np.linalg.det(A @ B))            # 30.0
print(np.linalg.det(A) * np.linalg.det(B))  # 30.0

# Singular matrix: det = 0
C = np.array([[1., 2.], [2., 4.]])
print(f"det(C) = {np.linalg.det(C):.4f}")  # ~0 (columns are dependent)

A nonzero vector v is an eigenvector of matrix A with eigenvalue λ if:

Av = λv

This says: when A acts on v, it does not rotate v at all — it only scales it by the factor λ. If λ = 2, A doubles the length. If λ = −1, A reflects it. If λ = 0, A collapses it to zero. These special directions are the 'bones' of the transformation.

Finding eigenvalues — the characteristic polynomial: Av = λv ⟺ (A − λI)v = 0. For a nonzero solution to exist, (A − λI) must be singular, i.e., det(A − λI) = 0. This gives the characteristic polynomial p(λ) = det(A − λI), a degree-n polynomial in λ. Its roots are the eigenvalues.

For a 2×2 matrix A = [[a,b],[c,d]]: p(λ) = λ² − (a+d)λ + (ad−bc) = λ² − tr(A)·λ + det(A) = 0

This gives two eigenvalues (possibly complex, possibly repeated). The sum of eigenvalues equals tr(A) (the trace, sum of diagonal entries), and their product equals det(A). These are beautiful invariants.

Finding eigenvectors: For each eigenvalue λᵢ, solve (A − λᵢI)v = 0. The solution space is the eigenspace E(λᵢ) = null(A − λᵢI). The eigenspace always has dimension ≥ 1 (since λᵢ is a root, the matrix is singular).

Geometric vs. algebraic multiplicity: If λ appears as a root of p(λ) with multiplicity k, k is the algebraic multiplicity. The dimension of E(λ) is the geometric multiplicity, which is always ≤ algebraic multiplicity. When they differ, the matrix cannot be fully diagonalized.

Analogy: Think of a drum. When you strike a drum, it vibrates at specific resonant frequencies — those frequencies are its eigenvalues, and the corresponding vibration shapes are its eigenvectors. Every linear system has its own set of 'natural frequencies'.

Important classes:

• If A is triangular, the eigenvalues are exactly the diagonal entries.
• Similar matrices (B = P⁻¹AP) have the same eigenvalues.
• Eigenvalues of Aᵏ are λᵢᵏ. Eigenvalues of A⁻¹ are 1/λᵢ.

import numpy as np

A = np.array([[4., 2.],
              [1., 3.]])

# Characteristic polynomial: det(A - lambda*I) = 0
# = (4-lambda)(3-lambda) - 2 = lambda^2 - 7*lambda + 10 = 0
# => lambda = 5 or lambda = 2

eigenvalues, eigenvectors = np.linalg.eig(A)
print("Eigenvalues:", eigenvalues)    # [5. 2.]
print("Eigenvectors (columns):\n", eigenvectors)

# Verify: A @ v = lambda * v
for i in range(2):
    lam = eigenvalues[i]
    v   = eigenvectors[:, i]
    print(f"A*v{i}: {A @ v}, lambda*v{i}: {lam * v}")  # Should match

# Trace = sum of eigenvalues, det = product of eigenvalues
print("Trace:", np.trace(A), "Sum λ:", np.sum(eigenvalues))  # 7, 7
print("Det:",  np.linalg.det(A), "Prod λ:", np.prod(eigenvalues))  # 10, 10

If an n×n matrix A has n linearly independent eigenvectors v₁, v₂, ..., vₙ with corresponding eigenvalues λ₁, λ₂, ..., λₙ, then A can be diagonalized:

A = P Λ P⁻¹

where P = [v₁ | v₂ | ... | vₙ] (eigenvectors as columns) and Λ = diag(λ₁, ..., λₙ).

This is a change-of-basis insight: in the eigenbasis, A acts by simply scaling each coordinate. P converts from the eigenbasis to the standard basis; P⁻¹ converts back.

When is A diagonalizable?

• If A has n distinct eigenvalues, it is always diagonalizable (distinct eigenvalues guarantee n linearly independent eigenvectors).
• A may be diagonalizable even with repeated eigenvalues, as long as the geometric multiplicity equals the algebraic multiplicity for every eigenvalue.
• If any eigenspace is 'deficient' (geo < alg multiplicity), A is not diagonalizable. The Jordan normal form is the closest you can get.

Why diagonalization is powerful — computing Aᵏ: A = PΛP⁻¹ ⟹ Aᵏ = PΛᵏP⁻¹, and Λᵏ = diag(λ₁ᵏ, ..., λₙᵏ). Matrix exponentiation that normally costs O(kn³) now costs O(n³) (just one diagonalization) plus O(n) scalar exponentiations. This is how you compute Fibonacci numbers in O(log n) or solve linear recurrences.

Solving linear ODEs: The system dx/dt = Ax has solution x(t) = e^{At} x(0). With diagonalization, e^{At} = P e^{Λt} P⁻¹, where e^{Λt} = diag(e^{λ₁t}, ..., e^{λₙt}). The solution is a sum of exponential modes e^{λᵢt} vᵢ — the eigenvalues determine stability (negative real part → stable decay).

Intuition: Diagonalization is like switching languages. The standard basis is like English — general-purpose but not perfectly suited to A. The eigenbasis is like A's native language, in which A speaks with perfect clarity: 'I scale direction 1 by λ₁, direction 2 by λ₂, ...'

import numpy as np

A = np.array([[4., 2.],
              [1., 3.]])

L, P = np.linalg.eig(A)   # L = eigenvalues, P = eigenvectors as columns
Lambda = np.diag(L)
P_inv  = np.linalg.inv(P)

# Verify A = P @ Lambda @ P_inv
print(np.allclose(A, P @ Lambda @ P_inv))  # True

# Compute A^10 efficiently via diagonalization
def matrix_power_via_diag(P, L, P_inv, k):
    return P @ np.diag(L**k) @ P_inv

print(np.allclose(np.linalg.matrix_power(A, 10),
                  matrix_power_via_diag(P, L, P_inv, 10)))  # True

# Fibonacci via matrix exponentiation
# F_{n+1}   =  [[1,1],[1,0]]^n  applied to [F_1, F_0]
Fib = np.array([[1.,1.],[1.,0.]])
Lf, Pf = np.linalg.eig(Fib)
def fib(n):
    Fn = Pf @ np.diag(Lf**n) @ np.linalg.inv(Pf)
    return round(Fn[1, 0])   # F_n element
print([fib(i) for i in range(1, 11)])  # [1,1,2,3,5,8,13,21,34,55]

The Spectral Theorem is one of the crown jewels of linear algebra. It says:

Every real symmetric matrix A (Aᵀ = A) has:

1. All real eigenvalues (even if the entries look like they might yield complex roots).
2. Eigenvectors corresponding to distinct eigenvalues are orthogonal.
3. There exists an orthonormal eigenbasis: A = Q Λ Qᵀ, where Q is orthogonal (Qᵀ = Q⁻¹) and Λ is diagonal with real entries.

This decomposition A = QΛQᵀ is called the eigendecomposition (or spectral decomposition) of a symmetric matrix.

Proof sketch of real eigenvalues: Suppose Av = λv with v ≠ 0. Take the inner product: v* A v = λ vv (where * is conjugate transpose). Since A is symmetric and real, v A v = (Av)* v = (λv)* v = λ̄ v*v. So λ = λ̄, meaning λ is real.

Proof sketch of orthogonality: If Av₁ = λ₁v₁ and Av₂ = λ₂v₂ with λ₁ ≠ λ₂, then λ₁(v₁·v₂) = (Av₁)·v₂ = v₁·(Av₂) = λ₂(v₁·v₂) (using symmetry). So (λ₁ − λ₂)(v₁·v₂) = 0; since λ₁ ≠ λ₂, v₁·v₂ = 0.

Outer product form: We can write A = λ₁q₁q₁ᵀ + λ₂q₂q₂ᵀ + ... + λₙqₙqₙᵀ, a sum of rank-1 projection matrices scaled by eigenvalues. Each qᵢqᵢᵀ projects onto the i-th eigendirection. This is the spectral decomposition — 'spectrum' refers to the set of eigenvalues.

Why this matters profoundly:

• Covariance matrices (Σ = XᵀX / n) are positive semidefinite symmetric: all eigenvalues ≥ 0. PCA finds their eigenvectors.
• Graph Laplacian matrices are symmetric: spectral clustering uses their eigenvalues to detect communities.
• Quadratic forms xᵀAx with symmetric A are fully characterized by eigenvalues: positive definite ⟺ all λᵢ > 0.

Positive (semi)definiteness: A symmetric matrix is positive definite (PD) if all eigenvalues > 0, positive semidefinite (PSD) if all ≥ 0. Equivalent characterization: xᵀAx > 0 for all nonzero x.

import numpy as np

# Construct a random symmetric matrix
rng = np.random.default_rng(42)
B   = rng.standard_normal((4, 4))
A   = B @ B.T          # Always PSD: xTAx = xT(BBT)x = ||BTx||^2 >= 0

L, Q = np.linalg.eigh(A)   # eigh: exploit symmetry, guarantees real eigenvalues
print("Eigenvalues (all >= 0 since PSD):", np.round(L, 4))
print("Q orthogonal? QTQ = I:", np.allclose(Q.T @ Q, np.eye(4)))

# Spectral decomposition: A = Q Lambda QT
A_reconstructed = Q @ np.diag(L) @ Q.T
print("Reconstruction error:", np.max(np.abs(A - A_reconstructed)))  # ~1e-14

# Outer product form: sum of rank-1 pieces
A_sum = sum(L[i] * np.outer(Q[:,i], Q[:,i]) for i in range(4))
print("Outer product form error:", np.max(np.abs(A - A_sum)))  # ~1e-14

# Use np.linalg.eig for non-symmetric and eigh for symmetric (faster, stabler)
A_nonsymm = rng.standard_normal((4, 4))
L_ns, V_ns = np.linalg.eig(A_nonsymm)
print("Non-symmetric eigenvalues may be complex:", L_ns.dtype)  # complex128

The SVD is the generalization of eigendecomposition that works for any m×n matrix — not just square, not just symmetric, not even just invertible. It is arguably the single most important matrix factorization in applied mathematics.

The theorem: Every m×n matrix A can be written as:

A = U Σ Vᵀ

where:

• U is m×m orthogonal (columns are left singular vectors, orthonormal in Rᵐ)
• Σ is m×n diagonal with nonneg entries σ₁ ≥ σ₂ ≥ ... ≥ σₚ ≥ 0 (p = min(m,n)) — the singular values
• V is n×n orthogonal (columns are right singular vectors, orthonormal in Rⁿ)

Geometric picture: A = UΣVᵀ says: any linear map factors into three steps:

1. Vᵀ: rotate/reflect in the input space Rⁿ (change basis to right singular vectors)
2. Σ: scale each dimension independently (stretch by σ₁, σ₂, ...)
3. U: rotate/reflect in the output space Rᵐ (change basis to left singular vectors)

A ball in Rⁿ maps under A to an ellipsoid in Rᵐ. The singular values are the semi-axis lengths; the right singular vectors are the input directions that become those axes; the left singular vectors are those axes in the output space.

Connection to eigendecomposition:

• AᵀA = V Σᵀ Uᵀ U Σ Vᵀ = V (ΣᵀΣ) Vᵀ: right singular vectors are eigenvectors of AᵀA, σᵢ² are eigenvalues.
• AAᵀ = U (ΣΣᵀ) Uᵀ: left singular vectors are eigenvectors of AAᵀ.

Both AᵀA and AAᵀ are symmetric PSD — so we can use eigh for numerical stability.

Rank via SVD: The rank of A equals the number of nonzero singular values. The column space of A is spanned by the first r left singular vectors; the row space by the first r right singular vectors; the null space by the last n−r right singular vectors.

The operator/spectral norms:

• ||A||₂ = σ₁ (largest singular value) — the maximum stretch factor
• ||A||_F = sqrt(Σσᵢ²) (Frobenius norm)
• ||A⁻¹||₂ = 1/σₙ → condition number κ(A) = σ₁/σₙ

import numpy as np

# SVD of a 4x3 matrix
A = np.array([[1., 2., 3.],
              [4., 5., 6.],
              [7., 8., 9.],
              [1., 0., 1.]])

U, s, Vt = np.linalg.svd(A, full_matrices=True)
print("U shape:", U.shape)    # (4, 4)
print("s:",       s)           # [1.53e+01, 3.50e+00, ...]
print("Vt shape:", Vt.shape)  # (3, 3)

# Reconstruct A from full SVD
Sigma = np.zeros_like(A)          # 4x3
np.fill_diagonal(Sigma, s)
A_reconstructed = U @ Sigma @ Vt
print("Reconstruction error:", np.max(np.abs(A - A_reconstructed)))  # ~1e-14

# Rank: number of significant singular values
rank = np.sum(s > 1e-10)
print("Rank of A:", rank)    # 2 (rows 1,2,3 are nearly linearly dependent)

# Operator (spectral) norm = largest singular value
print("||A||_2 =", s[0])                  # 15.32...
print("||A||_F =", np.linalg.norm(A, 'fro'))  # sqrt(sum of squared entries)

Perhaps the most practically powerful consequence of SVD is the Eckart-Young-Mirsky Theorem (1936):

The best rank-k approximation to A (in both the Frobenius norm and the spectral norm) is:

Aₖ = σ₁u₁v₁ᵀ + σ₂u₂v₂ᵀ + ... + σₖuₖvₖᵀ

— i.e., keep only the k largest singular values and their corresponding singular vectors. No other rank-k matrix is closer to A.

Proof intuition: We are writing A as a sum of rank-1 outer products σᵢuᵢvᵢᵀ, each contributing σᵢ² to the total squared Frobenius norm. Truncating to the k largest is optimal because no reweighting or rotation of rank-1 terms can do better.

Information compression: The approximation Aₖ requires storing k(m + n + 1) numbers instead of mn. For an image matrix with m = n = 1000 and k = 50, that is 50 × 2001 ≈ 100,000 numbers instead of 1,000,000 — a 10× compression.

Applications:

• Image compression: A grayscale image is an m×n matrix. Aₖ for k much smaller than min(m,n) retains visible structure with high compression ratio.
• PCA: Principal Component Analysis of a data matrix X (n samples × d features, mean-centered) finds the top-k right singular vectors of X. These are the 'principal components' — directions of maximum variance.
• Latent Semantic Analysis (LSA): A term-document matrix A; Aₖ captures latent semantic topics.
• LoRA (Low-Rank Adaptation): In LLM fine-tuning, the weight update ΔW is constrained to be low-rank: ΔW = UV (U ∈ Rᵐˣᵏ, V ∈ Rᵏˣⁿ, k << min(m,n)), massively reducing trainable parameters.
• Recommendation systems: The user-item rating matrix is low-rank; SVD factors find latent user preferences and item properties.

The reconstruction error is ||A − Aₖ||_F = sqrt(σₖ₊₁² + ... + σₚ²): the singular values tail tells you exactly how much information is lost.

import numpy as np
import matplotlib
matplotlib.use('Agg')           # non-interactive backend
import matplotlib.pyplot as plt

rng = np.random.default_rng(0)
# Simulate a 100x80 'image' with true rank 5 + noise
true_rank = 5
U0 = rng.standard_normal((100, true_rank))
V0 = rng.standard_normal((80,  true_rank))
A  = U0 @ V0.T + 0.5 * rng.standard_normal((100, 80))

U, s, Vt = np.linalg.svd(A, full_matrices=False)  # economy SVD

def low_rank_approx(U, s, Vt, k):
    """Best rank-k approximation via Eckart-Young."""
    return (U[:, :k] * s[:k]) @ Vt[:k, :]

# Reconstruction errors
for k in [1, 3, 5, 10, 20]:
    Ak  = low_rank_approx(U, s, Vt, k)
    err = np.linalg.norm(A - Ak, 'fro')
    compression = k * (100 + 80 + 1) / (100 * 80)
    print(f"k={k:2d}  ||error||_F={err:.3f}  compression={compression:.2%}")

# Variance explained (like PCA scree plot)
var_explained = np.cumsum(s**2) / np.sum(s**2)
print("Variance explained by top 5 singular values:",
      f"{var_explained[4]:.1%}")

Principal Component Analysis (PCA) is the most widely used dimensionality reduction technique in data science and machine learning. It is, at its core, an eigendecomposition problem.

Setup: Given an n×d data matrix X (n samples, d features), first mean-center: X̃ = X − 1·μᵀ (subtract column means). The empirical covariance matrix is C = X̃ᵀX̃ / (n−1), which is d×d symmetric PSD.

PCA via eigendecomposition: Find the orthonormal eigenvectors of C (i.e., do spectral decomposition C = QΛQᵀ). The eigenvectors sorted by decreasing eigenvalue are the principal components — the directions of maximum variance. The eigenvalue λᵢ equals the variance of the data projected onto the i-th principal component.

PCA via SVD (numerically preferred): The singular values of X̃ / sqrt(n−1) are the standard deviations along principal components; the right singular vectors of X̃ are the principal components. SVD is preferred because it avoids explicitly forming the potentially ill-conditioned matrix X̃ᵀX̃.

Choosing k: Plot the 'scree plot' (eigenvalues vs. index). Look for an 'elbow' or choose k such that the top-k eigenvalues explain ≥ 90% of total variance (Σᵢ≤ₖ λᵢ / Σᵢ λᵢ).

Projection and reconstruction:

• Project: Z = X̃ Q_k (n×k low-dimensional representation)
• Reconstruct: X̃_reconstructed = Z Q_kᵀ (Eckart-Young guarantees this is the best k-dim reconstruction)

Geometric interpretation: PCA finds the k-dimensional hyperplane (affine subspace) that best fits the data cloud in the least-squares sense. The principal components are the axes of an ellipsoid fitted to the data.

Limitations to be aware of: PCA is a linear method. It finds the best linear subspace. Non-linear structure (Swiss roll, circular manifolds) requires non-linear methods like t-SNE, UMAP, or kernel PCA. Also, PCA is not invariant to feature scaling — always standardize features if they have different units.

import numpy as np

rng = np.random.default_rng(7)
# 200 samples, 5 features, true structure is 2D
W  = rng.standard_normal((5, 2))   # 5x2 loading matrix
Z0 = rng.standard_normal((200, 2)) # 200x2 latent coords
X  = Z0 @ W.T + 0.1 * rng.standard_normal((200, 5))

# Step 1: mean-center
X_c = X - X.mean(axis=0)

# Step 2: SVD (numerically preferred over eig on covariance)
U, s, Vt = np.linalg.svd(X_c, full_matrices=False)
n = X_c.shape[0]
principal_components = Vt          # rows are PCs (right singular vectors)
standard_devs        = s / np.sqrt(n - 1)
variance_ratio       = (s**2) / np.sum(s**2)

print("Variance explained per PC:",
      np.round(variance_ratio[:3], 4))   # ~[0.88, 0.08, ...]  (2D structure)

# Project to top-2 PCs
k = 2
Z_pca = X_c @ Vt[:k].T            # 200 x 2
print("Reduced representation shape:", Z_pca.shape)  # (200, 2)

# Reconstruct and compute reconstruction error
X_recon = Z_pca @ Vt[:k]
loss = np.mean((X_c - X_recon)**2)
print(f"Reconstruction MSE with k=2: {loss:.5f}")

A quadratic form is a function Q: Rⁿ → R of the form Q(x) = xᵀAx, where A is symmetric. Expanding: Q(x) = Σᵢⱼ aᵢⱼ xᵢxⱼ. Examples: the squared length ||x||² = xᵀIx, and the Mahalanobis distance xᵀΣ⁻¹x in statistics.

Connection to eigenvalues: Since A = QΛQᵀ (spectral theorem), let y = Qᵀx (orthogonal change of basis). Then Q(x) = xᵀQΛQᵀx = yᵀΛy = Σᵢ λᵢyᵢ². The quadratic form is a sum of scaled squares. This is the key insight: in the eigenbasis, a quadratic form is perfectly transparent.

Definiteness classification:

• Positive definite (PD): xᵀAx > 0 for all x ≠ 0 ⟺ all eigenvalues > 0
• Positive semidefinite (PSD): xᵀAx ≥ 0 ⟺ all eigenvalues ≥ 0
• Negative definite (ND): xᵀAx < 0 ⟺ all eigenvalues < 0
• Indefinite: some positive, some negative eigenvalues

Why this matters in optimization: The second-order Taylor expansion of f: Rⁿ → R near a critical point x* is f(x* + δ) ≈ f(x*) + δᵀ∇²f(x*)δ (Hessian). If ∇²f is PD, x* is a local minimum. If ND, a local maximum. If indefinite, a saddle point. Eigenvalues of the Hessian tell you the curvature in each eigendirection — which is why Newton's method uses the Hessian for second-order optimization.

Cholesky decomposition: A matrix is PD if and only if it has a Cholesky factorization A = LLᵀ (L is lower triangular with positive diagonal). This is the standard computational test for PD and is twice as fast as LU for symmetric matrices.

Level sets: The level set {x : xᵀAx = c} for PD A is an ellipsoid whose axes are the eigenvectors of A with lengths proportional to 1/sqrt(λᵢ). Larger eigenvalues → more curvature → shorter axis (the ellipsoid is 'squeezed' in that direction).

import numpy as np

# Check positive definiteness via eigenvalues and Cholesky
def is_positive_definite(A):
    try:
        np.linalg.cholesky(A)
        return True
    except np.linalg.LinAlgError:
        return False

# PD matrix (covariance)
A_pd  = np.array([[4., 2.], [2., 3.]])
L_pd  = np.linalg.eigvalsh(A_pd)
print("Eigenvalues:", L_pd, "=> PD?", np.all(L_pd > 0))  # [1.59, 5.41] True
print("Cholesky test:", is_positive_definite(A_pd))       # True

# Indefinite matrix (saddle point in optimization)
A_indef = np.array([[1., 0.], [0., -1.]])
L_indef = np.linalg.eigvalsh(A_indef)
print("Eigenvalues:", L_indef, "=> indefinite")  # [-1, 1]

# Mahalanobis distance with covariance Sigma
Sigma = np.array([[4., 1.], [1., 2.]])
Sigma_inv = np.linalg.inv(Sigma)
x = np.array([1., 1.])
mahal_sq = x @ Sigma_inv @ x
print(f"Mahalanobis^2 distance: {mahal_sq:.4f}")  # Accounts for covariance

The SVD gives the definitive answer to the question: what is the 'best' solution to Ax = b when A is not square or not invertible?

The Moore-Penrose pseudoinverse A⁺ (also written A†) is defined via SVD: if A = UΣVᵀ, then:

A⁺ = V Σ⁺ Uᵀ

where Σ⁺ replaces each nonzero σᵢ with 1/σᵢ (and leaves zeros as zeros). The pseudoinverse is the unique matrix satisfying the four Moore-Penrose conditions: AA⁺A = A, A⁺AA⁺ = A⁺, (AA⁺)ᵀ = AA⁺, (A⁺A)ᵀ = A⁺A.

Least squares via A⁺: The minimum-norm least squares solution to Ax = b is x* = A⁺b. It minimizes ||Ax − b||₂ (the residual) and among all solutions achieving the minimum residual, it has the smallest ||x||₂. This is what numpy.linalg.lstsq computes internally.

Three cases:

1. m > n, full column rank (overdetermined, consistent → unique solution): A⁺ = (AᵀA)⁻¹Aᵀ (the classic normal equations least squares)
2. m < n, full row rank (underdetermined → infinitely many solutions): A⁺ = Aᵀ(AAᵀ)⁻¹ (minimum-norm solution)
3. Rank deficient: Only SVD handles this correctly. Numerical implementations threshold small singular values (below max(σ) × ε_machine × max(m,n)) as zero.

Numerical stability: Direct inversion via normal equations AᵀAx = Aᵀb squares the condition number (κ(AᵀA) = κ(A)²). For ill-conditioned A, use QR decomposition or SVD directly. This is why modern least-squares solvers do not form AᵀA.

Geometric projection: A⁺b = V(1/σ)Uᵀb. The matrix UᵀU_r = AA⁺ is the orthogonal projector onto the column space of A. A⁺A is the orthogonal projector onto the row space. These are the fundamental projection matrices of linear algebra.

import numpy as np

# Overdetermined system: more equations than unknowns
A = np.array([[1., 1.],
              [1., 2.],
              [1., 3.],
              [1., 4.]])
b = np.array([2., 3., 3., 5.])

# Solve via SVD pseudoinverse
U, s, Vt = np.linalg.svd(A, full_matrices=False)
# A+ = V Sigma+ UT
Sigma_plus = np.diag(1.0 / s)
x_star = Vt.T @ Sigma_plus @ U.T @ b
print("Least-squares solution:", x_star)   # [1.0, 1.0] approx
print("Residual:", np.linalg.norm(A @ x_star - b))

# Compare with numpy's lstsq (uses SVD internally)
x_np, _, _, _ = np.linalg.lstsq(A, b, rcond=None)
print("numpy lstsq:", x_np)  # Same

# Pseudoinverse directly
A_pinv = np.linalg.pinv(A)   # Also uses SVD
print("A+ @ b =", A_pinv @ b)
print("A+ @ A = I_2?\n", np.round(A_pinv @ A, 10))  # 2x2 identity

One of the most beautiful applications of eigenvalue theory is to graphs. Given a graph G with n vertices, define:

• Adjacency matrix A: Aᵢⱼ = 1 if edge (i,j) exists, 0 otherwise.
• Degree matrix D: diagonal, Dᵢᵢ = degree of vertex i.
• Graph Laplacian L = D − A: a symmetric PSD matrix.

Why PSD? For any vector f: xᵀLx = Σ_{(i,j)∈E} (xᵢ − xⱼ)² ≥ 0. This is a sum of squared differences over edges — always non-negative. The Laplacian measures the 'smoothness' of a function f over the graph.

The spectrum of L tells you about graph structure:

• λ₁ = 0 always (constant vector is always an eigenvector)
• The number of zero eigenvalues = number of connected components
• Algebraic connectivity (Fiedler value) λ₂ > 0 iff the graph is connected. Larger λ₂ means a more 'robustly' connected graph.
• Fiedler vector (eigenvector for λ₂) reveals the graph's natural bipartition — used in spectral clustering and graph partitioning.
• Cheeger inequality: h(G)/2 ≤ sqrt(2λ₂) where h is the conductance (edge expansion). The Fiedler value certifies or refutes good connectivity.

Applications:

• Spectral clustering: embed vertices using the k smallest eigenvectors of L, then run k-means in this spectral space. Communities map to clusters.
• Graph neural networks: convolution on graphs generalizes Fourier analysis using the Laplacian eigenbasis.
• Laplacian eigenmaps: non-linear dimensionality reduction using the low-frequency eigenvectors of L as coordinates.
• PageRank: related to the dominant eigenvector of the (row-stochastic) adjacency matrix.

import numpy as np

# Triangle graph: 3 nodes, edges (0,1), (1,2), (0,2)
A = np.array([[0, 1, 1],
              [1, 0, 1],
              [1, 1, 0]], dtype=float)
D = np.diag(A.sum(axis=1))    # degree matrix
L = D - A                      # graph Laplacian

vals, vecs = np.linalg.eigh(L)
print("Laplacian eigenvalues:", np.round(vals, 4))  # [0, 3, 3]
print("Num components:", np.sum(np.abs(vals) < 1e-10))  # 1 (connected)
print("Fiedler value (lambda_2):", vals[1])  # 3.0 -> highly connected

# Two disconnected triangles
A2 = np.zeros((6, 6))
A2[:3, :3] = A
A2[3:, 3:] = A
D2 = np.diag(A2.sum(axis=1))
L2 = D2 - A2
vals2 = np.linalg.eigvalsh(L2)
print("Two triangles eigenvalues:", np.round(vals2, 4))
print("Num components:", np.sum(np.abs(vals2) < 1e-10))  # 2