Singular Value Decomposition

What is the Singular Value Decomposition?

The Singular Value Decomposition (SVD) is the most powerful and general matrix factorization in linear algebra. Unlike LU decomposition (which requires a square matrix) or QR decomposition (which requires linearly independent columns), SVD works for any matrix of any shape and any rank. Every m × n matrix A — even a singular or rectangular one — can be written as A = UΣVᵀ, where U and V are orthogonal matrices and Σ is a diagonal matrix of non-negative numbers.

The entries on the diagonal of Σ are called the singular values of A, conventionally ordered from largest to smallest: σ₁ ≥ σ₂ ≥ … ≥ σₗ ≥ 0, where l = min(m, n) is the number of singular values. The columns of U are the left singular vectors and the columns of V are the right singular vectors. Together, they provide a complete geometric picture of what A does to any vector: it first rotates (or reflects) via Vᵀ, then stretches along r independent axes by the singular values, then rotates (or reflects) again via U.

SVD reveals the true structure of a matrix more completely than any other factorization. It exposes the rank, the condition number, the four fundamental subspaces (column space, row space, null space, and left null space), and the best low-rank approximation — all from a single decomposition. It is the mathematical foundation of principal component analysis, image compression, latent semantic analysis, and dimensionality reduction.

The Three Factors: U, Σ, and V

For an m × n matrix A, the full SVD writes A = UΣVᵀ where:

U is an m × m orthogonal matrix: UᵀU = UUᵀ = Iₘ. Its columns u₁, u₂, …, uₘ are the left singular vectors.
Σ is an m × n diagonal matrix (padded with zeros) with non-negative diagonal entries σ₁ ≥ σ₂ ≥ … ≥ σₗ ≥ 0 where l = min(m, n). These are the singular values of A.
V is an n × n orthogonal matrix: VᵀV = VVᵀ = Iₙ. Its columns v₁, v₂, …, vₙ are the right singular vectors.

Singular Value Decomposition

A = U \Sigma V^T, \quad U^T U = I_m, \quad V^T V = I_n, \quad \Sigma = \mathrm{diag}(\sigma_1, \sigma_2, \ldots, \sigma_l), \quad \sigma_1 \geq \sigma_2 \geq \cdots \geq \sigma_l \geq 0

A is m × n, U is m × m orthogonal, Σ is m × n diagonal with σ₁ ≥ σ₂ ≥ … ≥ 0 on the diagonal, and V is n × n orthogonal. The number of non-zero singular values equals the rank r of A. The "thin" or "economy" SVD retains only the r non-zero singular values and the corresponding r columns of U and V.

There is also a "thin" or "economy" SVD, in which U is trimmed to m × r (keeping only the r left singular vectors that correspond to non-zero singular values), Σ is trimmed to r × r (keeping only the non-zero diagonal block), and V is trimmed to n × r. The thin SVD satisfies A = UΣVᵀ just as well as the full SVD but wastes no space on the zero singular values. In practice, numpy.linalg.svd and scipy.linalg.svd both offer a full_matrices flag to choose between full and thin forms.

Geometric Interpretation: Rotation → Stretch → Rotation

The SVD A = UΣVᵀ has a beautiful geometric meaning. When A acts on a vector x, the transformation happens in three stages:

Vᵀ rotates (or reflects) the input. Because V is orthogonal, multiplication by Vᵀ is a rigid motion — it preserves all lengths and angles. The right singular vectors v₁, …, vₙ define a coordinate system in the input space.
Σ stretches along coordinate axes. After rotating into the right singular vector basis, each coordinate is scaled by the corresponding singular value σᵢ. Directions with large σᵢ are amplified; directions with σᵢ = 0 are collapsed to zero.
U rotates (or reflects) the output. The left singular vectors u₁, …, uₘ define a coordinate system in the output space. Multiplication by U is again a rigid motion.

Outer Product Expansion

A = \sum_{i=1}^{r} \sigma_i \, u_i v_i^T = \sigma_1 u_1 v_1^T + \sigma_2 u_2 v_2^T + \cdots + \sigma_r u_r v_r^T

Every matrix A is the sum of r rank-one matrices σᵢuᵢvᵢᵀ, each built from one pair of singular vectors scaled by the corresponding singular value. The first term σ₁u₁v₁ᵀ captures the dominant "direction" of A; successive terms add finer and finer detail. Truncating this sum at k terms gives the best rank-k approximation to A.

This outer product view is extremely useful. It says that A "acts on" any vector x by projecting x onto each right singular vector vᵢ (to get the component vᵢᵀx), scaling by σᵢ, and placing the result along the left singular vector uᵢ in the output space. The most important directions of action are those with the largest singular values.

SVD and the Four Fundamental Subspaces

The SVD directly reveals all four fundamental subspaces of A that Gilbert Strang famously organized into his "big picture" of linear algebra:

Column space (range) of A: spanned by the left singular vectors u₁, …, uᵣ corresponding to non-zero singular values. Dimension = r = rank(A).
Left null space of A: spanned by the left singular vectors uᵣ₊₁, …, uₘ corresponding to zero singular values. Dimension = m − r.
Row space of A: spanned by the right singular vectors v₁, …, vᵣ. Dimension = r.
Null space (kernel) of A: spanned by the right singular vectors vᵣ₊₁, …, vₙ. Dimension = n − r.

The rank-nullity theorem (r + (n − r) = n) is immediately visible: the row space and null space together partition ℝⁿ into two orthogonal complements. Similarly, the column space and left null space partition ℝᵐ. The SVD makes this explicit by naming orthonormal bases for all four subspaces simultaneously.

Singular Values, Eigenvalues, and the Condition Number

Singular values are related to — but different from — eigenvalues. For a general m × n matrix A, the singular values σᵢ are the square roots of the eigenvalues of the symmetric positive semidefinite matrix AᵀA (or AAᵀ). If A is itself symmetric positive semidefinite, then its singular values equal its eigenvalues, and the SVD reduces to the eigendecomposition.

Singular Values and Eigenvalues

A^T A = V \Sigma^T U^T U \Sigma V^T = V \Sigma^2 V^T, \quad \sigma_i = \sqrt{\lambda_i(A^T A)}, \quad \kappa_2(A) = \frac{\sigma_1}{\sigma_\ell}, \;\; \ell = \min(m,n)

The right singular vectors V are the eigenvectors of AᵀA; the left singular vectors U are the eigenvectors of AAᵀ; the singular values satisfy σᵢ² = λᵢ(AᵀA). The two-norm ‖A‖₂ = σ₁ (largest singular value), and the condition number κ₂(A) = σ₁/σₗ with l = min(m, n) (ratio of largest to smallest singular value). A rank-deficient matrix has σₗ = 0, so κ₂ = ∞; a well-conditioned matrix has σ₁/σₗ close to 1.

The condition number κ₂(A) = σ₁/σₗ, with l = min(m, n), is the key to understanding how sensitive a linear system Ax = b is to perturbations. If κ₂ is large, small changes in b (e.g., measurement noise) cause large changes in x. The SVD makes the condition number explicit and computable, and it tells you exactly which directions in the input and output space are well-conditioned (large σᵢ) versus ill-conditioned (small σᵢ). When A is rank-deficient, σₗ = 0 and κ₂ = ∞; in that case what people usually quote is the effective condition number σ₁/σᵣ, restricted to the r = rank(A) nonzero singular values. The two coincide exactly when A is square and full rank.

Low-Rank Approximation: The Eckart–Young Theorem

The most celebrated application of SVD is optimal low-rank approximation. Given a matrix A of rank r, the best approximation of A by a rank-k matrix (for any k ≤ r) in both the Frobenius norm and the two-norm is the truncated SVD:

Eckart–Young Theorem

A_k = \sum_{i=1}^{k} \sigma_i u_i v_i^T, \quad \|A - A_k\|_2 = \sigma_{k+1}, \quad \|A - A_k\|_F = \sqrt{\sigma_{k+1}^2 + \cdots + \sigma_r^2}

Aₖ = Σᵢ₌₁ᵏ σᵢuᵢvᵢᵀ is the best rank-k approximation to A. The approximation error in the Frobenius norm is ‖A − Aₖ‖_F = √(σ²_{k+1} + … + σ²_r), and in the two-norm is ‖A − Aₖ‖₂ = σ_{k+1}. No rank-k matrix is closer to A in either norm. This is optimal compression: you capture the most important "directions" of A with the fewest terms.

The Eckart–Young theorem (1936) guarantees that the truncated SVD is the best possible compression of A — you cannot do better with any other rank-k matrix. In practice, the "explained variance" of the truncated SVD is measured by the fraction of total squared singular-value mass retained: (σ₁² + … + σₖ²) / (σ₁² + … + σᵣ²). For real data matrices, the singular values often decay rapidly, so a small k can capture 95% or more of the total variance.

Image Compression via SVD

A grayscale image is naturally represented as a matrix A of pixel intensities. An m × n image has mn numbers. The rank-k SVD approximation Aₖ stores only k(m + n + 1) numbers (k left singular vectors, k singular values, k right singular vectors) — a compression ratio of mn / [k(m + n + 1)]. For a 1000 × 1000 image with k = 50, this is 1,000,000 vs. 100,050 numbers — a 10× compression — while retaining the bulk of visual information.

Computing the SVD

Computing the full SVD is a two-phase process. In the first phase, A is reduced to bidiagonal form B = UᵀAV using Householder reflections (similar to how Householder reduces a matrix to Hessenberg form for eigenvalue computation). The bidiagonal matrix B has non-zero entries only on the main diagonal and the superdiagonal. This phase costs O(mn²) flops for m ≥ n.

In the second phase, the singular values of B are computed by an iterative algorithm (the Golub–Reinsch algorithm, a variant of the QR algorithm applied to B). Each step of the iteration deflates one singular value pair until all are found. In practice, the algorithm converges rapidly with proper shifts, and the overall cost of the full SVD is O(mn² + n³) — dominated by the bidiagonalization.

