If one theorem deserves the title “fundamental theorem of data science,” it is the singular value decomposition. The SVD factors any matrix — square or rectangular, full-rank or not — into a rotation, a scaling, and another rotation. It generalizes eigendecomposition to all matrices, exposes the four fundamental subspaces, gives the best low-rank approximation, defines the pseudoinverse, and underlies PCA, recommender systems, latent semantic analysis, and low-rank fine-tuning of large models. This is the summit of the core theory.

The decomposition

Every matrix \(\mA\in\R^{m\times n}\) factors as \[\mA = \mU\mSig\mV\tp,\] where \(\mU\in\R^{m\times m}\) and \(\mV\in\R^{n\times n}\) are orthogonal (columns orthonormal), and \(\mSig\in\R^{m\times n}\) is “diagonal” with nonnegative entries \(\sigma_1\ge\sigma_2\ge\cdots\ge0\) called the singular values. The columns of \(\mU\) are the left singular vectors, those of \(\mV\) the right singular vectors.

Read \(\mA\vx=\mU\mSig\mV\tp\vx\) right to left: \(\mV\tp\) rotates the input to align with the right singular vectors, \(\mSig\) stretches each axis by a singular value (and changes dimension), and \(\mU\) rotates the result into the output space. Every linear map, no matter how complicated, is just a rotation, an axis-aligned scaling, and another rotation. The SVD finds the two special orthonormal bases (one for the input space, one for the output) in which the matrix is diagonal.

Relation to eigendecomposition

The SVD is built from two symmetric matrices you already understand: \[\mA\tp\mA = \mV\mSig\tp\mSig\,\mV\tp,\qquad \mA\mA\tp = \mU\mSig\mSig\tp\mU\tp.\] So the right singular vectors \(\mV\) are the eigenvectors of \(\mA\tp\mA\), the left singular vectors \(\mU\) are the eigenvectors of \(\mA\mA\tp\), and the singular values are the square roots of the (shared, nonnegative) eigenvalues: \(\sigma_i=\sqrt{\lambda_i}\). Because \(\mA\tp\mA\) is symmetric PSD, this always exists — which is why the SVD exists for every matrix even when eigendecomposition fails.

The four subspaces, revealed

The SVD lays out all four fundamental subspaces with orthonormal bases. For rank \(r\) (the number of nonzero singular values):

Subspace Orthonormal basis from SVD
Column space \(\col(\mA)\) first \(r\) columns of \(\mU\)
Left null space \(\nullsp(\mA\tp)\) last \(m-r\) columns of \(\mU\)
Row space \(\col(\mA\tp)\) first \(r\) columns of \(\mV\)
Null space \(\nullsp(\mA)\) last \(n-r\) columns of \(\mV\)

The number of nonzero singular values is the rank — and, numerically, the number of singular values above a small threshold is the practical (“numerical”) rank.

Low-rank approximation: the Eckart–Young theorem

Writing the SVD as a sum of rank-one pieces, \[\mA = \sum_{i=1}^{r}\sigma_i\,\vu_i\vv_i\tp,\] and truncating after \(k\) terms gives \(\mA_k=\sum_{i=1}^k\sigma_i\vu_i\vv_i\tp\). This is the best possible rank-\(k\) approximation.

Among all matrices \(\mB\) of rank at most \(k\), the truncated SVD \(\mA_k\) minimizes the approximation error, and \[\norm{\mA-\mA_k}_2=\sigma_{k+1},\qquad \norm{\mA-\mA_k}_F=\sqrt{\sigma_{k+1}^2+\cdots+\sigma_r^2}.\] The proof rests on the orthogonal invariance of these norms: multiplying by an orthogonal matrix does not change them, so the problem reduces to the transparent diagonal case.

