3.3 SVD, Low-Rank Approximation, and PCA

: 40 minutes

This is where the week pays off. Everything in 3.1 and 3.2 was setting up one factorization, and that factorization is what makes ridge regression legible on Nov 3 and what we run before clustering on Dec 8.

Singular Value Decomposition

The singular value decomposition (SVD) is the most widely used factorization in numerical linear algebra. It spawns applications across science and engineering, and it provides the foundation for Principal Component Analysis (PCA). PCA is in turn one technique in the much larger area of dimensionality reduction, which seeks low-dimensional approximations to very high-dimensional data.

In this chapter we describe the factorization, then two of its uses: low-rank approximation and PCA.

Definition

Suppose we are given a data set in the form of an m\times n matrix: \mathbf{X}=\begin{bmatrix}\begin{matrix}\Big|\\\mathbf{x}_1\\\Big|\end{matrix} & \begin{matrix}\Big|\\\mathbf{x}_2\\\Big|\end{matrix} &\cdots& \begin{matrix}\Big|\\\mathbf{x}_n\\\Big|\end{matrix}\end{bmatrix}. Each column \mathbf{x}_i of \mathbf{X} is a vector of \mathbb{R}^m.

A convention that flips halfway through this chapter

For the SVD itself we read \mathbf{X} column by column, as above. For PCA we will read it row by row, one observation per row, because that is how a spreadsheet and a DataFrame store data. Nothing mathematical changes; only which index you call an observation. Watch for the switch at the PCA section and re-read the shapes there.

The SVD is a matrix decomposition that exists for every matrix \mathbf{X}—square or not, singular or not: \mathbf{X}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^T, where

\mathbf{U}_{m\times m},\mathbf{V}_{n\times n} are square, orthogonal matrices¹
the columns of \mathbf{U} and \mathbf{V} are called left singular vectors and right singular vectors, respectively
\mathbf{\Sigma}_{m\times n} is rectangular diagonal: zero everywhere off the main diagonal, with non-negative entries on it
the diagonal entries of \mathbf{\Sigma} are the singular values, ordered from largest to smallest.

¹ \mathbf{U}\mathbf{U}^T=\mathbf{U}^T\mathbf{U}=\mathbf{I}.

The singular values are unique. The singular vectors are not: you may flip the sign of any \mathbf{u}_k provided you flip \mathbf{v}_k with it, and when two singular values coincide you may rotate freely within that block. This is why two SVD implementations can hand you visibly different \mathbf{U} and \mathbf{V} and both be right, and it is why comparisons in the exercises below use np.abs.

Singular Values and Rank

The rank of \mathbf{X} equals the number of its non-zero singular values.

In floating-point arithmetic “non-zero” is a judgement call, so np.linalg.matrix_rank counts the singular values above a small tolerance. A matrix with singular values (5, 3, 10^{-16}) is numerically rank 2, and the difference between 10^{-16} and 0 is the difference between an inverse that is merely enormous and one that does not exist.

Reduced or Economy SVD

Note that the sandwiched matrix \mathbf{\Sigma} holding the singular values is very sparse². To save space—without losing any information—numerical software usually returns the economy SVD.

² Most of the elements are zero; Wiki.

Short-fat

In this case, the matrix \mathbf{X} for decomposition has more columns than rows, i.e., m\leq n. \begin{bmatrix}\quad\\\quad\quad\quad{\Large\mathbf{X}}\quad\quad\quad\\\quad\end{bmatrix} =\begin{bmatrix}\\\quad\quad{\Large\mathbf{U}}\quad\quad\\\quad\end{bmatrix} \begin{bmatrix} \begin{array}{c|c} \begin{matrix} \ddots & &\\ & {\large\mathbf{\hat\Sigma}} & \\ & & \ddots \end{matrix} & \quad{\huge\mathbf{0}}\quad \end{array} \end{bmatrix} \begin{bmatrix} \begin{array}{c|c}{\Large\mathbf{\hat V}}&{\Large\mathbf{\hat V}}^\perp\end{array} \end{bmatrix}^T . The block \hat{\mathbf{V}}^\perp is multiplied by the zero block of \mathbf{\Sigma} and so contributes nothing. In the economy version the SVD therefore returns only \hat{\mathbf{V}}, whose size is n\times m (as opposed to n\times n).

Tall-skinny

