Feature maps, kernels and kernel PCA

Part 4 raised two issues with PCA. Issue 1 was cost. Finding the eigenvectors of the \(d \times d\) covariance matrix takes about \(O(d^3)\) time, which hurts when \(d \gg n\). The fix was to work with the \(n \times n\) matrix \(X^\top X\) instead and convert its eigenvectors into those of the covariance matrix.

These lectures take up issue 2, which Part 4 also flagged. PCA is good at finding linear relationships among features, and only those. What if the features are related non-linearly?

Take two-dimensional data in which every point lies on a circle centred at \((a, b)\). Run standard PCA. It first centres the data, which moves \((a, b)\) to the origin. Then it looks for the direction along which the projected points have the largest variance, which is the same as the smallest reconstruction error.

On a circle no direction stands out. Projecting onto any direction gives roughly the same variance. Which direction wins depends only on exactly how the points happen to sit around the circle. Say PCA picks some \(\mathbf{w}_1\); then \(\mathbf{w}_2\) is the perpendicular direction.

The eigenvalues are the variances along these directions. Here \(\mathbf{w}_1\) captures slightly more than 50% of the variance and \(\mathbf{w}_2\) slightly less. With the rule of thumb from Part 3 (keep enough directions to explain 95% of the variance), PCA says both directions are needed. There is no compression.

Yet the data is obviously simple. Do we really need two numbers per point to describe it? That question drives the rest of the lecture.

Call the two features \(f_1\) and \(f_2\). Every data point lies on the circle of radius \(r\) centred at \((a, b)\), so it satisfies \((f_1 - a)^2 + (f_2 - b)^2 = r^2\). Expanding the squares:

This relationship is not linear: it contains \(f_1^2\) and \(f_2^2\). The picture says the same. Now map each two-dimensional point to a six-dimensional one:

Also fix one vector in \(\mathbb{R}^6\):

\(\mathbf{u}\) is built from the circle's centre and radius only. No individual data point, no \(f_1\) or \(f_2\), was used to define it. Its dot product with a lifted point is exactly the left side of \((\star)\):

So a point satisfies \((\star)\) exactly when \(\phi(\mathbf{x})^\top\mathbf{u} = 0\). Every data point, once mapped, is orthogonal to \(\mathbf{u}\).

\(\phi\) is defined for every point of \(\mathbb{R}^2\). Only points that lie on this particular circle land orthogonal to \(\mathbf{u}\).

Turn the example into a method. Choose a map \(\phi : \mathbb{R}^d \to \mathbb{R}^D\) into a higher dimension, \(D > d\), in the hope that the mapped data lies in a low-dimensional linear subspace of \(\mathbb{R}^D\). Finding the important directions of such a subspace is exactly what PCA does.

Mapping up makes \(D\) large, possibly much larger than \(n\). That is issue 1 again, and Part 4 already solved it: instead of the \(D \times D\) matrix built from \(\phi(X)\phi(X)^\top\), work with the \(n \times n\) matrix \(\phi(X)^\top\phi(X)\).

Writing \(\phi(X)\) for a matrix is a slight abuse of notation: it means \(\phi\) applied to each column (each data point) of \(X\).

Where did the non-linearity go? Into \(\phi\). The map holds whatever non-linear functions of the features we want to allow. After mapping, PCA only has to find linear structure, which it does well.

Curved in \(\mathbb{R}^d\) can be flat in \(\mathbb{R}^D\). Map with \(\phi\), then run ordinary PCA through the \(n \times n\) matrix.

The idea is sound, but it runs into a bottleneck at once. Suppose there are four features (say height, weight, age and gender) and we want \(\phi\) to capture every cubic relationship, such as height × weight × age or height × weight², all combinations whose powers add up to at most 3.

The lecture tallies the coordinates by degree: a constant, the 4 features themselves, the \(\binom{4}{2} = 6\) pairs \(f_1f_2, f_1f_3, \dots, f_3f_4\), and the \(\binom{4}{3} = 4\) triples. That gives \(1 + 4 + 6 + 4 = 15\) coordinates. With \(d\) features and relationships up to power \(p\):

