Eigenfaces, the two issues with PCA, and computing it fast

Before looking at PCA's shortcomings, L8 shows a real application of it, called eigenfaces. The data is the extended Yale face data set: about 20,000 grayscale photographs of faces, all the same size and all facing the camera.

From an image to a vector L8 · 1:21

Each image is a grid of pixels, say \(a \times b\), and each pixel holds a grayscale value between 0 and 255. Read the grid in a fixed order (row after row, for instance) and stack the values into one long column. The image becomes a vector with \(ab\) entries. The Yale images have about 32,000 pixels, so each face becomes a point in \(\mathbb{R}^{32{,}000}\).

The lecture writes the image size as \(m \times n\). These notes keep \(n\) for the number of data points.

The data set is now 20,000 points in 32,000 dimensions, and standard PCA from Part 3 applies unchanged: build the covariance matrix and find its eigenvectors. The covariance matrix is \(d \times d\), here \(32{,}000 \times 32{,}000\), and each eigenvector is itself a vector in \(\mathbb{R}^{32{,}000}\). The best line through the data is spanned by a 32,000-dimensional vector.

What the principal components look like

Because the reading order is fixed, vectorizing can be undone: any 32,000-dimensional vector can be folded back into an \(a \times b\) picture. Doing this to the eigenvectors shows the principal components as images, and they look like faces.

• PC1 is the direction of greatest variance (Part 3). As an image it is a prototype face: the single template that best captures what all the images share, their "faceness".
• Later components add detail. By some \(k\)-th component the nose and the mouth are sharper. The common template already accounts for most of the variance, so the other directions pick up the remaining characteristics in which faces differ, and each adds more detail.

Suppose a company builds a face recognition system for its employees (the lecture says 1,000 employees, then 20,000 to match the data set). Storing 32,000 pixel values per person looks wasteful. Is there a lower-dimensional representation that still does the job? PCA can tell us.

Keep the top \(r\) principal components. In this setting, keeping \(r\) components, running the residue procedure of Part 2 for \(r\) rounds and keeping \(r\) eigendirections all mean the same thing. Then rebuild the proxy of each image from those components alone 6:47.

• One component. The coefficient is the dot product \(\mathbf{x}^\top\mathbf{w}_1\) and the proxy is \((\mathbf{x}^\top\mathbf{w}_1)\,\mathbf{w}_1\). Every proxy is a scaled copy of the PC1 image.
• Two components. Take the dot products with \(\mathbf{w}_1\) and \(\mathbf{w}_2\), and combine the two component images with those numbers as weights.

Standard practice, not discussed in the lecture: if the images were centred before PCA (Part 2), the mean image is added back to \(\tilde{\mathbf{x}}\).

The lecture shows this for a test image that is not in the data set, for example a photo taken by a CCTV camera as someone walks into the office. Its reconstructions use \(r = 25, 50, 100, \dots\) up to 1,600 components 8:26.

Now feed in an image unlike anything in the data set: a dog's face, also 32,000 pixels. With 25 components, its proxy is the projection onto the 25-dimensional subspace that best represents human faces. What comes back is the best human-face version of the dog. At \(r = 100\), 200 and 400 the result is still poor. Only at about 1,600 components does the dog look like the dog.

The reason 11:37 is that the components were learnt from human faces only, so they capture the directions that matter for human faces. Anything else needs many more of them. How many components are enough depends on what you want to do.

• Human or dog? Even 25 components suffice. On the 25-dimensional face subspace the dog looks like neither a human nor a dog, and that difference is already enough to tell the two apart.
• Which human? About 200 to 300 components. At that level the dog still looks nothing like a dog.
• With 1,600 components both the human faces and the dog's face are rebuilt well. The lecture's reading: these directions hold enough information to represent perhaps any face, anything with eyes, a nose and a mouth.
• A tree would probably need more, perhaps 2,000 or 3,000. The lecturer says this is a guess he has not tried on this data set.
• An image of random pixels has no structure, so it is not even clear how to judge whether a reconstruction resembles it. Rebuilding it exactly needs all 32,000 components.

