Singular Value Decomposition

Singular Value Decomposition (SVD) is a matrix factorization that decomposes any real or complex matrix A into three structured matrices — U, Σ, and Vᵀ — revealing the intrinsic geometric action of the linear map as a rotation, scaling, and second rotation. It is arguably the most important matrix decomposition in modern AI, underpinning dimensionality reduction, recommender systems, natural language semantics, and the theoretical analysis of neural network weight matrices. The counterintuitive fact: SVD exists and is unique (up to sign conventions) for every matrix, regardless of shape or rank — even a matrix of all zeros has a valid SVD.

The core insight

Every matrix — no matter how chaotic its entries — is secretly just a machine that stretches space along a set of orthogonal axes by specific amounts, with a rotation before and after. SVD makes this structure explicit: the columns of U are the output axes, the columns of V are the input axes, and the diagonal entries of Σ (the singular values) tell you exactly how much stretching happens along each axis. The "aha moment" is realizing that the rank of a matrix equals the number of nonzero singular values, and that by keeping only the top k singular values you get the best possible rank-k approximation of the original matrix — a profound compression guarantee no other decomposition provides.

How It Works

For a matrix A of shape m × n, SVD produces the factorization:

1
A=UΣV⊤\quad A = U \Sigma V^\top

Explanation: The matrix A is expressed as the product of an m×m orthogonal matrix U, an m×n diagonal matrix Σ with non-negative entries, and the transpose of an n×n orthogonal matrix V — together encoding a rotation, stretch, and counter-rotation. Here every symbol earns its place: U ∈ ℝ^{m×m} whose columns {u₁, …, u_m} are called the left singular vectors; Σ ∈ ℝ^{m×n} with diagonal entries σ₁ ≥ σ₂ ≥ … ≥ σ_r ≥ 0 (the singular values, ordered by magnitude, where r = rank(A)); and V ∈ ℝ^{n×n} whose columns {v₁, …, v_n} are the right singular vectors. Both U and V are orthogonal: UᵀU = I and VᵀV = I. The singular values are derived from the eigenvalues of the symmetric positive-semidefinite matrix AᵀA:

2
A⊤A vi=σi2 vi\quad A^\top A \, v_i = \sigma_i^2 \, v_i

Explanation: The right singular vectors vᵢ are eigenvectors of AᵀA, and the singular values σᵢ are the square roots of those eigenvalues — which are always non-negative because AᵀA is positive semidefinite. Similarly, the left singular vectors satisfy AAᵀ uᵢ = σᵢ² uᵢ, and the two sets are linked by:

3
ui=1σiAvi,vi=1σiA⊤ui(σi>0)\quad u_i = \frac{1}{\sigma_i} A v_i, \qquad v_i = \frac{1}{\sigma_i} A^\top u_i \quad (\sigma_i > 0)

Explanation: Once you have the right singular vectors from diagonalizing AᵀA, you can recover the corresponding left singular vectors by simply applying A and normalizing — the two matrix families are not independent. A crucial identity is the dyadic expansion (also called the outer-product form), which expresses A as a sum of rank-1 matrices:

4
A=∑i=1rσi uivi⊤\quad A = \sum_{i=1}^{r} \sigma_i \, u_i v_i^\top

Explanation: The full matrix is the weighted sum of r rank-1 "layers" uᵢvᵢᵀ, each scaled by its singular value σᵢ — truncating this sum at k < r terms gives the optimal rank-k approximation in both the Frobenius norm and the spectral norm. The optimality of truncated SVD is guaranteed by the Eckart–Young–Mirsky theorem:

5
A^k=∥A−B∥F=∑i=1kσi uivi⊤\quad \hat{A}_k = \underset{\text{rank}(B) \leq k}{\arg\min} \| A - B \|_F = \sum_{i=1}^{k} \sigma_i \, u_i v_i^\top