In this case, the matrix \mathbf{X} for decomposition has more rows than columns, i.e., m\geq n. \begin{bmatrix}\quad\\\quad\\\quad{\Large\mathbf{X}}\quad\\\quad\\\quad\end{bmatrix} =\begin{bmatrix} \begin{array}{c|c}{\Large\hat{\mathbf{U}}} & {\Large\hat{\mathbf{U}}}^\perp \end{array} \end{bmatrix} \begin{bmatrix} \begin{array}{cc} \begin{matrix} \ddots & & \\ & {\large\hat{\mathbf{\Sigma}}} & \\ & & \ddots \end{matrix}\\\hline \begin{matrix}\\ \quad{\huge\mathbf{0}}\quad \\\end{matrix} \end{array} \end{bmatrix} \begin{bmatrix}\quad\\\quad\\\quad\quad\quad{\Large\mathbf{V}}\quad\quad\quad\\\quad\\\quad\end{bmatrix}^T. Symmetrically, \hat{\mathbf{U}}^\perp is annihilated, so the economy version returns only \hat{\mathbf{U}}, whose size is m\times n (as opposed to m\times m).

This is the case that matters for data: a tall-skinny \mathbf{X} is a data set with many more observations than features.

Numpy Implementation

The SVD algorithm is involved, so we skip it and use the implementation in numpy.linalg.svd³.

³ read the documentation.

Two things to notice about NumPy’s interface, both of which trip people up:

svd returns S as a 1D array of singular values, not a matrix. Use np.diag(S) if you need \mathbf{\Sigma}.
svd returns \mathbf{V}^T, not \mathbf{V}. The rows of the third output are the right singular vectors; the columns of its transpose are.

Note that U.T @ U is the identity but U @ U.T is not, because the economy U is 7\times2: it has orthonormal columns, but they cannot span \mathbb{R}^7. Only the full m\times m U from full_matrices=True is an orthogonal matrix in both directions.

The relation to eigenvalues

The SVD and the eigendecomposition of 3.1 are not two unrelated tools. Substituting \mathbf{X} = \mathbf{U}\mathbf{\Sigma}\mathbf{V}^T,

\mathbf{X}^T\mathbf{X} = \mathbf{V}\mathbf{\Sigma}^T\mathbf{U}^T\mathbf{U}\mathbf{\Sigma}\mathbf{V}^T = \mathbf{V}\,\mathbf{\Sigma}^T\mathbf{\Sigma}\,\mathbf{V}^T,

which is a spectral decomposition of the symmetric matrix \mathbf{X}^T\mathbf{X}. So:

the right singular vectors of \mathbf{X} are the eigenvectors of \mathbf{X}^T\mathbf{X};
the eigenvalues of \mathbf{X}^T\mathbf{X} are \sigma_k^2;
and, by the same argument on \mathbf{X}\mathbf{X}^T, the left singular vectors are its eigenvectors.

This single identity is the whole bridge to PCA, because \mathbf{X}^T\mathbf{X} is a covariance matrix once \mathbf{X} has been centred. It also explains the condition number: \kappa(\mathbf{X}) = \sigma_{\max}/\sigma_{\min}, and \kappa(\mathbf{X}^T\mathbf{X}) = \kappa(\mathbf{X})^2—which is precisely why 3.2 warned against forming A^\top A.

Exercise: singular values are square roots of eigenvalues

Verify the identity above. X is a 20\times4 matrix and S holds its singular values; produce the eigenvalues of \mathbf{X}^T\mathbf{X} from S alone, in descending order. Once you have it graded, run np.linalg.eigvalsh(X.T @ X)[::-1] and see that the two agree.

Matrix Approximation

We now come to a significant application of the SVD: approximating a very large matrix by an optimal low-rank one. For any desired rank r, a rank-r approximation is obtained by keeping the leading r singular values and singular vectors and discarding the rest.

Recall that the SVD arranges the singular values in \mathbf{\Sigma} in decreasing order of importance: \sigma_1\geq\sigma_2\geq\ldots\geq\sigma_{\min\{m,n\}}\geq0. If the matrix \mathbf{X} is not full-rank to begin with, some singular values are exactly zero. More often the singular values decay so rapidly that the first few carry nearly all of the information in \mathbf{X}.

A notable property of the SVD is that it writes the input matrix as a sum of rank-1 matrices: \mathbf{X}=\sum_{k=1}^{\min\{m,n\}}\sigma_k\mathbf{u}_k\mathbf{v}_k^T, \tag{15.1} where \sigma_k is the k^{\text{th}} diagonal entry of \mathbf{\Sigma}, and \mathbf{u}_k and \mathbf{v}_k are the k^{\text{th}} columns of \mathbf{U} and \mathbf{V}, respectively.

