Appearance
4.6 — Singular Value Decomposition and Principal Components
Eigenvalues have two limitations. They need a square matrix, and even then a matrix may not have enough eigenvectors to be useful. Most real data matrices are not square — a thousand customers by fifty attributes, or a million pixels by three colour channels.
The singular value decomposition fixes both. It works for every matrix, of every shape, real or complex, always. It is often called the most useful matrix factorisation there is, and once you see what it says geometrically it is hard to disagree.
1. The statement
Every m \times n matrix can be written
A = U\Sigma V^\mathsf{T}
where
- V is n \times n and orthogonal — its columns are perpendicular unit vectors, so it is a rotation (possibly with a reflection).
- \Sigma is m \times n, zero everywhere except the diagonal, holding the singular values \sigma_1 \ge \sigma_2 \ge \cdots \ge 0.
- U is m \times m and orthogonal — another rotation.
And here is what it means geometrically. Every linear transformation, however complicated, is a rotation, then a stretch along the (new) axes, then another rotation. That is all any matrix can ever do.
Reading right to left, as Chapter 4.2 requires:
- V^\mathsf{T} rotates the input so that the important directions line up with the axes.
- \Sigma stretches each axis by its singular value. Pure scaling, one number per axis.
- U rotates the result into its final orientation.
A test that captures it: apply any matrix to a circle and you always get an ellipse. The singular values are the ellipse's semi-axes. If one singular value is zero, the ellipse is flat — the matrix squashed a dimension, which is the singular matrix of Chapter 4.4 seen from another angle.
Relationship to eigenvalues. The singular values of A are the square roots of the eigenvalues of A^\mathsf{T}A, which is always square and symmetric and therefore always well behaved by Chapter 4.5's spectral theorem. That is how the SVD sidesteps every difficulty eigenvalues have.
2. The rank-one decomposition, which is where the power is
Multiply out the factorisation and it becomes a sum:
A = \sigma_1 \mathbf{u}_1\mathbf{v}_1^\mathsf{T} + \sigma_2\mathbf{u}_2\mathbf{v}_2^\mathsf{T} + \cdots + \sigma_r\mathbf{u}_r\mathbf{v}_r^\mathsf{T}
Each term is a full-size matrix built from one column of U and one column of V, weighted by a singular value. Each is as simple as a matrix can be — rank one, in the sense of Chapter 4.2.
The singular values are ordered largest first, so the first terms carry the most. Cut the sum off after k terms and you get the best possible rank-k approximation of A — best in the precise sense of minimising the total squared error. That result is the Eckart–Young theorem, and it is what makes the SVD an engineering tool rather than a curiosity.
In plain terms: the SVD sorts a matrix's content by importance and lets you keep the top slice.
Image compression, worked
A greyscale photograph is a matrix of brightness values. Take one that is 1000 \times 1000 — a million numbers.
Compute the SVD and keep only the first 50 terms. You now store 50 columns of U (50,000 numbers), 50 singular values, and 50 columns of V (50,000 numbers). About 100,050 numbers instead of a million: a tenfold reduction.
And the picture still looks like the picture. The large singular values captured the broad structure — the shapes, the shading, the outlines — while the discarded small ones held fine texture and noise.
Push it further and the image gets progressively blurrier and blockier, degrading gracefully rather than falling apart. This is not how JPEG actually works (JPEG uses the discrete cosine transform of Chapter 9.3), but the principle is identical: transform into a basis where importance is concentrated, then throw away the unimportant part.
3. Principal component analysis
PCA is the SVD applied to data, and it is one of the most-used techniques in all of statistics and machine learning.
The problem. You have data with many columns — 50 measurements per patient, 200 features per customer. You cannot plot 50 dimensions, many of the columns are correlated and therefore redundant, and everything downstream is slower and less stable than it needs to be.
The question PCA asks. Along which direction does the data vary the most? Then, among directions perpendicular to that one, which varies most next? And so on.
Those directions are the principal components, and they are exactly the eigenvectors of the data's covariance matrix, which is the same as the right singular vectors of the centred data matrix.
The procedure, in five steps.
- Centre the data — subtract each column's mean, so the cloud sits around the origin. Skipping this is the most common mistake, and it makes the first component point at the mean rather than along the spread.
- Optionally standardise — divide each column by its standard deviation. Do this when columns are in different units, otherwise a column measured in millimetres will dominate one measured in metres purely because its numbers are bigger.
- Compute the SVD of the centred matrix.
- The right singular vectors are the components, and the squared singular values are proportional to how much variance each one explains.
- Keep enough components to cover the variance you want — often 90 or 95 percent — and project the data onto them.
Reading the result. If two components explain 92% of the variance in a 40-column dataset, then your data was really two-dimensional wearing a forty-dimensional costume. The other 38 directions were noise and redundancy.
Worked intuition. Measure a hundred people's height, arm span, leg length, shoe size, weight and collar size. Six columns. But these all move together — tall people are large in all of them. PCA would find that one component ("overall size") explains perhaps 85% of the variation, and a second ("build: lanky versus stocky") most of the rest. Six measurements, two real dimensions.
Honest limitations, because PCA is routinely over-applied:
- The components are combinations of the original columns, so they are often uninterpretable. "0.3 × height + 0.5 × weight − 0.2 × age" is a direction, not a meaning.
- It only finds linear structure. Data lying on a curved surface will defeat it.
- It maximises variance, and variance is not always the same as usefulness. A low-variance direction can be the one that separates two classes.
- It is sensitive to scaling, which is why step 2 matters.
4. Other things the SVD gives you for free
The rank is the number of nonzero singular values. On a computer nothing is exactly zero, so in practice you count the singular values above a small threshold — this is the numerically reliable way to determine rank, far better than Chapter 4.3's elimination.
The condition number is \frac{\sigma_{\max}}{\sigma_{\min}}. Chapter 4.3 met ill-conditioning as an ugly surprise; here it is a number you can compute. A large ratio means the matrix stretches enormously in one direction and barely at all in another, so inverting it amplifies error. Anything above about 10^{8} in single precision, or 10^{16} in double, means your answer has no reliable digits left.
Least squares, robustly. The pseudoinverse built from the SVD solves A\mathbf{x} = \mathbf{b} in the least-squares sense of Chapter 4.3, and unlike the normal equations it stays stable when the matrix is nearly rank-deficient. It also gives the smallest solution when there are infinitely many.
Noise removal. Noise tends to spread evenly across all directions, so it lands mostly in the small singular values. Zeroing those and reconstructing removes noise while keeping the structure. This is used on seismic data, astronomical images and medical scans.
5. Recommendation systems, the famous application
Build a matrix with users as rows and films as columns, holding ratings. It is enormous and almost entirely empty — a typical user has rated a few dozen of thousands of films.
Take a low-rank SVD approximation. Something remarkable happens: the factorisation is forced to explain millions of ratings with only, say, 50 numbers per user and 50 per film. It cannot memorise; it has to find structure.
The structure it finds is a set of hidden dimensions — nobody labelled them, but they turn out to correspond to things like "how much action", "how old", "how mainstream". Each user becomes a point in that 50-dimensional taste space and each film a point in the same space, and a predicted rating is the dot product of the two, exactly as in Chapter 4.1.
And the empty cells get filled in. Multiply the factors back together and every user-film pair has a number, including all the ones nobody ever rated. Those are the recommendations.
The Netflix Prize (2006–2009) offered a million dollars for a 10% improvement in rating prediction, and matrix factorisation of this kind was the core of every leading entry. Modern systems use neural networks, but they learn embeddings that play the same role — Volume I, 11.20 covers the engineering.
Every formula above, built from scratch
None of the results in this chapter are worth memorising, because each one can be rebuilt in under a minute from something simpler. What follows is that rebuilding, one result at a time, so the formula and the reason for it sit on the same page as the explanation that needed them.
Orthogonality and Gram–Schmidt
A set of vectors is orthonormal when each has length 1 and every pair has dot product 0. Such a set is the best possible coordinate system, because finding a vector's coordinates in it costs only dot products:
\mathbf{v} = (\mathbf{v}\cdot\mathbf{q}_1)\mathbf{q}_1 + (\mathbf{v}\cdot\mathbf{q}_2)\mathbf{q}_2 + \cdots
Gram–Schmidt turns any independent set into an orthonormal one, by repeatedly subtracting off the parts that are already accounted for:
\mathbf{w}_k = \mathbf{v}_k - \sum_{j<k}\text{proj}_{\mathbf{q}_j}\mathbf{v}_k, \qquad \mathbf{q}_k = \frac{\mathbf{w}_k}{\|\mathbf{w}_k\|}
Read it as a procedure: take the next vector, remove everything pointing along the directions you have already fixed, and what is left must be perpendicular to all of them. Normalise it and move on.
Orthogonal matrices, whose columns are orthonormal, satisfy
Q^TQ = I, \qquad Q^{-1} = Q^T, \qquad \|Q\mathbf{x}\| = \|\mathbf{x}\|, \qquad \det Q = \pm1
They are exactly the rotations and reflections: transformations that move things without stretching them. The inverse being the transpose is why they are computationally precious — no division, no loss of accuracy.
Least squares
A^TA\hat{\mathbf{x}} = A^T\mathbf{b} \quad \Rightarrow \quad \hat{\mathbf{x}} = (A^TA)^{-1}A^T\mathbf{b}
The problem. You have more equations than unknowns — 500 data points and a straight line with 2 parameters — so A\mathbf{x}=\mathbf{b} has no solution at all. \mathbf{b} simply is not in the space that A can reach.
The fix, and why it is this formula. If you cannot hit \mathbf{b}, hit the closest reachable point instead. The closest point is the projection of \mathbf{b} onto the space A can reach, and "closest" means the error vector \mathbf{b}-A\hat{\mathbf{x}} is perpendicular to everything reachable. Perpendicular to every column of A means
A^T(\mathbf{b}-A\hat{\mathbf{x}}) = \mathbf{0}
Multiply out and rearrange:
A^T\mathbf{b} = A^TA\hat{\mathbf{x}}
That is the whole derivation, and it is the projection formula of §1 wearing matrix clothing. Every line of best fit, every linear regression in every statistics package, is this equation being solved.
For a straight line y = mx+c through points (x_i,y_i), it works out to:
m = \frac{n\sum x_iy_i - \sum x_i\sum y_i}{n\sum x_i^2 - \left(\sum x_i\right)^2}, \qquad c = \frac{\sum y_i - m\sum x_i}{n}
Singular value decomposition
A = U\Sigma V^T
Every matrix, of any shape, splits into three: V^T is a rotation, \Sigma is a diagonal stretch, and U is another rotation. Both U and V are orthogonal, and \Sigma holds the singular values \sigma_1 \ge \sigma_2 \ge \cdots \ge 0 down its diagonal.
Where the pieces come from. A^TA is symmetric, so by §6 it has real eigenvalues and perpendicular eigenvectors. Those eigenvectors are the columns of V, and
\sigma_i = \sqrt{\lambda_i(A^TA)}
The columns of U are then \mathbf{u}_i = \frac{1}{\sigma_i}A\mathbf{v}_i.
What it says in plain terms. Any linear transformation, however complicated it looks, is a rotation followed by a stretch along perpendicular axes followed by another rotation. There is nothing else it can do.
Low-rank approximation is the payoff:
A \approx \sum_{i=1}^{k}\sigma_i\,\mathbf{u}_i\mathbf{v}_i^T
Keep the largest k singular values and throw the rest away. This is provably the best possible approximation of rank k, and it is how image compression, noise removal, and latent-factor recommendation all work.
Principal component analysis is the SVD of your data after subtracting the mean of each column. The principal directions are the columns of V, and the variance captured by direction i is \frac{\sigma_i^2}{\sum\sigma_j^2}.
6. Where this shows up in your life
Every recommendation you receive, from films to products to the next video.
Every search engine's semantic matching. Latent semantic analysis is the SVD of a document-by-word matrix, and it lets "car" and "automobile" land near each other without anyone writing a synonym list.
Every face-recognition system's ancestry. Eigenfaces — principal components of face images — were the first method that worked, and the modern deep-learning versions are doing a nonlinear version of the same compression.
Every genomics and neuroscience study. Thousands of genes or channels reduced to a handful of components before anything else is attempted.
Every noise-reduction and background-subtraction feature in image and audio software.
Every "explained variance" figure in a research paper. That is PCA reporting how much of the data its top components captured.
Part 4 is finished. We can now describe space, transform it, solve large systems within it, and find the directions that matter. What we still cannot do is handle change — what happens to a quantity in the instant, and how small changes accumulate into large ones. That is calculus, and it is the largest and most consequential idea in this volume.
More places these turn up
Every rotation on your phone's screen is a 2\times2 matrix; every 3D game frame multiplies millions of vertices by a 4\times4 one. The dot product is how a search engine measures whether your query points in the same direction as a document. The least-squares formula of §8 is behind every trend line in every spreadsheet. The SVD compresses the images you send and powers the "people also liked" list. And the eigenvalue equation, in a form Google called PageRank, is what sorted the entire web — the ranking of every page is the eigenvector of a very large matrix, computed by the repeated-multiplication trick of §6.
Next: 4.P — Worked Problems works through vectors, systems, determinants, eigenvalues and a real least-squares fit, step by step.