Explanation: No rank-k matrix is closer to A in Frobenius norm than the k-truncated SVD — this is the mathematical bedrock that makes PCA, latent semantic analysis, and matrix completion methods theoretically sound. ```python import torch

# Full SVD in PyTorch — returns U, S, Vh (note: Vh = V^T) A = torch.randn(100, 50) U, S, Vh = torch.linalg.svd(A, full_matrices=True)

# U: (100,100), S: (50,), Vh: (50,50)

# Economical (thin) SVD — only r = min(m,n) singular vectors U_thin, S, Vh_thin = torch.linalg.svd(A, full_matrices=False)

# U_thin: (100,50), S: (50,), Vh_thin: (50,50)

# Rank-k approximation using dyadic expansion k = 10 A_k = (U_thin[:, :k] * S[:k]) @ Vh_thin[:k, :]

# Verify optimality: Frobenius error equals sqrt(sum of dropped σ²) residual_theory = torch.sqrt((S[k:] ** 2).sum()) residual_actual = torch.norm(A - A_k, p='fro') print(f"Theory: {residual_theory:.4f} | Actual: {residual_actual:.4f}")

# These will match exactly (up to float precision)

# Fraction of variance captured by top-k components variance_explained = (S[:k]2).sum() / (S2).sum() print(f"Variance explained by k={k}: {variance_explained:.2%}") ```

researcher · Gene Howard Golub | 1932 – 2007 | Numerically Stable SVD Algorithm

Gene Golub was an American numerical analyst at Stanford University who spent his career turning mathematically elegant decompositions into computationally reliable algorithms. His landmark 1965 paper with William Kahan introduced the bidiagonalization approach that made computing the SVD of large matrices numerically stable and practical — the Golub–Reinsch algorithm remains the backbone of SVD implementations in LAPACK, NumPy, and PyTorch to this day. A surprising personal fact: Golub was famously known for calculating matrix decompositions mentally during seminars, and the term "Golub's rule" — that the SVD should be tried before any other decomposition — became an unofficial maxim in numerical linear algebra.

Figure
SVD Geometric Decomposition
A = U Σ Vᵀ — Geometric Interpretation Input Space ℝⁿ Rotate by Vᵀ Scale by Σ Output Space ℝᵐ v₁ v₂ Vᵀ ê₁ ê₂ Σ σ₁ σ₂ U u₁ u₂ Singular values σ₁ ≥ σ₂ ≥ 0 encode how much each axis is stretched
Figure
Truncated SVD — Rank-k Approximation
Rank-k Truncation: Optimal Low-Rank Approximation A m × n rank r = U m × r Σ r × r Vᵀ r × n ≈ U_k m × k Σ_k k × k Vₖᵀ k × n Best rank-k approx (Eckart–Young theorem) Error = ‖A − Â_k‖_F = √(σ²_{k+1} + … + σ²_r)

Intuition

The Mental Model

Imagine you are a photographer analyzing thousands of face images, each one a vector in a very high-dimensional space. SVD asks: "Is there a small number of prototypical face components — a component for lighting direction, one for face width, one for smile degree — such that every face is just a weighted mixture of these?" The right singular vectors (V columns) are those prototype components in the original feature space. The left singular vectors (U columns) tell you how each face scores on each prototype. The singular values tell you which prototypes matter most — a large σ₁ means the first prototype explains a huge chunk of variation across all faces, while a tiny σ₅₀ means the fiftieth prototype adds almost nothing. Truncating to top k prototypes is like saying: "I'll describe every face using only the 10 most important components and accept a small reconstruction error" — and SVD guarantees you cannot find any other 10 components that do it better. The further geometric picture: think of A as a machine that turns input unit sphere into an output ellipsoid; the singular values are the semi-axis lengths of that ellipsoid, and U, V are the axis orientations in output and input space respectively.

Historical Development

