Matrix Decompositions

A decomposition rewrites a matrix as a product of simpler factors, each carrying one kind of structure: rotations and stretches (SVD), orthogonal bases (QR), square roots of positive-definite matrices (Cholesky), or the spectrum (eigen). Decompositions are the workhorses behind least squares, dimensionality reduction, regularization, and conditioning analysis - they turn "what does this matrix do" into "which factors, in what order." NumPy ships the standard set through np.linalg, backed by LAPACK.

SVD: The Universal One

import numpy as np

A = np.random.default_rng(0).normal(size=(6, 4))
U, s, Vt = np.linalg.svd(A, full_matrices=False)
print(U.shape, s.shape, Vt.shape)   # (6,4) (4,) (4,4) - the reduced form

print((U @ np.diag(s) @ Vt) - A).max()  # ~0 - the factorization holds

# rank and pseudoinverse from the singular values
rank = (s > 1e-10).sum()
Ainv = np.linalg.pinv(A)
print(Ainv.shape)          # (4, 6) - the least-squares inverse
print((A @ Ainv @ A) - A).max()   # ~0 - it reconstructs A

The SVD writes A = U S Vt with orthogonal U and Vt and non-negative singular values s. The singular values are the matrix's stretch spectrum - the largest is the operator norm, the count above a tolerance is the numeric rank, and truncating to the top-k factors is the best rank-k approximation of A. pinv (built on the SVD) is the least-squares inverse for any shape - the tool for underdetermined or rank-deficient systems. If you learn one decomposition, make it this one.

Cholesky: The Positive-Definite Fast Path

# a symmetric positive-definite matrix - the common SPD shape
X = np.random.default_rng(1).normal(size=(30, 5))
C = X.T @ X + 0.1 * np.eye(5)     # Gram matrix plus jitter: SPD

L = np.linalg.cholesky(C)         # C = L L.T with L lower triangular
print(np.allclose(C, L @ L.T))    # True

# solving with an SPD matrix: factor once, solve twice - cheaper than LU
b = np.random.default_rng(2).normal(size=5)
x = np.linalg.solve(L.T, np.linalg.solve(L, b))
print((C @ x) - b).max()          # ~0

Cholesky factors a symmetric positive-definite matrix into L L.T at roughly half the cost of a general factorization, and it doubles as the positive-definiteness test: if cholesky raises, the matrix is not SPD. The double-triangular-solve pattern (solve L, then L.T) is the efficient way to apply the inverse - the standard move for covariance inversion in statistics. Adding a small diagonal jitter (0.1 * eye) is the usual fix for a near-singular Gram matrix.

QR and Eigen: The Others

A = np.random.default_rng(3).normal(size=(5, 3))
Q, R = np.linalg.qr(A)
print(Q.T @ Q)              # ~identity - orthonormal columns
print((Q @ R) - A).max()    # ~0

S = np.array([[2.0, 1.0], [1.0, 2.0]])
w, V = np.linalg.eigh(S)    # symmetric: real values, orthonormal V
print(w)                    # [1. 3.]
print((V.T @ (S @ V)) - np.diag(w)).max()   # ~0 - diagonalized

qr produces an orthogonal factor - the stable route to least squares and the Gram-Schmidt idea without its numerical instability. For symmetric matrices, eigh diagonalizes: S = V diag(w) V.T, which is the basis of variance decomposition (PCA is exactly this, applied to a covariance matrix). Use the general eig only for genuinely non-symmetric problems, where eigenvalues can be complex.

Which Decomposition, Which Question

Summary

Decompositions factor a matrix into structured pieces, and each answers a different question: SVD for rank, norms, pseudoinverses, and low-rank approximation (the universal tool); Cholesky as the fast path and the positive-definiteness test for SPD matrices; QR for orthonormal bases and stable least squares; eigh for the symmetric spectrum. The habit: identify the question first, pick the factor that encodes the relevant structure, and verify the factorization numerically - the product of the factors minus the original, near zero - because a decomposition is only as trustworthy as the roundoff it carries.