Prev
Unit 5 · Least Squares and Data Fitting
The Singular Value Decomposition
65m100 XP
Next
📐Least Squares and Data Fitting

The Singular Value Decomposition

A = UΣVᵀ, singular values as gains, rank and pseudoinverse, low-rank approximation, regularisation

The Netflix Prize, 2006–2009
In October 2006 Netflix releases 100 million ratings — 480,189 customers, 17,770 films, one to five stars — and offers $1,000,000 to the first team that predicts held-out ratings 10% more accurately than its own system, Cinematch. Of the 8.5 billion possible customer–film pairs, about 99% are blank. Two months in, a contestant posting as Simon Funk publishes a blog post describing a startlingly simple idea. Suppose each film can be described by a few hidden numbers — how much action, how much romance, how serious — and each customer by how much they like each of those. Then the whole ratings matrix is approximately the product of a tall thin matrix and a short wide one: a low-rank matrix. Fit the factors to the ratings you have, and the product predicts the ones you do not. Three years and thousands of teams later, on 21 September 2009, BellKor's Pragmatic Chaos — a merger of teams from AT&T Labs, Austria and Canada — is awarded the prize. Their blend combines hundreds of models; at its core is matrix factorisation. The mathematics underneath is the most informative decomposition in linear algebra.
💡The SVD
Every m×nm \times n matrix — any shape, any rank — can be written A=UΣVT=∑i=1rσi uiviT,A = U\Sigma V^{\mathsf T} = \sum_{i=1}^{r} \sigma_i\,\mathbf u_i\mathbf v_i^{\mathsf T}, with UU (m×mm \times m) and VV (n×nn \times n) orthogonal and Σ\Sigma diagonal with σ1≥σ2≥⋯≥σr>0\sigma_1 \ge \sigma_2 \ge \dots \ge \sigma_r > 0, the singular values; rr is the rank. Read it as a recipe: Avi=σiuiA\mathbf v_i = \sigma_i\mathbf u_i. The matrix takes the orthonormal directions vi\mathbf v_i to orthonormal directions ui\mathbf u_i, stretched by σi\sigma_i. Every matrix maps the unit sphere to an ellipsoid whose semi-axes are σiui\sigma_i\mathbf u_i.
Derivation
Where the SVD comes from
Start from the square, symmetric matrix ATAA^{\mathsf T}A (n×nn \times n). A theorem of linear algebra — the spectral theorem — says every symmetric matrix has nn orthonormal eigenvectors v1,…,vn\mathbf v_1, \dots, \mathbf v_n with real eigenvalues λi\lambda_i (Unit 6 shows how they are computed). These eigenvalues cannot be negative:
λi=viTATA vi=∥Avi∥2≥0.\lambda_i = \mathbf v_i^{\mathsf T}A^{\mathsf T}A\,\mathbf v_i = \|A\mathbf v_i\|^2 \ge 0.
So define σi=λi\sigma_i = \sqrt{\lambda_i}, ordered from largest to smallest, and for each σi>0\sigma_i > 0 the vector ui=Avi/σi\mathbf u_i = A\mathbf v_i/\sigma_i. These u\mathbf us are orthonormal too:
uiTuj=viTATA vjσiσj=λj viTvjσiσj={1i=j0i≠j.\mathbf u_i^{\mathsf T}\mathbf u_j = \frac{\mathbf v_i^{\mathsf T}A^{\mathsf T}A\,\mathbf v_j}{\sigma_i\sigma_j} = \frac{\lambda_j\,\mathbf v_i^{\mathsf T}\mathbf v_j}{\sigma_i\sigma_j} = \begin{cases} 1 & i = j \\ 0 & i \ne j \end{cases}.
By construction Avi=σiuiA\mathbf v_i = \sigma_i\mathbf u_i for every ii (for σi=0\sigma_i = 0 both sides are zero). Put the v\mathbf vs in the columns of VV, the u\mathbf us in UU (completed with further orthonormal columns if m>nm > n), and the σ\sigmas on the diagonal of Σ\Sigma: AV=UΣAV = U\Sigma. Multiply on the right by VT=V−1V^{\mathsf T} = V^{-1}:
A=UΣVT.A = U\Sigma V^{\mathsf T}.
A consequence for lesson 5-1. ATA=VΣTUTUΣVT=V Σ2 VTA^{\mathsf T}A = V\Sigma^{\mathsf T}U^{\mathsf T}U\Sigma V^{\mathsf T} = V\,\Sigma^2\,V^{\mathsf T}: the singular values of ATAA^{\mathsf T}A are σi2\sigma_i^2, so κ(ATA)=σ12/σn2=κ(A)2\kappa(A^{\mathsf T}A) = \sigma_1^2/\sigma_n^2 = \kappa(A)^2. (Libraries do not compute the SVD this way, for exactly that reason; Golub and Kahan's algorithm works on AA directly.)
Almost everything one wants to know about a matrix can be read off its SVD:
  • Norm and condition: ∥A∥2=σ1\|A\|_2 = \sigma_1 and κ2(A)=σ1/σn\kappa_2(A) = \sigma_1/\sigma_n — the definitions of lesson 3-2 made concrete.
  • Rank: the number of nonzero singular values. In floating point nothing is exactly zero, so the numerical rank counts σi\sigma_i above a tolerance such as max⁡(m,n) u σ1\max(m,n)\,u\,\sigma_1.
  • Null space: the vi\mathbf v_i with σi=0\sigma_i = 0 — input directions the matrix destroys. Range: the ui\mathbf u_i with σi>0\sigma_i > 0.
  • Pseudoinverse: A+=∑σi>0σi−1viuiTA^+ = \sum_{\sigma_i>0} \sigma_i^{-1}\mathbf v_i\mathbf u_i^{\mathsf T}. x+=A+b\mathbf x^+ = A^+\mathbf b is the least-squares solution of smallest norm — well defined even when the columns are dependent.
Derivation
The norm, and least squares, in singular coordinates
Measure every input in the basis of the vi\mathbf v_i and every output in the basis of the ui\mathbf u_i. Orthogonal matrices preserve length (lesson 5-2), so nothing is lost by doing so — and in these coordinates AA is just Σ\Sigma.
The 2-norm. Write x=∑civi\mathbf x = \sum c_i\mathbf v_i, so ∥x∥2=∑ci2\|\mathbf x\|^2 = \sum c_i^2. Then Ax=∑ciσiuiA\mathbf x = \sum c_i\sigma_i\mathbf u_i and
∥Ax∥2=∑ici2σi2≤σ12∑ici2=σ12∥x∥2,\|A\mathbf x\|^2 = \sum_i c_i^2\sigma_i^2 \le \sigma_1^2\sum_i c_i^2 = \sigma_1^2\|\mathbf x\|^2,
with equality for x=v1\mathbf x = \mathbf v_1. So ∥A∥2=σ1\|A\|_2 = \sigma_1. For a square invertible AA the same argument applied to A−1=VΣ−1UTA^{-1} = V\Sigma^{-1}U^{\mathsf T} gives ∥A−1∥2=1/σn\|A^{-1}\|_2 = 1/\sigma_n, hence κ2=σ1/σn\kappa_2 = \sigma_1/\sigma_n.
Least squares. With y=VTx\mathbf y = V^{\mathsf T}\mathbf x and βi=uiTb\beta_i = \mathbf u_i^{\mathsf T}\mathbf b,
∥b−Ax∥2=∥UTb−Σy∥2=∑σi>0(βi−σiyi)2+∑other iβi2.\|\mathbf b - A\mathbf x\|^2 = \|U^{\mathsf T}\mathbf b - \Sigma\mathbf y\|^2 = \sum_{\sigma_i > 0}(\beta_i - \sigma_iy_i)^2 + \sum_{\text{other } i}\beta_i^2.
Each term of the first sum is made zero by yi=βi/σiy_i = \beta_i/\sigma_i; the second sum cannot be changed at all — it is the part of b\mathbf b outside the range of AA. If some σi=0\sigma_i = 0, the matching yiy_i does not appear in the residual and can be anything: the solution is not unique. Choosing those yi=0y_i = 0 gives the shortest x\mathbf x, because ∥x∥=∥y∥\|\mathbf x\| = \|\mathbf y\|. The result is the pseudoinverse solution
x+=∑σi>0uiTbσi vi.\mathbf x^+ = \sum_{\sigma_i > 0}\frac{\mathbf u_i^{\mathsf T}\mathbf b}{\sigma_i}\,\mathbf v_i.
📐A rank-deficient matrix
A=[123202011112]A = \begin{bmatrix} 1&2&3 \\ 2&0&2 \\ 0&1&1 \\ 1&1&2 \end{bmatrix} has third column = first + second. Its singular values are 5.1961525.196152, 1.7320511.732051 and ∼10−17\sim 10^{-17}: numerical rank 2. The null-space direction is v3∝(1,1,−1)\mathbf v_3 \propto (1, 1, -1) — exactly the dependency. For b=(1,2,3,4)\mathbf b = (1, 2, 3, 4), QR would break down on the zero pivot; the pseudoinverse gives the minimum-norm solution (1/3, 1/3, 2/3)(1/3,\ 1/3,\ 2/3), and says precisely which combination of unknowns the data cannot determine.
💡The best low-rank approximation
Truncate the sum after kk terms: Ak=∑i=1kσiuiviTA_k = \sum_{i=1}^{k}\sigma_i\mathbf u_i\mathbf v_i^{\mathsf T}. The error is the rest of the sum, A−Ak=∑i>kσiuiviTA - A_k = \sum_{i>k}\sigma_i\mathbf u_i\mathbf v_i^{\mathsf T} — itself an SVD, whose largest singular value is σk+1\sigma_{k+1}. The Eckart–Young theorem (1936, quoted here without its proof) adds that no matrix of rank kk does better: ∥A−Ak∥2=σk+1,∥A−Ak∥F=σk+12+⋯+σr2.\|A - A_k\|_2 = \sigma_{k+1}, \qquad \|A - A_k\|_F = \sqrt{\sigma_{k+1}^2 + \dots + \sigma_r^2}. Storage drops from mnmn numbers to k(m+n+1)k(m + n + 1). When the singular values fall quickly — as they do for smooth data, images, and anything governed by a few underlying factors — a small kk captures almost everything.
📐Hidden tastes in a ratings matrix
Simulate 300 users rating 120 films, where tastes come from three hidden factors, plus rating noise. The singular values of the observed matrix are 219.8, 181.5, 170.8219.8,\ 181.5,\ 170.8 — then a cliff to 13.6, 13.6, 13.3, …13.6,\ 13.6,\ 13.3,\ \dots: three factors, and noise. The rank-3 approximation reproduces the true tastes with 5.3% error using 3.5% of the storage. Keeping more singular values makes the fit to the observed matrix better but the estimate of the true tastes worse — rank 10 gives 11.8% error — because the extra layers model the noise. Choosing kk is choosing what to believe.
The combination of some data and an aching desire for an answer does not ensure that a reasonable answer can be extracted from a given body of data.
– John Tukey, Sunset salvo, The American Statistician (1986)
Ill-posed problems and regularisation. In Ax=bA\mathbf x = \mathbf b, the SVD writes the solution as x=∑(uiTb/σi)vi\mathbf x = \sum (\mathbf u_i^{\mathsf T}\mathbf b/\sigma_i)\mathbf v_i. Noise in b\mathbf b has components along every ui\mathbf u_i, and dividing by tiny σi\sigma_i multiplies them enormously. Deblurring a signal, reconstructing a load from measured strains, inverting seismic or gravity data — all have singular values that decay towards zero, and the naive solution is pure amplified noise. Two cures:
  • Truncated SVD: keep only the terms with σi\sigma_i above a threshold.
  • Tikhonov regularisation (Andrey Tikhonov, 1963): minimise ∥Ax−b∥2+λ2∥x∥2\|A\mathbf x - \mathbf b\|^2 + \lambda^2\|\mathbf x\|^2, which replaces 1/σi1/\sigma_i by the filter σi/(σi2+λ2)\sigma_i/(\sigma_i^2 + \lambda^2) — gently suppressing small singular values rather than cutting them.
For a blurred step-and-bump signal with 10−310^{-3} noise, the naive solve is off by a factor of more than 10310^3; truncated SVD with 40 terms gives 15% error and Tikhonov with λ=10−2\lambda = 10^{-2} gives 14%. Too little regularisation lets noise through, too much blurs the answer.
Derivation
Where Tikhonov's filter factors come from
Minimise ∥Ax−b∥2+λ2∥x∥2\|A\mathbf x - \mathbf b\|^2 + \lambda^2\|\mathbf x\|^2. The calculus of lesson 5-1 applies with one extra term, whose gradient is 2λ2x2\lambda^2\mathbf x:
(ATA+λ2I) x=ATb.(A^{\mathsf T}A + \lambda^2I)\,\mathbf x = A^{\mathsf T}\mathbf b.
Insert the SVD: ATA=VΣ2VTA^{\mathsf T}A = V\Sigma^2V^{\mathsf T} and ATb=VΣUTbA^{\mathsf T}\mathbf b = V\Sigma U^{\mathsf T}\mathbf b. In the coordinates y=VTx\mathbf y = V^{\mathsf T}\mathbf x the system is diagonal, one equation per singular value:
(σi2+λ2) yi=σi uiTb⟹yi=σiσi2+λ2 uiTb.(\sigma_i^2 + \lambda^2)\,y_i = \sigma_i\,\mathbf u_i^{\mathsf T}\mathbf b \quad\Longrightarrow\quad y_i = \frac{\sigma_i}{\sigma_i^2 + \lambda^2}\,\mathbf u_i^{\mathsf T}\mathbf b.
Compare with the unregularised yi=uiTb/σiy_i = \mathbf u_i^{\mathsf T}\mathbf b/\sigma_i. When σi≫λ\sigma_i \gg \lambda the factor is almost 1/σi1/\sigma_i: nothing changes. When σi≪λ\sigma_i \ll \lambda it is about σi/λ2\sigma_i/\lambda^2, which goes to zero instead of blowing up. The single parameter λ\lambda marks the singular value below which the data are no longer trusted.
🕰From Beltrami to Golub–Kahan
Eugenio Beltrami (1873) and Camille Jordan (1874) discovered the decomposition for square matrices while studying bilinear forms; Eckart and Young (1936) proved the approximation theorem for rectangular ones in a psychometrics journal. It was a theoretical tool until 1965, when Gene Golub and William Kahan published a stable way to compute it — bidiagonalisation by Householder reflections followed by an implicit QR iteration — and Golub and Christian Reinsch turned it into the algorithm (1970) that LAPACK still descends from.
💡Where you will use this
Attitude Determination: the optimal rotation in Wahba's problem comes straight from the SVD of a 3 × 3 matrix. Structural health monitoring and Vibration: proper orthogonal decomposition — the SVD of a matrix of sensor time histories — extracts dominant mode shapes. Optimal Guidance and Control: model reduction (balanced truncation) keeps the large Hankel singular values. Data-driven engineering: principal component analysis, recommendation systems, image compression and every 'reduced-order model'.
📚Further reading
  • Lloyd N. Trefethen and David Bau, Numerical Linear Algebra (SIAM, 1997), Lectures 4–5, 'The Singular Value Decomposition' and 'More on the SVD', and Lecture 31, 'Computing the SVD'.
  • Gilbert Strang, Introduction to Linear Algebra (5th ed., Wellesley–Cambridge, 2016), chapter 7, 'The Singular Value Decomposition' — the geometry, gently.
  • Per Christian Hansen, Discrete Inverse Problems: Insight and Algorithms (SIAM, 2010) — truncated SVD and Tikhonov regularisation for real ill-posed problems.
↓But there is a problem
Everything in this unit so far has assumed the model is linear in its unknowns. Many of the most important measurement models are not. A GNSS receiver measures the time a signal takes to arrive from each satellite. Multiplied by the speed of light, that is the distance ∣si−r∣|\mathbf s_i - \mathbf r| — plus an unknown receiver clock error. The receiver position r\mathbf r sits inside a square root, and the clock error is a fourth unknown. In 1944, at the US Army's Frankford Arsenal, a mathematician named Kenneth Levenberg wrote a short paper on fitting such models. Nineteen years later a statistician at DuPont, Donald Marquardt, rediscovered and improved his idea. Every curve-fitting tool you will ever use runs their algorithm.