For very large matrices, computing all singular values is prohibitive. Randomized SVD algorithms (e.g., the Halko–Martinsson–Tropp method) find an approximate truncated SVD in O(mn log k + (m + n)k²) time by randomly projecting A to a lower-dimensional subspace, computing the exact SVD of the small projected matrix, and lifting back. For a 10,000 × 10,000 matrix and k = 100, randomized SVD is roughly 100× faster than the full SVD and produces approximation errors that are nearly optimal with high probability.

Applications of SVD

Principal Component Analysis (PCA)

PCA is SVD applied to centered data. Given a data matrix X with m observations and n features, center the data to get X̃ = X − mean(X). The SVD X̃ = UΣVᵀ gives the principal components as the columns of V (right singular vectors). The variance explained by each component is σᵢ²/(m − 1). PCA finds the directions of maximum variance in the data — equivalently, the directions along which the data "spreads out" the most. It is the go-to method for exploratory data analysis, visualization in 2D/3D, feature extraction, and denoising.

Recommender Systems

A user–item ratings matrix R (rows = users, columns = movies) is typically very sparse and low-rank: most movies are rated by few users, and user preferences cluster into a small number of "latent factors" (genres, moods, etc.). The truncated SVD Rₖ = UₖΣₖVₖᵀ finds these latent factors and can predict missing ratings by filling in the low-rank approximation. This is the mathematical core of collaborative filtering — the algorithm behind Netflix and Spotify recommendations.

Latent Semantic Analysis

In natural language processing, a term–document matrix T has one row per word and one column per document, with entry T_{ij} equal to the frequency of word i in document j. The low-rank SVD Tₖ finds "latent semantic" topics — directions in the word-space that capture recurring co-occurrence patterns. Documents with similar semantic content (even if they use different words) map to nearby vectors in the latent space. This is the foundation of latent semantic indexing, a precursor to modern word embeddings.

Pseudoinverse and Least Squares

The Moore–Penrose pseudoinverse A⁺ of any matrix A is defined via SVD as A⁺ = VΣ⁺Uᵀ, where Σ⁺ replaces each non-zero σᵢ with 1/σᵢ and leaves zeros unchanged. The pseudoinverse gives the minimum-norm least-squares solution x = A⁺b: among all vectors x that minimize ‖Ax − b‖, it picks the one with smallest ‖x‖. For full-rank square matrices, A⁺ = A⁻¹. For over- or underdetermined systems, A⁺ provides the "best possible" solution in a well-defined sense.

Moore–Penrose Pseudoinverse

A^+ = V \Sigma^+ U^T, \quad \Sigma^+_{ii} = \begin{cases} 1/\sigma_i & \sigma_i > 0 \\ 0 & \sigma_i = 0 \end{cases}, \quad x^* = A^+ b

A⁺ = VΣ⁺Uᵀ where Σ⁺ inverts non-zero singular values. The minimum-norm least-squares solution to Ax = b is x* = A⁺b = VΣ⁺Uᵀb. This generalizes matrix inversion to any matrix, square or not, full-rank or not. When A has full column rank, A⁺ = (AᵀA)⁻¹Aᵀ — the normal equations solution, computed stably via SVD.

Numerical Properties of SVD

SVD is numerically the most reliable of all matrix factorizations. The singular values are computed to near machine precision even for ill-conditioned matrices, and the computed U and V are orthogonal to within machine precision. This reliability comes at a price — SVD is about 6× more expensive than QR decomposition for the same matrix — but for problems where accuracy is paramount, the extra cost is justified.

One important use of SVD is determining the numerical rank of a matrix: rather than counting the number of non-zero singular values (which is fragile in floating point), one counts the number of singular values above a threshold τ = ε · σ₁, where ε is the machine epsilon. This "numerical rank" is the appropriate notion of rank for real data corrupted by floating-point errors or measurement noise.

Engineering Implication

SVD is ubiquitous in engineering: numpy.linalg.svd and scipy.linalg.svd compute it in Python; MATLAB's svd is a standard tool. sklearn.decomposition.TruncatedSVD and sklearn.decomposition.PCA use it under the hood. In signal processing, SVD-based methods appear in MUSIC and ESPRIT direction-of-arrival algorithms, Wiener filtering, and adaptive beamforming. In control theory, the Hankel singular values of a dynamical system (computed via SVD) determine which states can be truncated in model order reduction. In machine learning, SVD underlies word2vec and GloVe word embeddings, matrix factorization recommender systems, and the attention mechanism in transformers (which computes QKV projections — effectively a learned low-rank decomposition of the attention matrix).

Key Takeaways

SVD factors any matrix A = UΣVᵀ into two orthogonal matrices and a diagonal matrix of singular values. It reveals rank, condition number, and all four fundamental subspaces. The Eckart–Young theorem guarantees that the truncated SVD is the best low-rank approximation in any unitarily invariant norm. The pseudoinverse A⁺ = VΣ⁺Uᵀ solves any linear system in the minimum-norm least-squares sense. Applications span image compression, PCA, recommender systems, latent semantic analysis, and pseudoinverse computation. SVD is the most reliable and most informative matrix factorization in numerical linear algebra.