For any rank r we form the rank-r approximation \tilde{\mathbf{X}} by truncating the sum (Equation 15.1) after r terms.

The Eckart-Young theorem

That truncation is not merely a good rank-r approximation; it is the best possible one. Among all matrices \mathbf{B} of rank at most r, \tilde{\mathbf{X}} = \operatorname*{argmin}_{\operatorname{rank}(\mathbf{B})\le r}\|\mathbf{X}-\mathbf{B}\|, simultaneously for the Frobenius norm and the spectral norm (more). The error it achieves is a formula in the discarded singular values: \|\mathbf{X}-\tilde{\mathbf{X}}\|_F=\sqrt{\sum_{k>r}\sigma_k^2},\qquad \|\mathbf{X}-\tilde{\mathbf{X}}\|_2=\sigma_{r+1}. So the singular values do not merely rank the components—they tell you in advance exactly how much you lose by dropping each one. You can choose r before computing anything.

Exercise: measure the truncation error

Build the rank-3 truncation of X and report its Frobenius distance from X. Then compare your number to \sqrt{\sum_{k>3}\sigma_k^2} and see that Eckart-Young was not exaggerating.

Xr = U[:, :r] @ np.diag(S[:r]) @ Vt[:r, :]
answer = np.linalg.norm(X - Xr, 'fro')

The slicing pattern is worth memorising: columns of U, entries of S, rows of Vt. Since Vt is \mathbf{V}^T, its rows are the right singular vectors. A common bug is slicing Vt[:, :r] instead, which silently produces a matrix of the wrong shape or, worse, the right shape and the wrong answer.

Image Compression

We demonstrate matrix approximation using image compression. A gray-scale image is modeled as an m\times n matrix, where m and n are vertical and horizontal pixel directions, respectively. Each of the mn pixels contains a gray-scale value, depending on the chosen color-depth: 8-bit, 16-bit, etc.

The image below is 750\times1125, so it has 843{,}750 numbers and at most 750 non-zero singular values. Storing a rank-r approximation costs only r(m+n+1) numbers: at r=10 that is about 2\% of the original. Move the slider and watch how few components a recognisable rose needs.

viewof r = Inputs.range([1, 120], {step: 1, label: tex`r`, value: 10})

The reason this works at all is that natural images are not random. A random matrix has singular values that decay slowly, so no small r suffices; a photograph has enormous redundancy between neighbouring rows and columns, and its singular values fall off a cliff. Plotting them makes the point better than any argument:

The left panel is the whole story of dimensionality reduction in one picture: compressibility is exactly slow rank, and slow rank is exactly fast singular-value decay. The right panel says how many components you need to keep a chosen fraction of the total “energy” \sum\sigma_k^2, which is the same computation as the explained variance we are about to meet in PCA.

Low-rank approximation

You compute the SVD of a 1000\times800 matrix \mathbf{X} and find that \sigma_{26}/\sigma_1 < 10^{-10}, with every later singular value smaller still. Which conclusion is best supported?

\mathbf{X} has exactly 25 non-zero rows.\mathbf{X} is numerically close to a matrix of rank 25, and truncating the SVD after 25 terms gives the best rank-25 approximation in both the Frobenius and the spectral norm.Truncating after 25 terms discards most of \sum_k\sigma_k^2, so the approximation will be poor.The condition number of \mathbf{X} is close to 1.

Rank is the count of non-zero singular values, so a sharp drop after \sigma_{25} says \mathbf{X} sits very near the set of rank-25 matrices. The Eckart-Young theorem then says the truncation is not just a good approximation but the optimal one at that rank, simultaneously for both norms—an unusually strong guarantee, since a minimiser for one norm is rarely a minimiser for another.

The other options fail for concrete reasons:

Non-zero rows. Rank is a statement about linear dependence among rows or columns, not about how many entries are non-zero. A matrix of all ones has rank 1 and no zero entries at all.
A poor approximation. The error is \|\mathbf{X}-\tilde{\mathbf{X}}\|_F=\sqrt{\sum_{k>25}\sigma_k^2}. Every term in that sum is bounded by (10^{-10}\sigma_1)^2, so the discarded energy is negligible next to \sigma_1^2. The approximation is superb, which is exactly why the compression is worth doing.
Condition number. \kappa(\mathbf{X})=\sigma_{\max}/\sigma_{\min}, and the tail of tiny singular values makes that ratio at least 10^{10}. Rapid decay and good conditioning are opposites: this matrix is beautifully compressible and terrible to solve with.