Origins · 1 minhow this was discovered
History
1873
Eugenio Beltrami publishes the first known treatment of the SVD for square bilinear forms, establishing the core factorization idea independently for real matrices.
1874
Camille Jordan independently derives the same decomposition for real square matrices, working from a different geometric tradition in France.
1907
Erhard Schmidt extends the decomposition to integral operators in Hilbert spaces, giving the continuous functional-analysis generalization that would later underpin kernel methods.
1936
Carl Eckart and Gale Young prove the landmark approximation theorem bearing their names — demonstrating that truncated SVD gives the provably optimal low-rank approximation in Frobenius and spectral norms.
1965
Gene Golub and William Kahan introduce the bidiagonalization algorithm, making SVD numerically stable and computationally practical for large matrices for the first time.
1970
Golub and Christian Reinsch refine this into the Golub–Reinsch algorithm, which is standardized in EISPACK and later LAPACK and becomes the workhorse implementation used worldwide.
1988
Scott Deerwester and colleagues apply truncated SVD to term-document matrices under the name Latent Semantic Analysis (LSA), demonstrating that SVD can capture semantic similarity between words — an early precursor to word embeddings.
2006
Yoshua Bengio, Yann LeCun, and Geoffrey Hinton's work on deep representations spurs interest in SVD as a lens for understanding weight matrix rank and intrinsic dimensionality of learned representations.
2019
The paper "The Lottery Ticket Hypothesis" (Frankle & Carlin) and subsequent work on neural network pruning formalizes SVD-based low-rank approximation as a principled neural compression technique, used in production at companies including Google and Meta.
2023
LoRA (Low-Rank Adaptation) and its variants — which decompose weight matrix updates as products of two thin matrices, exactly approximating a rank-k SVD update — become the dominant technique for fine-tuning large language models such as LLaMA and GPT-4 class models.

Why It Matters

Real-World Impact
Principal Component Analysis (PCA)PCA is exactly the SVD of the mean-centered data matrix — the right singular vectors V are the principal components, and the singular values σᵢ divided by √(n−1) are the standard deviations along each principal direction. Every call to `sklearn.decomposition.PCA` internally calls a truncated SVD routine.
Recommender Systems (Netflix, Spotify)Matrix factorization for collaborative filtering decomposes the user–item rating matrix using approximate SVD, identifying latent factors like "this user prefers jazz" or "this item is acoustic". The Netflix Prize–winning solution used SVD-inspired factorization achieving a 10% RMSE reduction over their baseline.
LoRA Fine-Tuning of LLMs (LLaMA-2, Mistral, GPT-4)Low-Rank Adaptation (Hu et al., 2022) represents weight matrix updates ΔW ≈ BA where B ∈ ℝ^{m×r} and A ∈ ℝ^{r×n} with r ≪ min(m,n) — a trainable rank-r SVD-like factorization. This reduces trainable parameters by 10,000× for 7B-parameter models while preserving 95%+ of full fine-tuning quality.
Latent Semantic Analysis & Word EmbeddingsLSA applies truncated SVD to a term–document TF-IDF matrix (vocabulary × documents), producing dense word vectors where cosine similarity captures semantic relatedness — predating Word2Vec by 25 years. The singular vectors encode latent topics without any supervision.
Numerical Linear Algebra & Least SquaresThe Moore–Penrose pseudoinverse A⁺ = VΣ⁺Uᵀ (where Σ⁺ inverts each nonzero singular value) gives the minimum-norm least-squares solution to any linear system Ax = b, even rank-deficient ones. PyTorch's `torch.linalg.lstsq` uses this internally.
Neural Network Compression (BERT, Vision Transformers)Production inference systems compress transformer attention and FFN layers by replacing weight matrices with their rank-k SVD approximations, achieving 2–4× model size reduction with under 1% accuracy degradation — deployed in Google's on-device ML pipelines.
Signal Denoising and Image CompressionHard thresholding small singular values of an image matrix (treated as a 2D array of pixel intensities) removes noise while preserving structure — the basis for classical image compression techniques and still used in scientific imaging (MRI, astronomical imaging).

Prerequisites

What This Leads To

See Also

Principal Component Analysis · Matrix Factorization · Eigendecomposition · Low-Rank Approximation · Latent Semantic Analysis

Where this leads