SVD and PCA: The Shape of Data

Discover: a few numbers for a face

In 1991, Matthew Turk and Alex Pentland at the MIT Media Lab built a system that could recognise a face in about 0.4 seconds. That was a hard problem. A $256 \times 256$ greyscale photo is a list of 65,536 numbers: one point in 65,536-dimensional space (Chamber 1's data as vectors). Comparing faces pixel by pixel is slow, and it's thrown off by small changes in lighting or position. Turk and Pentland described each face with a handful of numbers instead, computed like this:

Read it aloud: “omega k equals u k transpose, times gamma minus psi”. $\Gamma$ (capital gamma) is a new face image, unrolled into one long column of pixels. $\Psi$ (capital psi) is the average face of the training set, so $\Gamma - \Psi$ is “how this face differs from the average face”. Each $\mathbf{u}_k$ is an eigenface, a ghostly face-shaped pattern, and the weight $\omega_k$ (omega k) says how much of pattern $k$ this face contains. In many of their tests, trained on $M = 16$ images, they kept just seven eigenfaces ($M' = 7$, “M prime”): seven numbers per face instead of tens of thousands of pixels.

That raises three questions. Where do the patterns $\mathbf{u}_k$ come from? Why those patterns and not others? And how can so few numbers carry so much? The answers are the singular value decomposition, which says that every matrix is built from a few strong layers and many weak ones, and principal component analysis, which is the same idea applied to a cloud of data.

Learn: rotate, stretch, rotate

In Chamber 4 you watched matrices move the plane. Apply a $2 \times 2$ matrix to the unit circle and, unless it flattens the plane, you get an ellipse. An ellipse has a longest axis and a shortest axis, at right angles to each other. Here's the surprise: there is always a pair of perpendicular unit vectors $\mathbf{v}_1, \mathbf{v}_2$ that the matrix sends exactly onto those axes. It sends $\mathbf{v}_1$ to the long semi-axis, of length $\sigma_1$ (“sigma one”) in the direction of a unit vector $\mathbf{u}_1$, and $\mathbf{v}_2$ to the short one, of length $\sigma_2$ along $\mathbf{u}_2$:

Put the $\mathbf{v}$'s in the columns of a matrix $V$, the $\mathbf{u}$'s in the columns of $U$, and the $\sigma$'s on the diagonal of $\Sigma$ (capital sigma). The two equations become one, $AV = U\Sigma$. Since $V$ has orthonormal columns, $V^\top V = I$ and $V^{-1} = V^\top$ (Chamber 7), so multiplying on the right by $V^\top$ gives

read “A equals U sigma V transpose”. This is the singular value decomposition (SVD), and the remarkable fact is that every real matrix has one: square or rectangular, invertible or not.

Read the factorisation from right to left, as three steps applied to a vector $\mathbf{x}$:

1. Rotate. $V^\top\mathbf{x}$ lists the dot products $\mathbf{v}_i^\top\mathbf{x}$: the coordinates of $\mathbf{x}$ in the basis of $\mathbf{v}$'s (Chamber 3). Because $V$ is orthogonal, this step is a rotation, possibly with a reflection. No length changes.
2. Stretch. $\Sigma$ multiplies coordinate $i$ by $\sigma_i$. This is the only step that changes lengths. When $m \ne n$ it also pads with zeros or drops coordinates, which is how you get from $\mathbb{R}^n$ to $\mathbb{R}^m$.
3. Rotate. $U$ turns the stretched axes to point along $\mathbf{u}_1, \mathbf{u}_2, \ldots$ in the output space.

So every linear map, however tangled it looks, is a rotation, a stretch along perpendicular axes, and another rotation. A sphere always becomes an ellipsoid (possibly flattened), with semi-axes $\sigma_1\mathbf{u}_1, \sigma_2\mathbf{u}_2, \ldots$. In particular, $\sigma_1$ is the biggest stretch $A$ applies to any unit vector. (Write a unit vector as $\mathbf{x} = \sum_i c_i\mathbf{v}_i$ with $\sum_i c_i^2 = 1$. Then $\|A\mathbf{x}\|^2 = \sum_i \sigma_i^2 c_i^2$, which is at most $\sigma_1^2$.)

A worked example. Take $A = \begin{pmatrix} 3 & 0 \\ 4 & 5 \end{pmatrix}$. Its SVD is $A = U\Sigma V^\top$ with

Multiply the three out and you get $A$ back.

So $\mathbf{v}_1 = \tfrac{1}{\sqrt{2}}(1, 1)$, $\mathbf{v}_2 = \tfrac{1}{\sqrt{2}}(-1, 1)$, $\mathbf{u}_1 = \tfrac{1}{\sqrt{10}}(1, 3)$ and $\mathbf{u}_2 = \tfrac{1}{\sqrt{10}}(-3, 1)$. Check one column by hand: $A\mathbf{v}_1 = \tfrac{1}{\sqrt{2}}(3, 9)$, and $\sigma_1\mathbf{u}_1 = \tfrac{3\sqrt{5}}{\sqrt{10}}(1, 3) = \tfrac{3}{\sqrt{2}}(1, 3)$. The same vector.

This $A$ turns the unit circle into an ellipse with semi-axes $3\sqrt{5} \approx 6.71$ and $\sqrt{5} \approx 2.24$. The circle's area $\pi$ becomes $\pi\sigma_1\sigma_2 = 15\pi$, and indeed $\det A = 15$ (Chamber 5). For any square matrix, $|\det A| = \sigma_1\sigma_2\cdots\sigma_n$: the rotations don't change volume, only the stretch does.

Learn: where the SVD comes from

Why should every matrix have such perpendicular directions? The derivation is short, and it rests on the matrix $A^\top A$ (“A transpose A”), an $n \times n$ matrix that you met in the normal equations of Chamber 7.

Theorem. Every real $m \times n$ matrix $A$ has a singular value decomposition.

Proof. We build $V$, then $\Sigma$, then $U$.

1. 1

   A transpose A is symmetric

   $(A^\top A)^\top = A^\top (A^\top)^\top = A^\top A$. By the spectral theorem (the last section of this chamber), a symmetric matrix has real eigenvalues and an orthonormal basis of eigenvectors. Call the eigenvectors $\mathbf{v}_1, \ldots, \mathbf{v}_n$, with eigenvalues ordered $\lambda_1 \ge \cdots \ge \lambda_n$.

2. 2

   Its eigenvalues are never negative

   Since $\mathbf{v}_i$ is a unit vector,

   $$\begin{aligned} \lambda_i = \lambda_i\,\mathbf{v}_i^\top\mathbf{v}_i &= \mathbf{v}_i^\top A^\top A\,\mathbf{v}_i \\ &= (A\mathbf{v}_i)^\top(A\mathbf{v}_i) \\ &= \|A\mathbf{v}_i\|^2 \ge 0. \end{aligned}$$

   So we can take square roots. Define $\sigma_i = \sqrt{\lambda_i} = \|A\mathbf{v}_i\|$.

3. 3

   Divide to get the u's

   Let $r$ be the number of non-zero $\sigma_i$. For $i \le r$, define $\mathbf{u}_i = A\mathbf{v}_i / \sigma_i$. These are orthonormal:

   $$\mathbf{u}_i^\top\mathbf{u}_j = \frac{\mathbf{v}_i^\top A^\top A\,\mathbf{v}_j}{\sigma_i\sigma_j} = \frac{\lambda_j\,\mathbf{v}_i^\top\mathbf{v}_j}{\sigma_i\sigma_j} = \delta_{ij},$$

   the Kronecker delta (1 if $i = j$, 0 otherwise), because $\mathbf{v}_i^\top\mathbf{v}_j$ is 0 when $i \ne j$, and when $i = j$ the fraction is $\lambda_i / \sigma_i^2 = 1$.

4. 4

   Fill in the rest

   For $i > r$, $\|A\mathbf{v}_i\| = \sigma_i = 0$, so $A\mathbf{v}_i = \mathbf{0}$. Extend $\mathbf{u}_1, \ldots, \mathbf{u}_r$ to an orthonormal basis $\mathbf{u}_1, \ldots, \mathbf{u}_m$ of $\mathbb{R}^m$ (Gram–Schmidt, Chamber 7). Now $A\mathbf{v}_i = \sigma_i\mathbf{u}_i$ for every $i$ up to $\min(m, n)$, and $A\mathbf{v}_i = \mathbf{0}$ for any $i$ beyond $m$. Column by column, that is exactly $AV = U\Sigma$.

5. 5

   Multiply by V transpose

   $V$ is orthogonal, so $VV^\top = I$, and $A = AVV^\top = U\Sigma V^\top$. $\blacksquare$

The proof tells you more than existence:

• The singular values are $\sigma_i = \sqrt{\lambda_i(A^\top A)}$, and the $\mathbf{v}_i$ are eigenvectors of $A^\top A$. By the same argument with $AA^\top = U\Sigma\Sigma^\top U^\top$, the $\mathbf{u}_i$ are eigenvectors of $AA^\top$, which has the same non-zero eigenvalues.
• The rank of $A$ is $r$, the number of non-zero singular values. $\mathbf{u}_1, \ldots, \mathbf{u}_r$ are an orthonormal basis for the column space of $A$, and $\mathbf{v}_{r+1}, \ldots, \mathbf{v}_n$ for its null space (Chamber 6). The SVD hands you orthonormal bases for the subspaces Chamber 6 introduced.

For the worked example, $A^\top A = \begin{pmatrix} 25 & 20 \\ 20 & 25 \end{pmatrix}$ has eigenvalues 45 and 5, with eigenvectors along $(1, 1)$ and $(-1, 1)$. Their square roots, $\sqrt{45} = 3\sqrt{5}$ and $\sqrt{5}$, are the singular values we found.

Learn: a sum of layers, and the best low-rank copy

Chamber 4 gave four views of matrix multiplication. The fourth, columns times rows, turns $U\Sigma V^\top$ into a sum:

Each $\mathbf{u}_i\mathbf{v}_i^\top$ is an outer product, a rank-one matrix (Chamber 6): every row is a multiple of the same row $\mathbf{v}_i^\top$. Drawn as an image it looks like a plaid, stripes in two directions. The SVD writes any matrix as a sum of plaids, sorted from strongest to weakest. Keep the first $k$ and drop the rest, and you get the truncated SVD, “A sub k”:

$A_k$ has rank $k$. How close is it to $A$? To answer, we need the size of a matrix. The Frobenius norm unrolls $A$ into one long vector and takes its length (Chamber 2):

read “the Frobenius norm of A”. Then $\|A - B\|_F$ is the distance between two matrices, entry by entry. For images, it's the square root of the total squared pixel error.

Claim. $\|A\|_F^2 = \sigma_1^2 + \sigma_2^2 + \cdots + \sigma_r^2$.

Proof. First, multiplying by an orthogonal matrix $Q$ doesn't change the Frobenius norm. For any vector $\mathbf{b}$, $\|Q\mathbf{b}\|^2 = \mathbf{b}^\top Q^\top Q\,\mathbf{b} = \mathbf{b}^\top\mathbf{b} = \|\mathbf{b}\|^2$ (Chamber 7). The squared Frobenius norm of a matrix $B$ is the sum of its squared column lengths, so $\|QB\|_F^2 = \sum_j \|Q\mathbf{b}_j\|^2 = \sum_j \|\mathbf{b}_j\|^2 = \|B\|_F^2$. Multiplying on the right is the same argument applied to rows: transposing only rearranges entries, so $\|BQ\|_F = \|Q^\top B^\top\|_F = \|B^\top\|_F = \|B\|_F$, since $Q^\top$ is orthogonal too. Now peel the orthogonal factors off the SVD:

the last step because the only non-zero entries of $\Sigma$ are the $\sigma_i$. $\blacksquare$

The same argument measures the truncation error. $A - A_k = \sum_{i > k}\sigma_i\mathbf{u}_i\mathbf{v}_i^\top$ is again of the form $U(\text{diagonal})V^\top$, with $\sigma_1, \ldots, \sigma_k$ replaced by zeros, so

The error is exactly the size of the singular values you threw away. You can read the quality of every truncation straight off the list of singular values, the spectrum.

We won't prove Eckart–Young here, since the proof takes a page. The readings below include a proof of the spectral-norm version and an intuitive argument for the Frobenius one. But notice what it means for machine learning. Whenever a method approximates a known matrix by a low-rank one, say to compress the weights of a trained layer, the truncated SVD is the benchmark: no rank-$k$ matrix gets closer in the Frobenius norm. (LoRA, from Chamber 6, is different: it learns a low-rank update instead of approximating a known one.)

Storage. To store $A_k$ you keep $k$ columns of $U$ ($m$ numbers each), $k$ columns of $V$ ($n$ each) and $k$ singular values: $k(m + n + 1)$ numbers instead of $mn$. A $1000 \times 1000$ photo at rank 50 needs $50 \times 2001 = 100{,}050$ numbers, about a tenth of the original million. Try it on two images, each a $64 \times 64$ matrix of brightness values between 0 (black) and 1 (white).

Things to try:

1. At rank 1 the picture is a plaid: every row is a multiple of $\mathbf{v}_1^\top$. Step up one layer at a time and watch the third panel, the newest layer $\sigma_k\mathbf{u}_k\mathbf{v}_k^\top$, gold where it adds brightness and blue where it removes it. Every layer is another plaid.
2. Watch the spectrum while you go. The error readout is computed from the bars you haven't kept, exactly as the formula above says. The “energy kept” is $1$ minus the relative error squared.
3. Switch to the Greek key. At some rank the error drops to exactly zero and the spectrum falls off a cliff. Why? (Hint: how many different rows does the image have, and are any two of them multiples of each other? Think back to Chamber 6.)
4. The challenge: for either image, find the smallest rank that rebuilds it with under 10% relative error, and leave the slider there.

Learn: principal component analysis

Now switch from one matrix to a dataset: $n$ points $\mathbf{x}_1, \ldots, \mathbf{x}_n \in \mathbb{R}^d$, stored as the rows of an $n \times d$ data matrix $X$. Which directions matter?

In 1901 Karl Pearson asked for the line (and the plane) of closest fit to a cloud of points, measuring each point's distance perpendicular to the line rather than vertically, as regression does. Pythagoras turns this into a question about variance. Suppose the data are centred (their mean is $\mathbf{0}$) and the line runs along a unit vector $\mathbf{w}$. Each point splits into a part along the line, of length $\mathbf{w}^\top\mathbf{x}_i$, and a perpendicular part:

The left side doesn't depend on $\mathbf{w}$. So making the total squared distance small is the same as making $\sum_i (\mathbf{w}^\top\mathbf{x}_i)^2$ large: the best-fitting line points in the direction along which the data vary the most. That direction is the first principal component.

1. 1

   Centre

   Compute the mean $\bar{\mathbf{x}} = \frac{1}{n}\sum_i \mathbf{x}_i$ (“x bar”) and subtract it from every row: $X_c = X - \mathbf{1}\bar{\mathbf{x}}^\top$ (“X centred”). The best-fitting line passes through the mean, so this loses nothing.

2. 2

   Form the covariance matrix

   $C = \frac{1}{n}X_c^\top X_c$, a $d \times d$ matrix. Entry $C_{jk}$ is the covariance of features $j$ and $k$, and the diagonal holds the variances. Statisticians often divide by $n - 1$ instead of $n$: the eigenvectors are the same and the eigenvalues scale.

3. 3

   Find its eigenvectors

   $C$ is symmetric, so $C = Q\Lambda Q^\top$ with orthonormal eigenvectors $\mathbf{q}_1, \ldots, \mathbf{q}_d$ (the principal directions) and eigenvalues $\lambda_1 \ge \cdots \ge \lambda_d \ge 0$ (the variances along them).

4. 4

   Measure the explained variance

   Component $i$ explains the fraction $\lambda_i / (\lambda_1 + \cdots + \lambda_d)$ of the total variance. Plotting these fractions (a “scree plot”) shows how many components are worth keeping.

5. 5

   Project

   Keep the top $k$ directions, $Q_k = [\mathbf{q}_1 \cdots \mathbf{q}_k]$, and replace each point by its scores $\mathbf{z}_i = Q_k^\top(\mathbf{x}_i - \bar{\mathbf{x}}) \in \mathbb{R}^k$: $k$ numbers instead of $d$.

Why the top eigenvector? The variance of the data along a unit vector $\mathbf{w}$ is

using $(\mathbf{w}^\top\mathbf{a})^2 = \mathbf{w}^\top\mathbf{a}\mathbf{a}^\top\mathbf{w}$, and the fact that a sum of outer products of the rows is $X_c^\top X_c$ (columns times rows again).

Claim. Over all unit vectors $\mathbf{w}$, the variance $\mathbf{w}^\top C\mathbf{w}$ is largest at $\mathbf{w} = \mathbf{q}_1$, where it equals $\lambda_1$.

Proof. Write $\mathbf{w}$ in the orthonormal eigenbasis: $\mathbf{w} = a_1\mathbf{q}_1 + \cdots + a_d\mathbf{q}_d$, with $a_1^2 + \cdots + a_d^2 = \|\mathbf{w}\|^2 = 1$. Then $C\mathbf{w} = \sum_j a_j\lambda_j\mathbf{q}_j$, and because $\mathbf{q}_i^\top\mathbf{q}_j = 0$ for $i \ne j$, all the cross terms vanish:

The inequality holds because every $\lambda_j \le \lambda_1$, and $\mathbf{w} = \mathbf{q}_1$ reaches the bound, since $\mathbf{q}_1^\top C\mathbf{q}_1 = \lambda_1$. $\blacksquare$

Repeat the argument among directions perpendicular to $\mathbf{q}_1$ and you get $\mathbf{q}_2$, and so on. This is word for word the definition in the eigenfaces paper, whose Equation 1 asks for the $\mathbf{u}_k$ for which $\lambda_k = \frac{1}{M}\sum_n(\mathbf{u}_k^{T}\Phi_n)^2$ is a maximum, subject to the $\mathbf{u}$'s being orthonormal. It also tells you what PCA throws away. By Pythagoras again, the average squared distance from the points to the best $k$-dimensional subspace is $\lambda_{k+1} + \cdots + \lambda_d$: exactly the variance in the directions you dropped.

A worked example. Take four points in $\mathbb{R}^2$: $(5, 4)$, $(-1, -2)$, $(3, 0)$ and $(1, 2)$.

1. The mean is $\bar{\mathbf{x}} = (2, 1)$. The centred points are $(3, 3)$, $(-3, -3)$, $(1, -1)$ and $(-1, 1)$.
2. The covariance matrix is $C = \frac{1}{4}\begin{pmatrix} 9 + 9 + 1 + 1 & 9 + 9 - 1 - 1 \\ 9 + 9 - 1 - 1 & 9 + 9 + 1 + 1 \end{pmatrix} = \begin{pmatrix} 5 & 4 \\ 4 & 5 \end{pmatrix}$.
3. $\det(C - \lambda I) = (5 - \lambda)^2 - 16 = 0$ gives $\lambda_1 = 9$ and $\lambda_2 = 1$, with eigenvectors $\mathbf{q}_1 = \tfrac{1}{\sqrt{2}}(1, 1)$ and $\mathbf{q}_2 = \tfrac{1}{\sqrt{2}}(1, -1)$.
4. The first component explains $9 / (9 + 1) = 90\%$ of the variance.
5. The scores along $\mathbf{q}_1$ are $3\sqrt{2}$, $-3\sqrt{2}$, $0$ and $0$. Their variance is $(18 + 18 + 0 + 0)/4 = 9 = \lambda_1$, as promised. Keeping one component squashes the last two points onto the mean, and the average squared distance lost is $(0 + 0 + 2 + 2)/4 = 1 = \lambda_2$.

PCA through the SVD. Take the SVD of the centred data, $X_c = U\Sigma V^\top$. Then

That is an eigendecomposition of $C$. The principal directions are the right singular vectors of $X_c$, the variances are $\lambda_i = \sigma_i^2/n$, and the scores come for free: $Z = X_cV_k = U_k\Sigma_k$. In the worked example, $X_c$ has singular values 6 and 2, and indeed $6^2/4 = 9$ and $2^2/4 = 1$. This is how PCA is computed in practice (scikit-learn's PCA centres the data and takes an SVD), because it's more accurate than forming $C$, and because it's cheap when there are far fewer points than features. That was Turk and Pentland's situation: they found eigenvectors $\mathbf{v}_l$ of the small $M \times M$ matrix $A^{T}A$ and formed the eigenfaces as $\mathbf{u}_l = A\mathbf{v}_l$, which is step 3 of our derivation of the SVD without the division by $\sigma_l$.

Learn: the spectral theorem

You've now leaned on one fact three times: $A^\top A$ and $C$ are symmetric, so they have real eigenvalues and perpendicular eigenvectors. Time to name it.

Read it “S equals Q lambda Q transpose”. Each $\mathbf{q}_i\mathbf{q}_i^\top$ projects onto the line through $\mathbf{q}_i$ (Chamber 7), so a symmetric matrix simply stretches space by $\lambda_i$ along $n$ perpendicular axes. That's the SVD's picture with the same axes going in and coming out. For a symmetric matrix with no negative eigenvalues, such as $A^\top A$ or $C$, the eigendecomposition is an SVD: $U = V = Q$ and $\Sigma = \Lambda$. When some eigenvalues are negative, the singular values are $|\lambda_i|$ and the signs move into $U$.

The proof has three parts.

1. The eigenvalues are real. The proof needs complex numbers, and it's in the box below.
2. Eigenvectors for different eigenvalues are perpendicular. You proved this yourself in Chamber 8's Your turn: compute $\mathbf{v}_1^\top S\mathbf{v}_2$ in two ways and use $S^\top = S$.
3. Even when eigenvalues repeat, there are enough perpendicular eigenvectors to fill a basis. This needs an induction argument that we'll skip; Theorem 4.15 of Mathematics for Machine Learning states it, and linear algebra textbooks prove it.

Read beyond the course

The SVD sits at a crossroads of numerical computing, statistics and machine learning, so every field has written about it. These are free, and each one shows the ideas from a different side.

Papers and lectures

The eigenfaces paper is one of the most readable classics in machine learning, and you now know everything its central sections use. Read Calculating Eigenfaces (Equations 1–6) and then Using Eigenfaces to Classify a Face Image (Equations 7–9). One thing to watch: the paper stores each face $\Phi_n$ as a column, so its matrix $A = [\Phi_1\ \Phi_2 \cdots \Phi_M]$ is the transpose of our $X_c$. On a first read, skip the later sections on locating faces in a scene and tracking moving heads.

Eckart and Young's paper is eight pages long and written in the notation of 1936. Its abstract is worth reading for one sentence: the normal equations of least squares (Chamber 7) “cannot be immediately written down”, because the entries of a low-rank matrix aren't independent of one another. They solve it by “first expressing the matrices in a canonic form”, which is the SVD. Pearson's paper is where PCA began: fitting a line or a plane by perpendicular distance, rather than by the vertical errors of regression.

Now take apart the equation at the heart of the eigenfaces paper.

Watch

Brunton's overview follows Chapter 1 of his book, where each data point is a column (as in the eigenfaces paper), so his $U$ holds the patterns and $V$ the weights: the transpose of this chamber's rows convention. StatQuest works PCA out slowly on a small example, from finding PC1 to the percentage of variation per component and the scree plot. Watch how he scales things: he works with sums of squared distances, where this chamber divides by $n$ to get variances. The directions and the percentages come out the same.

Your turn

You've already found the smallest rank in the lab above. Now match the symbols, assemble one proof and write another, and then compress, decompose and recognise in code.