Principal Component Analysis (PCA)

PCA is the central application of the SVD. From here on, let \mathbf{X} be an m\times n matrix in which each row is an observation and each column is a feature⁴. The principal components provide an alternative, orthogonal coordinate system centred on the mean of the data. Along these directions the data reveal their maximum variation.

⁴ This is how data are entered in a spreadsheet, and how a DataFrame is laid out. It is the opposite of the column-wise reading we used for the SVD definition above.

The algorithm

Step 1: centre the data. Compute the mean of each column—the average row— \overline{x}_j=\frac{1}{m}\sum_{i=1}^m x_{ij}, and assemble it into a matrix by repeating that row m times: \overline{\mathbf{X}}=\begin{bmatrix}1\\1\\\vdots\\1\end{bmatrix}_{m\times 1}\overline{\mathbf{x}}_{1\times n}= \begin{bmatrix} \overline{x}_1 & \overline{x}_2 & \ldots & \overline{x}_n\\ \overline{x}_1 & \overline{x}_2 & \ldots & \overline{x}_n\\ \vdots & \vdots & \ddots & \vdots\\ \overline{x}_1 & \overline{x}_2 & \ldots & \overline{x}_n\\ \end{bmatrix}_{m\times n}. Subtracting gives the mean-centred data \mathbf{B}=\mathbf{X}-\overline{\mathbf{X}}.

Centring is week 2 coming back

You never build \overline{\mathbf{X}} in NumPy. X.mean(axis=0) has shape (n,), and broadcasting stretches it up the rows for free:

B = X - X.mean(axis=0)      # (m, n) - (n,)  ->  (m, n)

That one line is the outer product \mathbf{1}\overline{\mathbf{x}} above, computed without ever materialising it. And by 3.2, \overline{\mathbf{X}} is the projection of \mathbf{X} onto the span of the all-ones vector, so centring is subtracting a projection: \mathbf{B} = (\mathbf{I} - \tfrac{1}{m}\mathbf{1}\mathbf{1}^\top)\mathbf{X}. Every column of \mathbf{B} is orthogonal to \mathbf{1}, which is just the statement that its entries sum to zero.

Step 2: factor the centred matrix. \mathbf{B}=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^T. In applications \mathbf{B} is tall-skinny (m\geq n: more observations than features). If the singular values are \sigma_1\geq\sigma_2\geq\ldots\geq\sigma_n, then the principal components are the columns of \mathbf{V}, and the variance captured along the k^{\text{th}} of them—the k^{\text{th}} principal value—is \lambda_k = \frac{\sigma_k^2}{m-1}.

Why \sigma_k^2/(m-1) is a variance

This is the step that makes PCA more than a recipe, and it uses nothing you have not already seen. The sample covariance matrix of the data is \mathbf{C}=\frac{1}{m-1}\mathbf{B}^T\mathbf{B}, an n\times n symmetric matrix whose (j,k) entry is the covariance between features j and k. This formula requires centred data—that is the only reason step 1 exists. Now substitute the SVD, exactly as we did earlier in this chapter: \mathbf{C}=\frac{1}{m-1}\mathbf{V}\mathbf{\Sigma}^T\mathbf{U}^T\mathbf{U}\mathbf{\Sigma}\mathbf{V}^T =\mathbf{V}\left(\frac{\mathbf{\Sigma}^T\mathbf{\Sigma}}{m-1}\right)\mathbf{V}^T. The right-hand side is a spectral decomposition of \mathbf{C} with an orthonormal \mathbf{V}. Reading it off:

The eigenvectors of the covariance matrix are the right singular vectors of the centred data. They are the principal components.
The eigenvalues of the covariance matrix are \sigma_k^2/(m-1). They are the variances along those directions.
Because \mathbf{V} is orthonormal, the principal components are mutually orthogonal: PCA hands you a coordinate system in which the features are uncorrelated.

So there are two routes to the same answer—eigh on \mathbf{C}, or svd on \mathbf{B}—and they are the same computation written twice. Prefer the SVD route. Forming \mathbf{B}^T\mathbf{B} squares the condition number (3.1), so the small principal values, which are the ones you were hoping to identify as negligible, are exactly the ones you compute worst.

The explained variance ratio of component k is \frac{\sigma_k^2}{\sum_{j}\sigma_j^2}, and its running total is what you plot to choose how many components to keep. It is the same “cumulative energy” curve we drew for the rose.

Skipping the centring does not give you “PCA without centring”