The lecture's illustrations: with \(d = 10\) and 15th-power relationships, \(\phi\) would have about \(10^{15}\) coordinates instead of 10. Even two features with 20th-power relationships would need about \(2^{20}\). Spaces like that are simply too large to work with, and the problem gets worse as either \(d\) or \(p\) grows.

Standard counts, not from the lecture. The tally above counts products of distinct features only. Counting every monomial of degree at most \(p\), including powers such as weight² that the lecture also mentions, gives \(\binom{d+p}{p}\). \(d^p\) counts ordered products (\(f_1f_2\) and \(f_2f_1\) separately). For fixed \(p\) and large \(d\), all three grow like a constant times \(d^p\), which is the lecture's point. When \(p\) exceeds \(d\), as in the two extreme examples, \(d^p\) overstates the count by a lot.

The first row is the circle example: the map \(\phi\) above, with its 6 coordinates, is exactly the list of all monomials of degree at most 2 in two features. With a thousand features and cubic relationships, each point would need over \(10^8\) coordinates.

Whichever way the coordinates are counted, computing \(\phi(\mathbf{x})\) explicitly soon becomes impractical. Lecture 11 starts from this bottleneck.

\(\phi(\mathbf{x}) \in \mathbb{R}^D\) may be too expensive to compute. The obstacle is size: there are simply too many coordinates. Without \(\phi(\mathbf{x})\), is PCA in the feature space doomed?

Look at what the \(n \times n\) route actually uses. Its eigenvectors come from \(\phi(X)^\top\phi(X)\), and entry \((i, j)\) of that matrix is a dot product of two mapped points:

So we need the pairwise dot products, not the mapped vectors themselves. The obvious route is to map \(\mathbf{x}_i\) and \(\mathbf{x}_j\) with \(\phi\) and then take the dot product. The question is whether there is a direct route from \(\mathbf{x}_i, \mathbf{x}_j\) to the number \(\phi(\mathbf{x}_i)^\top\phi(\mathbf{x}_j)\) that never computes \(\phi\). The lecturer admits this "sounds like magic". Without such a route, capturing complicated \(p\)-th power relationships would be out of reach.

Let \(\mathbf{x} = [f_1, f_2]^\top\) and \(\mathbf{x}' = [g_1, g_2]^\top\), and consider a function pulled out of a hat: \((\mathbf{x}^\top\mathbf{x}' + 1)^2\). Expand it:

Each term is a product of something that depends only on \(\mathbf{x}\) and something that depends only on \(\mathbf{x}'\). Collect the \(\mathbf{x}\) parts into one vector and the \(\mathbf{x}'\) parts into another. A factor of 2 is split as \(\sqrt2 \cdot \sqrt2\), one \(\sqrt2\) on each side:

The two vectors have the same functional form: the same recipe applied to \(\mathbf{x}\) and to \(\mathbf{x}'\). That recipe is a feature map.

Mapping from 2 to 6 dimensions is not much of a saving; we could just compute \(\phi\) and work in \(\mathbb{R}^6\). The trick matters when \(d\) is large and the power \(p\) is high. The calculator below does both routes for any two points in \(\mathbb{R}^2\) and powers up to 4.

The same works in general. For \(\mathbf{x}, \mathbf{x}' \in \mathbb{R}^d\) and any power \(p \ge 1\), consider

The lecture states, without proof, that this is a valid function. Valid has a precise meaning here:

Validity is what lets us trust \(k\). If such a \(\phi\) exists, the matrix with entries \(k(\mathbf{x}_i, \mathbf{x}_j)\) really is \(\phi(X)^\top\phi(X)\), and its eigenvectors give PCA in the feature space. If no such \(\phi\) exists, that matrix means nothing for PCA. Not every function computes dot products in some space; an arbitrary formula usually does not.

One way to convince yourself that \(k\) is valid is to write the map down. An explicit \(\phi\) is a certificate of validity, exactly as for \(p = 2\) with two features.

Exercise from the lecture 12:31. Write out \(\phi\) explicitly for \(p = 3\) and \(p = 4\), and watch how its dimension grows with \(d\) and \(p\). A check (standard result): \((\mathbf{x}^\top\mathbf{x}'+1)^p\) needs one coordinate per monomial of degree at most \(p\), so \(\binom{d+p}{p}\) coordinates. For \(d = 2\) that is 3, 6, 10 and 15 for \(p = 1, \dots, 4\), which the calculator shows.

A second famous example, widely used in machine learning:

It looks like the Gaussian density, so it is sometimes called the Gaussian kernel; the standard name is the radial basis function (RBF) kernel. It is also valid: there is a \(\phi\) for which it computes dot products.

What is surprising is the dimension of that \(\phi\). For the polynomial kernel, \(D\) grew like \(d^p\). For the RBF kernel, \(\phi\) maps each point into an infinite-dimensional space. Infinity brings technicalities that the course sets aside. The working picture:

• A point becomes a function. A function takes any real input and returns a number, so you can think of it as an infinitely long vector, one coordinate per input (uncountably many). The RBF map sends each data point to such a function.
• Dot products become integrals. A dot product sums coordinate-by-coordinate products. For functions the sum over coordinates becomes an integral.
• Another way to see it. \(e^z\) has a power series containing every power of \(z\), with shrinking weights. So the RBF kernel mixes in polynomial relationships of every degree, each contributing a little less than the last.

Standard expansion, not worked in the lecture. For \(d = 1\) and \(\sigma = 1\): \(e^{-(x-x')^2/2} = e^{-x^2/2}\,e^{-x'^2/2}\,e^{xx'}\) and \(e^{xx'} = \sum_{m=0}^{\infty} \frac{x^m}{\sqrt{m!}}\cdot\frac{x'^m}{\sqrt{m!}}\). So \(\phi(x) = e^{-x^2/2}\,\big[1,\ x,\ \tfrac{x^2}{\sqrt{2!}},\ \tfrac{x^3}{\sqrt{3!}},\ \dots\big]\) works, with one coordinate for every power \(m\).

Valid functions get a name. "Kernel" means several things in mathematics; in machine learning it means this:

These are two examples among many. A different non-linear function may suit a given problem better. So the general question is: given some \(k : \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}\) and the claim that it is a kernel, how would you check?

Method 1: exhibit \(\phi\). Find a map \(\phi\) explicitly and show \(k(\mathbf{x}, \mathbf{x}') = \phi(\mathbf{x})^\top\phi(\mathbf{x}')\), as we did for the quadratic kernel. This is the standard method, but coming up with the map can be hard. Since not every function is a kernel, we want a test that does not require finding \(\phi\).

Method 2 is Mercer's theorem. The course states an informal, machine-learning version rather than the fully rigorous one, and does not prove it.

The "if and only if" is what makes the theorem powerful: it characterizes validity completely. \(K\) is symmetric automatically once \(k\) is, so condition (b) is about the signs of its eigenvalues.

Why symmetry is necessary

If \(k\) is valid, \(k(\mathbf{x}, \mathbf{x}') = \phi(\mathbf{x})^\top\phi(\mathbf{x}')\). The dot product is symmetric, so this equals \(\phi(\mathbf{x}')^\top\phi(\mathbf{x}) = k(\mathbf{x}', \mathbf{x})\). So if you can find a single pair with \(k(\mathbf{x}, \mathbf{x}') \ne k(\mathbf{x}', \mathbf{x})\), the claim that \(k\) is a kernel is dead. A symmetric function, though, is not automatically a kernel; condition (b) must also hold.

Why non-negative eigenvalues are necessary

If \(k\) is valid then \(K_{ij} = \phi(\mathbf{x}_i)^\top\phi(\mathbf{x}_j)\), so \(K = \phi(X)^\top\phi(X)\) for some \(\phi\), even if we do not know which. From Part 4, the non-zero eigenvalues of \(\phi(X)^\top\phi(X)\) are the same as those of \(\phi(X)\phi(X)^\top\). That matrix is \(n\) times the covariance matrix in the feature space, and its eigenvalues are \(n\) times the variance along the corresponding eigenvectors. A variance cannot be negative. So every eigenvalue of \(K\) is \(\ge 0\). Zero is allowed; negative is not.

The same fact in one line (standard): for any \(\mathbf{v} \in \mathbb{R}^n\), \(\mathbf{v}^\top K\mathbf{v} = \mathbf{v}^\top\phi(X)^\top\phi(X)\mathbf{v} = \lVert\phi(X)\mathbf{v}\rVert^2 \ge 0\), and \(\mathbf{v}^\top K\mathbf{v} \ge 0\) for all \(\mathbf{v}\) is the usual definition of positive semi-definite.

So both conditions are necessary, and that direction is easy. Mercer's theorem adds that together they are also sufficient. That direction is more involved and is not proved in the course, but it is the part that lets you certify a kernel without ever finding \(\phi\).

To rule a function out, one counterexample is enough: one asymmetric pair, or one data set whose kernel matrix has a negative eigenvalue.

Lecture 12 puts everything together into kernel PCA: PCA that can capture non-linear relationships between features, with kernels doing the mapping. The input is the data set \(\{\mathbf{x}_1, \dots, \mathbf{x}_n\}\), \(\mathbf{x}_i \in \mathbb{R}^d\), together with a kernel \(k : \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}\) already known to be valid (say, checked with Mercer's theorem). The kernel fixes some \(\phi : \mathbb{R}^d \to \mathbb{R}^D\), which we never compute.

Why the eigenvectors of \(K\) are enough

This is the argument of Part 4 with \(\phi(\mathbf{x}_i)\) in place of \(\mathbf{x}_i\). For now assume the mapped points are centred (fixed below). PCA in the feature space wants the eigenvectors of the feature-space covariance matrix:

1. \(\mathbf{w}_k\) is a combination of the mapped points. Substitute the sum for \(C\) and move \(\lambda_k \ne 0\) across: \(\mathbf{w}_k = \sum_i \frac{\phi(\mathbf{x}_i)^\top\mathbf{w}_k}{n\lambda_k}\,\phi(\mathbf{x}_i)\). So \(\mathbf{w}_k = \phi(X)\,\boldsymbol\alpha_k\) for some weight vector \(\boldsymbol\alpha_k \in \mathbb{R}^n\).
2. The weights solve an eigen-equation in \(K\). Put \(\mathbf{w}_k = \phi(X)\boldsymbol\alpha_k\) into \(\phi(X)\phi(X)^\top\mathbf{w}_k = n\lambda_k\mathbf{w}_k\) and multiply on the left by \(\phi(X)^\top\): \(K^2\boldsymbol\alpha_k = n\lambda_k K\boldsymbol\alpha_k\) with \(K = \phi(X)^\top\phi(X)\). Any \(\boldsymbol\alpha_k\) with \(K\boldsymbol\alpha_k = n\lambda_k\boldsymbol\alpha_k\) satisfies it.
3. Unit length fixes the scale. \(1 = \mathbf{w}_k^\top\mathbf{w}_k = \boldsymbol\alpha_k^\top\phi(X)^\top\phi(X)\boldsymbol\alpha_k = \boldsymbol\alpha_k^\top K\boldsymbol\alpha_k\).
4. Convert the solver's output. An eigen-solver returns unit vectors \(\boldsymbol\beta_k\) with \(K\boldsymbol\beta_k = n\lambda_k\boldsymbol\beta_k\) (the non-zero eigenvalues of \(K\) and \(nC\) agree). Then \(\boldsymbol\beta_k^\top K\boldsymbol\beta_k = n\lambda_k\), so dividing by \(\sqrt{n\lambda_k}\) gives \(\boldsymbol\alpha_k^\top K\boldsymbol\alpha_k = n\lambda_k / n\lambda_k = 1\).
5. The kernel supplies \(K\). \(K_{ij} = \phi(\mathbf{x}_i)^\top\phi(\mathbf{x}_j) = k(\mathbf{x}_i, \mathbf{x}_j)\). Nothing so far needs \(\phi\).

1. Kernel PCA, steps 1 and 2
2. Kernel matrix. Compute \(K \in \mathbb{R}^{n\times n}\) with \(K_{ij} = k(\mathbf{x}_i, \mathbf{x}_j)\) for all \(i, j\).
3. Eigen-decompose and normalize. Find unit eigenvectors \(\boldsymbol\beta_1, \dots, \boldsymbol\beta_l\) of \(K\) with eigenvalues \(n\lambda_1 \ge \dots \ge n\lambda_l\), for as many directions \(l\) as you want. Set \(\boldsymbol\alpha_k = \boldsymbol\beta_k / \sqrt{n\lambda_k}\).

The lecture briefly renames the component index to avoid a clash with the kernel \(k\). Here \(k\) as a subscript always counts components; \(k(\cdot, \cdot)\) is the kernel.

In Part 4, step 3 rebuilt the directions as \(\mathbf{w}_k = X\boldsymbol\alpha_k\). Here it would be \(\mathbf{w}_k = \phi(X)\boldsymbol\alpha_k\): combine the mapped points with the weights \(\boldsymbol\alpha_k\). That needs \(\phi(\mathbf{x})\), which defeats the purpose. The whole point of going through the kernel matrix was to avoid computing \(\phi\). When \(\phi\) maps into an infinite-dimensional space, as for the RBF kernel, it cannot be computed at all. So we do not compute \(\mathbf{w}_k\), and the eigenvectors of the feature-space covariance matrix stay out of reach.

What we can still compute is the compressed representation. In PCA, compressing a point means storing its projection on each important direction, \(\phi(\mathbf{x}_i)^\top\mathbf{w}_k\) for the \(k\)-th direction. Expand \(\mathbf{w}_k\) and move \(\phi(\mathbf{x}_i)\) inside the sum:

The projection depends only on dot products among the data points, and the kernel computes those. The direction itself stays unknown, yet the value of every projection is computable.

There is a nice reading of this 10:44. The numbers \(K_{i1}, \dots, K_{in}\) are row \(i\) of the kernel matrix: how similar point \(i\) is to each data point, as measured by the kernel. Each new coordinate weights those similarities by one \(\boldsymbol\alpha_k\) and adds them up.

A quick check, not made in the lecture: since \(K\boldsymbol\alpha_k = n\lambda_k\boldsymbol\alpha_k\), the coordinate simplifies to \(\sum_j \alpha_{kj}K_{ij} = n\lambda_k\,\alpha_{ki} = \sqrt{n\lambda_k}\,\beta_{ki}\).

What the projections are for L12 · 7:21

The usual reason to run (kernel) PCA is to compute these projections, throw away the original points, and keep only the top few projections per point. A 100-dimensional point becomes, say, the 5 numbers giving its projections on the 5 most important directions. That is where the dimensionality reduction happens. The projections then serve as the data for a downstream task, typically a supervised learning problem later in the course. For that use, the projections are usually good enough.

What is lost is reconstruction 8:59. A proxy is built by multiplying each projection by its eigenvector and adding up, and that needs \(\mathbf{w}_k\), hence \(\phi\). So kernel PCA cannot rebuild proxies in the feature space. In practice step 3 is simply replaced by the projection formula above.

The dimension can go up L12 · 11:17

The final representation has \(l\) numbers per point, and \(l\) may be larger than the original \(d\). For example, 2-D data might be mapped (implicitly) to a 100-dimensional space and reduced to 10 dimensions there. That looks like an increase from 2 to 10. It is still a low-dimensional representation of the mapped data, and the extra room is what lets it capture the non-linear structure.

One important detail is still missing. PCA centres the data before forming the covariance matrix, because the directions only make sense for centred data (Part 2). Even if the original \(\mathbf{x}_i\) are centred, the mapped points \(\phi(\mathbf{x}_1), \dots, \phi(\mathbf{x}_n)\) may not be. That depends on the map, and nothing about a valid kernel guarantees it. So the data must be centred in the feature space, where we cannot compute anything except kernel values.

The question becomes: can one function do both jobs, mapping with \(\phi\) and then centring? We want a new kernel matrix \(K^c\) whose entries are dot products of the centred mapped points. The good news is that it can be built from \(K\) alone.

There are four terms. The first is \(K_{ij}\). The lecture leaves the rest as routine algebra. Written out, each is an average of kernel entries:

The lecture writes the formula with \(\theta_i\), the average row sum of \(K\) for point \(i\), and indicator vectors \(\mathbf{1}_i\), and does not derive it. The captions do not capture that expression cleanly, so the form above is derived from the four-term expansion the lecture gives; it is the same subtraction.

The important news is that a kernel can be centred without computing \(\phi\). From here on, use \(K^c\) wherever \(K\) appeared: its eigenvectors in step 2 and its rows in the projection step.

1. Kernel PCA
2. Input: \(\{\mathbf{x}_1, \dots, \mathbf{x}_n\}\), \(\mathbf{x}_i \in \mathbb{R}^d\), and a valid kernel \(k : \mathbb{R}^d \times \mathbb{R}^d \to \mathbb{R}\).
3. Kernel matrix. \(K_{ij} = k(\mathbf{x}_i, \mathbf{x}_j)\) for all \(i, j\); \(K\) is \(n \times n\).
4. Centre. Form \(K^c\) from \(K\) with the formula above.
5. Eigen-decompose \(K^c\). Unit eigenvectors \(\boldsymbol\beta_1, \dots, \boldsymbol\beta_l\) with eigenvalues \(n\lambda_1 \ge \dots \ge n\lambda_l\) (the top \(l\), all non-zero).
6. Normalize. \(\boldsymbol\alpha_k = \boldsymbol\beta_k / \sqrt{n\lambda_k}\) for \(k = 1, \dots, l\).
7. Project. Do not build \(\mathbf{w}_k = \phi(X)\boldsymbol\alpha_k\). Represent each point by \(\mathbf{z}_i = \big(\sum_j \alpha_{1j}K^c_{ij}, \dots, \sum_j \alpha_{lj}K^c_{ij}\big)\).
8. Output: \(\mathbf{z}_1, \dots, \mathbf{z}_n \in \mathbb{R}^l\), the kernel PCA representation of the data.

How the three versions of PCA in this block compare:

Kernel PCA settles both issues raised in Part 4. High dimensionality is handled by working with an \(n \times n\) matrix. Non-linearity is handled by mapping to a higher dimension through a kernel. Surprisingly, the same object does both jobs: the kernel matrix.

The trick works because the PCA algorithm depends only on pairwise dot products between data points, not on the points themselves. Given the dot products, we can run it, and a kernel supplies dot products in a high-dimensional space.

To kernelize an algorithm: solve the easier linear problem first, check that the algorithm uses the data only through pairwise dot products, then replace each \(\mathbf{x}_i^\top\mathbf{x}_j\) by \(k(\mathbf{x}_i, \mathbf{x}_j)\).

The kernel matrix comes back again and again in the course. Later algorithms will be checked for the same property and kernelized when they have it: kernel regression in Part 11 and support vector machines in Part 16.

This closes the representation-learning block: PCA and its kernel version, kernel PCA. Next come other kinds of unsupervised learning, starting with clustering in Part 6 and then estimation in Part 8.

Non-linear relationship

Features tied by an equation with powers or products, like a circle; PCA's issue 2.

Feature map \(\phi\)

A map \(\mathbb{R}^d \to \mathbb{R}^D\) that adds non-linear functions of the features as new coordinates.

Feature space

The space \(\mathbb{R}^D\) that \(\phi\) maps into, where the relationship should become linear.

Kernel function

A valid \(k : \mathbb{R}^d\times\mathbb{R}^d \to \mathbb{R}\): \(k(\mathbf{x},\mathbf{x}') = \phi(\mathbf{x})^\top\phi(\mathbf{x}')\) for some \(\phi\).

Valid

Having such a \(\phi\). An explicit \(\phi\) is a certificate of validity.

Kernel matrix \(K\)

The \(n\times n\) matrix \(K_{ij} = k(\mathbf{x}_i, \mathbf{x}_j)\) of pairwise feature-space dot products.

Polynomial kernel

\((\mathbf{x}^\top\mathbf{x}' + 1)^p\): all relationships up to power \(p\), with a finite but fast-growing \(\phi\).

RBF (Gaussian) kernel

\(\exp(-\lVert\mathbf{x}-\mathbf{x}'\rVert^2/(2\sigma^2))\), \(\sigma > 0\); its \(\phi\) is infinite-dimensional.

Positive semi-definite

A symmetric matrix whose eigenvalues are all \(\ge 0\).

Mercer's theorem

\(k\) is a valid kernel exactly when it is symmetric and every kernel matrix it produces is positive semi-definite.

Kernel PCA

PCA in the feature space computed from \(K^c\) alone; outputs projections, not directions.

Centred kernel \(K^c\)

Dot products of mean-subtracted mapped points, computed from \(K\) by subtracting row and column averages and adding back the overall average.

Kernel trick

Replacing every dot product in an algorithm by a kernel value, so it runs in the feature space without \(\phi\).