Keeping all 32,000 components reconstructs any image exactly, because no information has been thrown away. The real question is when far fewer components are enough. For human faces, 200 to 300 do the job.

Eigenfaces in a recognition system L8 · 14:22

In practice PCA is used as a preprocessing step. Store each training face only as its projections onto the top few eigendirections, and train a supervised learning algorithm on those short vectors. A new test image is projected onto the same directions, and the trained algorithm makes its prediction from the projections.

This works because natural images have structure. Most of the information about human faces lies in an underlying linear subspace of far lower dimension than 32,000.

The rest of L8 raises two concerns about PCA, to be fixed over the next lectures. Understanding and fixing them leads to an idea that is useful far beyond PCA.

The first concern is the running time; the lecture calls it a concern more than a problem. PCA builds the covariance matrix and then finds its eigenvectors and eigenvalues, and the eigendecomposition is the expensive step 16:37. Special-purpose numerical linear algebra solvers can do better in particular cases. As a general figure, though, the eigendecomposition of a \(d \times d\) matrix takes time of order \(d^3\).

This becomes a problem when \(d\) is large. In eigenfaces \(d \approx 32{,}000\), while \(n\) is only about 20,000: there are fewer data points than features. With \(d\) in the tens of thousands or the millions, \(d^3\) is enormous.

Face recognition is the lecture's example of where this happens. The fix comes in L9. First, the second issue.

PCA looks for important linear combinations of features: weight plus height, say, or two times weight plus three times height. Features need not be related linearly. Add body mass index to weight and height: it depends on them through a quadratic, weight divided by height squared. (The lecturer warns he may have the formula wrong; that is the standard one.)

Picture three features plotted in \(\mathbb{R}^3\), with all the data points lying on a curved surface, like a blanket draped in space. A relationship such as \(f_3 = f_1 f_2 + f_2^2\) between the features could produce this; it is the lecture's example. The lecture calls such a curved space inside \(\mathbb{R}^3\) a manifold and skips the formal definition. The data has real structure, but the structure is curved, because the features themselves are related non-linearly.

PCA only looks for relationships in which each feature gets a weight and the weighted features are added up. It still produces an output on the curved surface, and the output misleads:

1. PCA returns the best plane for the data. The surface is curved, so the reconstruction error with respect to that plane is still large, far from insignificant.
2. That large error suggests a third direction matters as well. Add it, and the error drops to 0, but only because three dimensions can rebuild any data set in \(\mathbb{R}^3\).
3. PCA never learns that the data is essentially two-dimensional.

Data need not lie in a low-dimensional linear subspace. It can be low-dimensional and curved. 23:44

Here a subspace means a line or a plane through the origin, or the higher-dimensional analogue. The widget below shows the same failure one dimension lower: points on a curve in the plane, where PCA can only offer lines.

The two issues look unrelated: one is about computation, the other about the shape of the data. The punch line of L8 is that they have the same solution 24:51. Solve issue 1 intelligently, and the solution can be reused for issue 2.

The lecture keeps the suspense. The plan is to solve issue 1 first (L9, below), make some observations along the way, and use them to solve issue 2 in Part 5.

Issue 1 is a large \(d\), in particular \(d \gg n\): many more features than data points. The first step is some notation that makes the problem easier to handle. Stack the data points as the columns of one matrix:

The covariance matrix is \(C = \frac1n\sum_i \mathbf{x}_i\mathbf{x}_i^\top\) (Part 3). That sum is exactly the matrix product \(XX^\top\): a \(d \times n\) matrix times an \(n \times d\) matrix whose rows are the data points.

The lecture leaves the identity as an exercise. Write \((\mathbf{x}_i)_p\) for the \(p\)-th coordinate of \(\mathbf{x}_i\). Entry \((p, q)\) of \(XX^\top\) is row \(p\) of \(X\) dotted with row \(q\) of \(X\), which is \(\sum_i (\mathbf{x}_i)_p(\mathbf{x}_i)_q\). That is also entry \((p, q)\) of \(\sum_i \mathbf{x}_i\mathbf{x}_i^\top\).