It gives you a different decomposition that answers a different question. Without centring, the first right singular vector chases the direction of the mean rather than the direction of greatest spread, because the mean offset is usually far larger than any deviation around it. The exercise below makes this concrete: shifting one feature by 50 leaves the centred first component untouched and rotates the uncentred one onto that axis almost exactly.

Standardising—dividing each centred column by its standard deviation—is a further, optional choice. It matters whenever features are measured in different units, since otherwise “maximum variance” just means “measured in the smallest units”. Income in dollars will dominate age in years every time.

Python Implementation

We generate a two-dimensional cloud with a known shape and recover its principal axes.

The two principal values come back at roughly 4 and 0.25, which are \sigma^2 for the axes we built the cloud with. PCA recovered the generating structure from the data alone, and the arrows point along the axes of the cloud rather than along the coordinate axes, because those are the directions in which the data actually vary.

Exercise: covariance eigenvalues from singular values

B is the centred data, S its singular values, and C = np.cov(X, rowvar=False) the covariance matrix. Produce the eigenvalues of C—in descending order—using S and m only. Check yourself afterwards against np.linalg.eigvalsh(C)[::-1].

Exercise: what centring is actually for

Y is the same data with feature 0 shifted by 50—a change of units, nothing more. The variable v1_centred holds the first right singular vector of the centred Y; it comes out as (0.944, 0.330, -0.015, -0.001), unchanged by the shift. Now compute the first right singular vector of the uncentred Y and compare.

PCA and centring

You run PCA by taking the SVD of a data matrix \mathbf{X} (rows = observations, columns = features) but forget to subtract the column means. Feature 1 is recorded in a unit that puts its values near 1000, while every feature varies by only about \pm1 around its own mean. What happens?

Nothing changes: the SVD is unaffected by shifting a column by a constant.np.linalg.svd raises a LinAlgError, since the decomposition requires centred input.The first right singular vector points almost exactly along the feature-1 axis, because \sigma_1 is now measuring distance from the origin instead of spread about the mean.The principal values come out right, but every component has its sign flipped.

The SVD always returns the directions in which \|\mathbf{X}\mathbf{v}\| is largest, and that is a statement about the raw vectors, measured from the origin. When one column sits near 1000, the dominant feature of the uncentred cloud is simply that it is a long way from the origin in that direction, and \mathbf{v}_1 locks onto it. The genuine structure—deviations of size \pm1—is buried three orders of magnitude below.

Centring moves the origin to the centroid, so \|\mathbf{B}\mathbf{v}\|^2 becomes a sum of squared deviations, i.e. a variance. Only then is “maximum \sigma” the same as “maximum spread”, and only then does \mathbf{B}^T\mathbf{B}/(m-1) deserve the name covariance matrix.

Why the others are wrong: the SVD is emphatically not shift-invariant (that is the whole point of step 1); np.linalg.svd runs happily on any finite matrix and will never warn you that you skipped centring, which is what makes this bug so easy to ship; and the damage is to the directions and the variances themselves, not merely to a sign convention—signs are arbitrary in every SVD, centred or not.

Note that centring is mandatory while standardising (dividing each centred column by its standard deviation) is a modelling choice. Standardise when the features carry different units, or “maximum variance” will just report whichever feature was measured in the smallest units.

Where this goes

Three later lectures are cashing cheques written in this chapter.

Nov 3—ridge regression. Write the design matrix as A=\mathbf{U}\mathbf{\Sigma}\mathbf{V}^T. Ordinary least squares weights the k^{\text{th}} direction by 1/\sigma_k, so a tiny \sigma_k blows the coefficients up. Ridge replaces that weight by \frac{\sigma_k}{\sigma_k^2+\lambda}, which agrees with 1/\sigma_k when \sigma_k is large and collapses to nearly zero when \sigma_k is small. Ridge does not shrink all coefficients equally; it shrinks hardest along the directions the data explored least, and the singular values are the measuring stick. That sentence is unavailable to anyone who has not seen the SVD.

Nov 17—the curse of dimensionality. As dimension grows, the distances between random points concentrate: the ratio of the spread of pairwise distances to their average shrinks toward zero, so “nearest” stops being meaningfully different from “farthest” and k-NN degrades. Problem 10 in 3.4 has you measure this.

Dec 8—clustering. k-means is defined entirely in terms of Euclidean distance, so it inherits that problem. Projecting onto the top few principal components first discards the low-variance directions—which are mostly noise—and restores some meaning to the distances that k-means is about to trust.