Eigendecomposition works only for square matrices, and nicely only for symmetric ones. The Singular Value Decomposition (SVD) works for every matrix. Gilbert Strang calls it the climax of linear algebra, and in machine learning it is everywhere: PCA, latent semantic analysis, recommender systems, model compression, numerical least squares and the analysis of neural network weights.
The theorem#
Any real matrix $\mathbf{A} \in \mathbb{R}^{m \times n}$ of rank $r$ can be written as
where
- $\mathbf{U} \in \mathbb{R}^{m \times m}$ is orthogonal — its columns $\mathbf{u}_i$ are left singular vectors;
- $\mathbf{V} \in \mathbb{R}^{n \times n}$ is orthogonal — its columns $\mathbf{v}_i$ are right singular vectors;
- $\boldsymbol{\Sigma} \in \mathbb{R}^{m \times n}$ is diagonal with singular values $\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_r > 0$.
Equivalently, as a sum of rank-one pieces:
Geometric meaning#
Every linear map, however complicated, does three simple things in sequence: rotate (by $\mathbf{V}^\top$), stretch along the axes (by $\boldsymbol{\Sigma}$), and rotate again (by $\mathbf{U}$). The unit sphere is mapped to an ellipsoid whose semi-axes have lengths $\sigma_i$ in directions $\mathbf{u}_i$.
Connection to eigenvalues#
So the right singular vectors are eigenvectors of $\mathbf{A}^\top\mathbf{A}$, the left singular vectors are eigenvectors of $\mathbf{A}\mathbf{A}^\top$, and $\sigma_i^2$ are their eigenvalues. (In practice we never form $\mathbf{A}^\top\mathbf{A}$ explicitly — it squares the condition number and loses precision.)
The Eckart–Young theorem: best low-rank approximation#
Truncate the SVD to the top $k$ terms:
Theorem (Eckart–Young–Mirsky). Among all matrices of rank at most $k$, $\mathbf{A}_k$ is the closest to $\mathbf{A}$ in both the spectral and Frobenius norms:
This is the mathematical basis of every "keep the important directions, drop the noise" technique.
Application 1: image compression#
A grayscale image is a matrix. Storing a rank-$k$ approximation needs $k(m + n + 1)$ numbers instead of $mn$.
import numpy as np
rng = np.random.default_rng(0)
# A synthetic "image": smooth structure plus noise
x = np.linspace(0, 1, 200)
img = np.outer(np.sin(6 * x), np.cos(4 * x)) + 0.5 * np.outer(x, x**2) + 0.05 * rng.normal(size=(200, 200))
U, s, Vt = np.linalg.svd(img, full_matrices=False)
for k in [1, 2, 5, 20]:
approx = (U[:, :k] * s[:k]) @ Vt[:k]
err = np.linalg.norm(img - approx) / np.linalg.norm(img)
print(f"rank {k:>2}: relative error {err:.3f}, storage {k * (200 + 200 + 1)} vs {200 * 200}")Application 2: PCA#
Centre the data matrix $\mathbf{X}$ (subtract column means). Its SVD $\mathbf{X} = \mathbf{U}\boldsymbol{\Sigma}\mathbf{V}^\top$ directly gives the principal directions (columns of $\mathbf{V}$) and the variance explained by each ($\sigma_i^2/(n-1)$). This is how scikit-learn implements PCA.
Application 3: recommender systems#
A user–item rating matrix is approximately low rank: tastes are driven by a few latent factors (genre, mood, price sensitivity). Low-rank factorisation — conceptually an SVD of a matrix with missing entries — powered the winning approaches in the Netflix Prize.
Application 4: least squares and the pseudoinverse#
The Moore–Penrose pseudoinverse is
where $\boldsymbol{\Sigma}^+$ inverts the non-zero singular values. Then $\mathbf{x} = \mathbf{A}^+\mathbf{b}$ is the minimum-norm least-squares solution of $\mathbf{Ax} \approx \mathbf{b}$, even when $\mathbf{A}$ is rank-deficient.
Application 5: conditioning and neural networks#
The condition number $\kappa(\mathbf{A}) = \sigma_{\max}/\sigma_{\min}$ measures sensitivity of $\mathbf{Ax} = \mathbf{b}$ to perturbations. In deep learning, the singular values of weight matrices govern how signals and gradients grow through layers. Spectral normalisation, used to stabilise GAN training, divides a weight matrix by $\sigma_{\max}$. And LoRA fine-tuning exploits the empirical observation that useful weight updates are approximately low rank.