Dimension check: \(d \times n\) times \(n \times d\) gives \(d \times d\), the size of \(C\) 2:47. Some definitions divide by \(n - 1\) instead of \(n\); it does not matter, and this course uses \(n\). As in Parts 2 and 3, the data is taken to be centred.

PCA wants the eigenvectors and eigenvalues of \(C\). Arrange the eigenvalues in decreasing order, \(\lambda_1 \ge \lambda_2 \ge \dots\), and let \(\mathbf{w}_k\) be the eigenvector for the \(k\)-th largest, \(\lambda_k\). An eigenvector is a direction that the matrix only scales. The scale may be negative, which flips the vector but keeps its line.

Where do these eigenvectors live? They are \(d\)-dimensional vectors, but more can be said 5:02. Replace \(C\) by its definition. Since \(\mathbf{x}_i^\top\mathbf{w}_k\) is a number, \(\mathbf{x}_i\mathbf{x}_i^\top\mathbf{w}_k = (\mathbf{x}_i^\top\mathbf{w}_k)\,\mathbf{x}_i\). Then divide both sides by \(\lambda_k\):

Dividing by \(\lambda_k\) needs \(\lambda_k \neq 0\), which the lecture assumes from here on. The right-hand side takes each data point \(\mathbf{x}_i\), multiplies it by a number and adds the results. So \(\mathbf{w}_k\) is a linear combination of the data points. Each eigenvector combines them with its own weights.

Collect the weights for \(\mathbf{w}_k\) into a vector \(\boldsymbol\alpha_k\) with one entry per data point. Multiplying \(X\) by it combines the columns of \(X\), which are the data points, with those weights:

So finding the eigenvectors is the same as finding the right weights. Once we have \(\mathbf{w}_k\) we know what to do: project each point onto it to get the compressed representation. And once we have \(\boldsymbol\alpha_k\), we have \(\mathbf{w}_k = X\boldsymbol\alpha_k\).

A chicken-and-egg problem L9 · 10:14

The derivation even says what the weights are: \(\alpha_{ki} = \mathbf{x}_i^\top\mathbf{w}_k / (n\lambda_k)\). Computing them that way needs \(\mathbf{w}_k\), the very thing we are trying to find. So set this formula aside and look for a different way to get \(\boldsymbol\alpha_k\), one for each \(k\). The new route is worth it only if it avoids the \(O(d^3)\) eigendecomposition of \(C\).

The lecture works through the algebra for completeness; what comes out of it matters more. Start from the eigenvector equation \(C\mathbf{w}_k = \lambda_k\mathbf{w}_k\) and substitute the two facts \(C = \frac1n XX^\top\) and \(\mathbf{w}_k = X\boldsymbol\alpha_k\). Then move \(n\) to the right:

Now pre-multiply both sides by \(X^\top\) 13:38. Matrix multiplication is not commutative, so factors cannot be swapped. It is associative, so the brackets can be placed anywhere:

The matrix \(X^\top X\) is \(n \times d\) times \(d \times n\), so it is only \(n \times n\). Give it a name:

It is enough to find an \(\boldsymbol\alpha_k\) with \(K\boldsymbol\alpha_k = (n\lambda_k)\,\boldsymbol\alpha_k\) 15:54: multiplying both sides of it by \(K\) gives the equation above. And that is an eigen-equation. Applying \(K\) to \(\boldsymbol\alpha_k\) must only scale it, by \(n\lambda_k\). So \(\boldsymbol\alpha_k\) should be an eigenvector of \(K\) with eigenvalue \(n\lambda_k\).

Which length? L9 · 17:01

An eigen-equation does not fix the length of its solution. If \(K\mathbf{u} = (n\lambda_k)\,\mathbf{u}\), then also \(K(2\mathbf{u}) = (n\lambda_k)(2\mathbf{u})\), and the same holds for any multiple of \(\mathbf{u}\). These multiples cannot all be \(\boldsymbol\alpha_k\): we need the specific weights that give a \(\mathbf{w}_k\) of length 1. The length condition on \(\mathbf{w}_k\) turns into a condition on \(\boldsymbol\alpha_k\) 18:08:

One gap remains. The equation involves \(\lambda_k\), an eigenvalue of \(C\), but we only intend to decompose \(K\). If we solve the eigen-equation of \(K\), how are its eigenvalues related to the ones we would have got from \(XX^\top\)?

The answer is a theorem of linear algebra, which the lecture states as a fact.

The two matrices have different sizes, but both are built from the same \(d \times n\) matrix \(X\), which is why their eigenvalues are related. Some details 22:37:

• Each matrix has at most \(\min(d, n)\) non-zero eigenvalues (another linear algebra fact).
• So when \(d\) is large, \(XX^\top\) has many zero eigenvalues, with no counterpart in \(X^\top X\).
• If \(n < d\) and \(X\) has full rank, \(X^\top X\) has \(n\) non-zero eigenvalues, and they match the top \(n\) eigenvalues of \(XX^\top\).

The lecture suggests proving the fact in two steps with the singular value decomposition, and leaves it to you. A shorter standard argument: suppose \(X^\top X\boldsymbol\beta = \mu\boldsymbol\beta\) with \(\mu \neq 0\) and \(\boldsymbol\beta \neq \mathbf{0}\). Then

So \(X\boldsymbol\beta\) is a non-zero vector, and \(\mu\) is an eigenvalue of \(XX^\top\) with eigenvector \(X\boldsymbol\beta\). Swapping the roles of \(X\) and \(X^\top\) gives the other direction. Notice that this eigenvector, \(X\boldsymbol\beta\), is again a combination of the data points.

Now put everything together. We want the unit eigenvectors \(\mathbf{w}_1, \dots, \mathbf{w}_l\) of \(C\), with eigenvalues \(\lambda_1 > \lambda_2 > \dots > \lambda_l\). Compare the three matrices involved:

So \(K\boldsymbol\beta_k = (n\lambda_k)\,\boldsymbol\beta_k\): \(\boldsymbol\beta_k\) satisfies the eigen-equation that \(\boldsymbol\alpha_k\) must satisfy. Is \(\boldsymbol\beta_k = \boldsymbol\alpha_k\)? Check the length condition. Using \(K\boldsymbol\beta_k = n\lambda_k\boldsymbol\beta_k\) and \(\boldsymbol\beta_k^\top\boldsymbol\beta_k = 1\):

That is \(n\lambda_k\), not 1 in general, so \(\boldsymbol\beta_k\) has the right direction and the wrong length. Rescale it. A multiple of an eigenvector is still an eigenvector, so \(\boldsymbol\alpha_k = \boldsymbol\beta_k/\sqrt{n\lambda_k}\) still satisfies \(K\boldsymbol\alpha_k = (n\lambda_k)\,\boldsymbol\alpha_k\), and the square root is applied twice in the length condition:

The square root is real. The eigenvalues of \(C\) are never negative, because \(\mathbf{w}^\top C\mathbf{w} = \frac1n\sum_i(\mathbf{x}_i^\top\mathbf{w})^2 \ge 0\) (Part 2), so the non-zero ones are positive.

1. PCA when d ≫ n
2. Input \(\mathcal{D} = \{\mathbf{x}_1, \dots, \mathbf{x}_n\}\), \(\mathbf{x}_i \in \mathbb{R}^d\), with \(d \gg n\): too many features, the setting where the cost of PCA hurts.
3. Compute \(K = X^\top X \in \mathbb{R}^{n \times n}\).
4. Eigendecompose \(K\): unit eigenvectors \(\boldsymbol\beta_1, \dots, \boldsymbol\beta_l\) with eigenvalues \(n\lambda_1 \ge \dots \ge n\lambda_l\). This costs \(O(n^3)\).
5. Set \(\boldsymbol\alpha_k = \boldsymbol\beta_k / \sqrt{n\lambda_k}\) for \(k = 1, \dots, l\).
6. If you want the directions themselves, \(\mathbf{w}_k = X\boldsymbol\alpha_k\).