This one theorem is the mathematical core of data compression and dimensionality reduction. “Keep the largest singular values, discard the small ones” is provably the optimal way to approximate a matrix with a simpler one. The discarded tail (small \(\sigma_i\)) is usually noise; the retained head captures the signal. PCA, image compression, recommender systems, latent semantic analysis, and embedding compression are all this single move.

A grayscale image is a matrix of pixel intensities with a quickly-decaying singular-value spectrum. Reconstructing from the top \(k\) singular triplets stores only \(k(m+n+1)\) numbers instead of \(mn\), yet looks almost identical for modest \(k\) — a vivid demonstration that real data lives near a low-dimensional subspace.

The pseudoinverse

The SVD defines an inverse for any matrix, even non-square or singular ones.

The Moore–Penrose pseudoinverse of \(\mA=\mU\mSig\mV\tp\) is \[\mA^{+}=\mV\mSig^{+}\mU\tp,\] where \(\mSig^{+}\) inverts each nonzero singular value (and transposes the shape). The pseudoinverse gives the minimum-norm least-squares solution \(\hat\vx=\mA^{+}\vb\) of \(\mA\vx=\vb\): when there is no exact solution it returns the least-squares fit, and when there are infinitely many it returns the smallest one.

The pseudoinverse unifies the whole book. For an invertible square matrix it is the ordinary inverse. For a tall full-rank matrix it equals \((\mA\tp\mA)\inv\mA\tp\) — the least-squares operator of Chapter [ch:orthogonality]. For a wide or rank-deficient matrix it picks the minimum-norm solution from the infinite family — exactly what L2 regularization approaches as \(\lambda\to0\). One formula, every case.

Compute least squares and the pseudoinverse via the SVD (, ), not by forming \((\mA\tp\mA)\inv\). Forming \(\mA\tp\mA\) squares the condition number \(\sigma_1/\sigma_r\), so you lose up to twice as many digits of accuracy. This is the same reason production PCA uses the SVD of the data matrix rather than the eigendecomposition of the covariance matrix.

The SVD is the most useful single tool in applied linear algebra for ML. PCA is the SVD of centered data: principal components are right singular vectors, explained variance comes from \(\sigma_i^2\). Recommender systems factor the user–item matrix as a low-rank product (the Netflix-Prize approach). Latent semantic analysis applies truncated SVD to term–document matrices. Embedding compression for vector databases is a truncated SVD/PCA step before quantization (8–32\(\times\) smaller). And LoRA, the dominant way to fine-tune large language models, freezes the pretrained weights and learns a low-rank update \(\Delta\mW=\mB\mA\) — a direct bet that the needed change is low-rank, exactly the Eckart–Young intuition.

import numpy as np
A = np.random.randn(8, 5)
U, s, Vt = np.linalg.svd(A, full_matrices=False)   # reduced SVD
print(np.allclose(A, (U * s) @ Vt))                # True
k = 2
A_k = (U[:, :k] * s[:k]) @ Vt[:k, :]               # best rank-2 approx
print(np.linalg.norm(A - A_k, 2), s[k])            # equals sigma_{k+1}
x = np.linalg.pinv(A) @ np.random.randn(8)         # min-norm least squares

Chapter summary

  • Every matrix factors as \(\mA=\mU\mSig\mV\tp\): rotate, stretch by singular values, rotate.

  • \(\mV\) are eigenvectors of \(\mA\tp\mA\), \(\mU\) of \(\mA\mA\tp\), and \(\sigma_i=\sqrt{\lambda_i}\); the SVD always exists.

  • Nonzero singular values \(=\) rank; \(\mU,\mV\) give orthonormal bases for all four subspaces.

  • Eckart–Young: the truncated SVD \(\mA_k\) is the best rank-\(k\) approximation, with error \(\sigma_{k+1}\).

  • The pseudoinverse \(\mA^+=\mV\mSig^+\mU\tp\) gives minimum-norm least squares — one formula for all cases.

  • SVD powers PCA, recommender systems, LSA, embedding compression, and LoRA.