The old route computed \(C\) and eigendecomposed it in \(O(d^3)\). Step 2 is still an eigendecomposition, but of an \(n \times n\) matrix, so it takes \(O(n^3)\). When \(d \gg n\) that can be far cheaper. We reach the \(\mathbf{w}_k\) by a different road: decompose a related, smaller matrix, then convert its eigenvectors into the weights that combine the data points into the eigenvectors of \(C\) 32:26.

Two remarks from the lecture 33:29:

• It helps only when \(d \gg n\). If \(n\) is larger than \(d\), you might as well decompose the covariance matrix directly.
• An eigendecomposition cannot be avoided. The problem has two sizes, \(d\) and \(n\), and the trick is simply to decompose the matrix whose side is the smaller of the two.

Not part of the lecture's count: forming \(C\) takes about \(nd^2\) multiplications and forming \(K\) about \(n^2 d\), so building \(K\) is also the cheaper step when \(d \gg n\).

L9 ends with an observation that the next lecture builds on. Look at one entry of \(K = X^\top X\). Row \(i\) of \(X^\top\) is \(\mathbf{x}_i^\top\) and column \(j\) of \(X\) is \(\mathbf{x}_j\), so:

In a sense \(K_{ij}\) captures how similar \(\mathbf{x}_i\) and \(\mathbf{x}_j\) are, in the dot-product sense. More importantly, after step 1 the algorithm works with \(K\) alone. To solve the PCA problem you need only the pairwise dot products between the data points, not necessarily the points themselves.

A one-line check of this claim (the lecture promises to make it precise next time): the coefficient of a data point \(\mathbf{x}_j\) along \(\mathbf{w}_k\), which is what compression stores, also comes from \(K\) alone.

Carry this into Part 5: PCA can be run from the matrix of pairwise dot products \(K\) alone.

Why would a speed-up for issue 1 help with curved data in issue 2? This observation, that only similarities between points are needed, is what lets the solution to issue 1 be turned into a solution to issue 2 35:46. It takes some creative thought, and the lecture calls it one of the very important ideas of classical machine learning. That is the subject of Part 5: feature maps, kernels and kernel PCA.

In summary, L8 identified two issues with PCA, and L9 solved the time-complexity issue by converting the original problem into a computationally simpler one whose solution gives back the solution to the original.

Eigenfaces

PCA on vectorized face images; the principal components, folded back into images, look like faces.

Vectorizing an image

Reading an \(a \times b\) pixel grid in a fixed order into a vector with \(ab\) entries. The same order reverses it.

Components kept, \(r\)

The number of top eigendirections retained; the same as running the residue procedure for \(r\) rounds.

Reconstruction from \(r\) components

\(\tilde{\mathbf{x}} = \sum_{k=1}^r(\mathbf{x}^\top\mathbf{w}_k)\mathbf{w}_k\), stored as \(r\) numbers per point.

Linear subspace

A line or plane through the origin, or its higher-dimensional analogue. PCA can only find these.

Manifold

The lecture's informal word for a curved space inside \(\mathbb{R}^d\) on which the data lies, like a draped blanket in \(\mathbb{R}^3\).

Time complexity \(O(d^3)\)

The general cost of an eigendecomposition of a \(d \times d\) matrix: time grows like the cube of \(d\).

Data matrix \(X\)

The \(d \times n\) matrix with the data points as columns, so that \(C = \frac1n XX^\top\).

\(K = X^\top X\)

The \(n \times n\) matrix of pairwise dot products, \(K_{ij} = \mathbf{x}_i^\top\mathbf{x}_j\). Standard name: the Gram matrix.

Weights \(\boldsymbol\alpha_k\)

The \(n\) numbers that combine the data points into the eigenvector, \(\mathbf{w}_k = X\boldsymbol\alpha_k\).

Shared non-zero eigenvalues

\(XX^\top\) and \(X^\top X\) have the same non-zero eigenvalues, at most \(\min(d, n)\) of them.

Rescaling \(\boldsymbol\beta_k\)

\(\boldsymbol\alpha_k = \boldsymbol\beta_k/\sqrt{n\lambda_k}\) turns a unit eigenvector of \(K\) into weights that give a unit \(\mathbf{w}_k\).