Part 11 · Kernel and probabilistic regression

At a glance

Kernel regression. A least-squares solution can always be taken in the span of the data, \(\mathbf{w}^* = X\boldsymbol{\alpha}^*\). Then \(\boldsymbol{\alpha}^* = K^{-1}\mathbf{y}\) with the \(n \times n\) kernel matrix \(K\), and a prediction is \(\sum_i \alpha^*_i\, k(\mathbf{x}_i, \mathbf{x}_{\text{test}})\). No feature map is ever computed.
Probabilistic view. Model labels as \(y \mid \mathbf{x} = \mathbf{w}^\top\mathbf{x} + \epsilon\) with \(\epsilon \sim \mathcal{N}(0, \sigma^2)\). Maximum likelihood is then least squares, so \(\hat{\mathbf{w}}_{\text{ML}} = \mathbf{w}^*\). Choosing a noise model means choosing a loss: Laplacian noise gives absolute error.
How good is \(\hat{\mathbf{w}}_{\text{ML}}\)? It is random because \(\mathbf{y}\) is. Its average squared distance from the true \(\mathbf{w}\) is \(\sigma^2\,\mathrm{tr}\big((XX^\top)^{-1}\big) = \sigma^2\sum_i 1/\lambda_i\). The noise is out of our hands, so any improvement has to come from the feature term.
A new estimator. \(\hat{\mathbf{w}}_{\text{new}} = (XX^\top + \lambda I)^{-1}X\mathbf{y}\) adds \(\lambda\) to every eigenvalue, which tames the small ones. Some \(\lambda\) is guaranteed to beat \(\hat{\mathbf{w}}_{\text{ML}}\) in mean squared error, but the theorem does not say which.
Cross-validation. Choose \(\lambda\) by holding data out. K-fold cross-validation averages the validation error over \(K\) held-out folds. Leave-one-out is the extreme case \(K = n\).

L32 · 0:00 L32 · 1:03 the kernel PCA idea

When the relationship is not linear

So far the label has been explained by \(\mathbf{w}^\top\mathbf{x}\). Suppose the features are height and weight and the label is best explained by something like height² + 3 × weight − weight², or needs a cross term such as height × weight. Linear regression can still be run, but no \(\mathbf{w}\) can produce squared or cross terms. Even the best \(\mathbf{w}\) will then make large errors on the training set and on the validation set.

The way out is the one that took us from PCA to kernel PCA in Part 5. There, the PCA problem was rewritten so that the data entered only through dot products between data points, and those dot products were then replaced by a kernel matrix. The effect is the same as mapping every point to a higher-dimensional space and working there, without ever computing the mapping. Here we want a linear regressor in that higher-dimensional space.

To kernelize regression, express its solution through the dot products \(\mathbf{x}_i^\top\mathbf{x}_j\) alone.

L32 · 2:06 L32 · 4:51 the argument

The solution lies in the span of the data

The closed form from Part 10 is \(\mathbf{w}^* = (XX^\top)^{\dagger}X\mathbf{y}\), where \(X = [\mathbf{x}_1 \cdots \mathbf{x}_n] \in \mathbb{R}^{d\times n}\) holds the data points as columns and \(\mathbf{y} \in \mathbb{R}^n\) the labels. The lecture writes "inverse" throughout and means the pseudo-inverse whenever the inverse does not exist.

As written, this does not depend on dot products directly. \(XX^\top\) is a \(d \times d\) matrix. The matrix of pairwise dot products, \(X^\top X\), is \(n \times n\), and it does not appear. The step that unlocks a rewrite is this observation:

Key observation

A minimizer \(\mathbf{w}^*\) of the squared error can always be taken to lie in the span of the data points.

The lecture argues it with a picture 3:42. Take 3-D data points that all lie on a 2-D plane through the origin. Suppose a minimizer \(\mathbf{w}^*\) sticks out of the plane. Let \(\tilde{\mathbf{w}}\) be its projection onto the plane, so that \(\mathbf{w}^* = \tilde{\mathbf{w}} + \tilde{\mathbf{w}}_\perp\), where \(\tilde{\mathbf{w}}_\perp\) is perpendicular to the plane and therefore to every data point:

\[ \mathbf{w}^{*\top}\mathbf{x}_i = (\tilde{\mathbf{w}} + \tilde{\mathbf{w}}_\perp)^\top\mathbf{x}_i = \tilde{\mathbf{w}}^\top\mathbf{x}_i + \underbrace{\tilde{\mathbf{w}}_\perp^\top\mathbf{x}_i}_{=\,0} = \tilde{\mathbf{w}}^\top\mathbf{x}_i \qquad \text{for every } i \]

So \(\tilde{\mathbf{w}}\) makes exactly the same prediction as \(\mathbf{w}^*\) on every training point, and suffers exactly the same squared error \(\sum_i(\tilde{\mathbf{w}}^\top\mathbf{x}_i - y_i)^2\). It cannot be strictly better, since \(\mathbf{w}^*\) is a minimizer, and it is not worse. We may as well use \(\tilde{\mathbf{w}}\). In general, without loss of generality, the minimizer is a linear combination of the data points:

\[ \mathbf{w}^* = X\boldsymbol{\alpha}^* = \sum_{i=1}^n \alpha^*_i\,\mathbf{x}_i \qquad \text{for some } \boldsymbol{\alpha}^* \in \mathbb{R}^n \]

\(\boldsymbol{\alpha}^*\) has one entry per data point: \(\alpha^*_i\) is the weight given to \(\mathbf{x}_i\) when the data points are added up to form \(\mathbf{w}^*\).

Leaving the span gains nothing. The part of \(\mathbf{w}^*\) perpendicular to the data's plane is orthogonal to every \(\mathbf{x}_i\), so dropping it changes no prediction on the training data and no error. A minimizer inside the span always exists.

L32 · 9:04 L32 · 11:21 the kernel matrix

Solving for \(\boldsymbol{\alpha}^*\)

We now have two descriptions of the same vector: \(\mathbf{w}^* = X\boldsymbol{\alpha}^*\), and the closed form from the gradient analysis of Part 10. Setting them equal:

\[ X\boldsymbol{\alpha}^* = (XX^\top)^{\dagger}X\mathbf{y} \]

The goal is to get rid of \(XX^\top\) and see \(X^\top X\) instead. First pre-multiply both sides by \(XX^\top\), which cancels the (pseudo-)inverse on the right:

\[ XX^\top X\boldsymbol{\alpha}^* = XX^\top(XX^\top)^{\dagger}X\mathbf{y} = X\mathbf{y} \]

Strictly, \(XX^\top(XX^\top)^{\dagger}\) is the projection onto the column space of \(X\), and \(X\mathbf{y}\) already lies in that space, so the right side is exactly \(X\mathbf{y}\). The lecture simply says the two cancel.

The left side now contains \(X^\top X\), but an \(X\) also sits alone on the right. Pre-multiply both sides by \(X^\top\):

\[ X^\top X\,X^\top X\,\boldsymbol{\alpha}^* = X^\top X\,\mathbf{y} \quad\Longleftrightarrow\quad (X^\top X)^2\,\boldsymbol{\alpha}^* = X^\top X\,\mathbf{y} \]

\(X^\top X\) is the matrix of all pairwise dot products, \((X^\top X)_{ij} = \mathbf{x}_i^\top\mathbf{x}_j\): the kernel matrix \(K\) we met in kernel PCA. In terms of \(K\):

Result \[ K^2\boldsymbol{\alpha}^* = K\mathbf{y}, \qquad \text{so if } K \text{ is invertible,} \qquad \boldsymbol{\alpha}^* = K^{-1}\mathbf{y}, \qquad K = X^\top X \in \mathbb{R}^{n\times n} \]

The weights that combine the data points into \(\mathbf{w}^*\) are the inverse of the kernel matrix applied to the labels. The data now enters only through dot products.

L32 · 11:55 L32 · 14:49 the algorithm L32 · 15:20 reading it

Predicting with a kernel

What is \(\boldsymbol{\alpha}^*\) good for? Map every point to a higher-dimensional space with a feature map \(\phi\) and solve linear regression there:

\[ \mathbf{w}^* = \arg\min_{\mathbf{w}} \sum_{i=1}^n \big(\mathbf{w}^\top\phi(\mathbf{x}_i) - y_i\big)^2, \qquad \text{prediction for a new point: } \mathbf{w}^{*\top}\phi(\mathbf{x}_{\text{test}}) \]

Computing \(\phi\) explicitly is exactly what we cannot afford; that is the whole point of kernels (Part 5). But the argument above works unchanged with \(\phi(\mathbf{x}_i)\) in place of \(\mathbf{x}_i\). So \(\mathbf{w}^*\) is a linear combination of the mapped data points, and the kernel matrix becomes \(K_{ij} = \phi(\mathbf{x}_i)^\top\phi(\mathbf{x}_j) = k(\mathbf{x}_i, \mathbf{x}_j)\):

\[ \mathbf{w}^* = \sum_{i=1}^n \alpha^*_i\,\phi(\mathbf{x}_i) \quad\Longrightarrow\quad \mathbf{w}^{*\top}\phi(\mathbf{x}_{\text{test}}) = \sum_{i=1}^n \alpha^*_i\,\phi(\mathbf{x}_i)^\top\phi(\mathbf{x}_{\text{test}}) = \sum_{i=1}^n \alpha^*_i\,k(\mathbf{x}_i, \mathbf{x}_{\text{test}}) \]

The kernel function evaluates dot products in the high-dimensional space without going there. \(K\) is the same \(n \times n\) matrix as in kernel PCA, evaluated on the data set. The one new ingredient is \(\mathbf{y}\), through \(\boldsymbol{\alpha}^* = K^{-1}\mathbf{y}\). It has to appear: PCA is unsupervised, while regression is supervised, so the labels must matter.

Same prediction, no feature map. The greyed-out route needs \(\phi\) and a vector \(\mathbf{w}^*\) in the huge space. Kernel regression stores only \(n\) numbers \(\boldsymbol{\alpha}^*\) and needs only kernel evaluations, and it gives the same \(\hat y\) because \(\mathbf{w}^* = \sum_i \alpha^*_i\phi(\mathbf{x}_i)\).

Kernel regression
Input: data \(\{(\mathbf{x}_1, y_1), \dots, (\mathbf{x}_n, y_n)\}\) and a kernel function \(k\).
Evaluate the kernel on every pair of data points to get the \(n \times n\) matrix \(K\), \(K_{ij} = k(\mathbf{x}_i, \mathbf{x}_j)\).
Compute \(\boldsymbol{\alpha}^* = K^{-1}\mathbf{y}\).
For any test point, predict \(\hat y = \sum_{i=1}^n \alpha^*_i\, k(\mathbf{x}_i, \mathbf{x}_{\text{test}})\).

Reading the prediction

Think of the kernel as a measure of similarity. \(k(\mathbf{x}_i, \mathbf{x}_{\text{test}})\) says how similar the test point is to training point \(i\). \(\alpha^*_i\) says how important point \(i\) is to \(\mathbf{w}^*\). A test point can be very similar to a training point that matters little for \(\mathbf{w}^*\), so neither number decides alone. The prediction adds up their product over all training points.

Picking a kernel (quadratic, cubic, Gaussian, which is the RBF kernel of Part 5, or any other) picks the higher-dimensional space in which the regression is solved 16:30. This is kernel regression: regression is to kernel regression as PCA is to kernel PCA.

Kernel linear, 1 + xx′ Test point x_test 0.35

Prediction ŷ at x_test–

Training error (MSE)–

Feature-space dimension–

Kernel regression on synthetic data (10 points generated for this page from \(y = 1.6x^2 - 0.6\) plus noise). Every curve comes from the same three lines of the algorithm: build \(K\), solve for \(\boldsymbol{\alpha}^*\), add up \(\alpha^*_i k(\mathbf{x}_i, x_{\text{test}})\). The bars are those products (rescaled so the tallest fills the panel) and add up to \(\hat y\). The linear kernel can only fit a line; the polynomial kernels bend; the Gaussian kernels make only nearby points count. For the polynomial kernels \(K\) is singular (10 points, at most 4 feature dimensions), so the widget uses the pseudo-inverse \(K^{\dagger}\), which still solves \(K^2\boldsymbol{\alpha} = K\mathbf{y}\).

L32 · 17:06

Three views so far, and the next one

The lecture takes stock of what we know about linear regression:

Geometry. The predictions \(X^\top\mathbf{w}^*\) are the projection of the label vector onto the subspace spanned by the features (Part 10).
Computation. Gradient descent and stochastic gradient descent avoid the expensive inverse (Part 10).
Non-linear relationships. Kernel regression, from this lecture.

Everything so far is deterministic. The next aspect brings in probability, and it leads to several interesting variants of linear regression 18:14.

L33 · 0:00 L33 · 1:05 the model L33 · 2:49 Gaussian noise

A probabilistic model for the labels

In Part 8 we estimated parameters in an unsupervised setting with maximum likelihood and Bayesian methods. Now linear regression itself is treated as an estimation problem. That needs a probabilistic mechanism that we assume generated the data we see.

The data set is \(\{(\mathbf{x}_1, y_1), \dots, (\mathbf{x}_n, y_n)\}\) with \(\mathbf{x}_i \in \mathbb{R}^d\) and \(y_i \in \mathbb{R}\). The model:

\[ y \mid \mathbf{x} \;=\; \underbrace{\mathbf{w}^\top\mathbf{x}}_{\text{structure}} \;+\; \underbrace{\epsilon}_{\text{noise}}, \qquad \epsilon \sim \mathcal{N}(0, \sigma^2) \]

We do not model how the features are generated, only how a label arises once a feature vector is given.
\(\mathbf{w} \in \mathbb{R}^d\) is unknown but fixed. It is the parameter of the problem, and \(\mathbf{w}^\top\mathbf{x}\) is the structure part.
The label is not exactly \(\mathbf{w}^\top\mathbf{x}\): noise \(\epsilon\) is added explicitly. It is Gaussian with mean 0 and a known variance \(\sigma^2\).

So each \(y_i\) arose like this: someone gave us \(\mathbf{x}_i\), the unknown \(\mathbf{w}\) produced \(\mathbf{w}^\top\mathbf{x}_i\), noise was added, and we see only the noisy result. The noise statistics are known and only \(\mathbf{w}\) is not, so finding \(\mathbf{w}\) is an estimation problem 3:54. The simplest method we know is maximum likelihood.

Structure plus noise. Synthetic points (seeded, generated for this page) from \(y = 0.7x + \epsilon\) with \(\sigma = 0.7\). At every \(x\) the label is a Gaussian centred on the line \(\mathbf{w}^\top\mathbf{x}\) with the same spread \(\sigma\); each point's vertical gap to the line is its noise \(\epsilon_i\).

L33 · 4:28 L33 · 7:06 log-likelihood L33 · 8:39 least squares

Maximum likelihood gives back least squares

The likelihood is a function of the parameter \(\mathbf{w}\). It also depends on the observed data, but we treat the data as fixed and look for the \(\mathbf{w}\) that makes it most likely. The usual i.i.d. assumption holds: each \(y_i\) is generated independently given its \(\mathbf{x}_i\), with the same noise distribution. So the likelihood is a product:

\[ \mathcal{L}(\mathbf{w};\ \mathbf{x}_1, \dots, \mathbf{x}_n, y_1, \dots, y_n) = \prod_{i=1}^n p(y_i \mid \mathbf{x}_i;\ \mathbf{w}) \]

What is the density of \(y_i\) given \(\mathbf{x}_i\)? \(\mathbf{w}^\top\mathbf{x}_i\) is a fixed number, nothing random. Adding a constant to a zero-mean Gaussian moves the mean and leaves the spread unchanged: add 5 to a zero-mean Gaussian and you get a Gaussian with mean 5 and the same variance. So \(y_i \mid \mathbf{x}_i \sim \mathcal{N}(\mathbf{w}^\top\mathbf{x}_i, \sigma^2)\) 6:04, and

\[ \mathcal{L}(\mathbf{w}) = \prod_{i=1}^n \frac{1}{\sqrt{2\pi}\,\sigma}\exp\!\left(-\frac{(\mathbf{w}^\top\mathbf{x}_i - y_i)^2}{2\sigma^2}\right) \]

Take the logarithm, since sums are easier to handle than products. The log cancels the exponential:

\[ \log\mathcal{L}(\mathbf{w}) = \sum_{i=1}^n \left[ -\frac{(\mathbf{w}^\top\mathbf{x}_i - y_i)^2}{2\sigma^2} + \log\frac{1}{\sqrt{2\pi}\,\sigma} \right] \]

Only \(\mathbf{w}\) varies; the \(\mathbf{x}_i\), \(y_i\) and the known \(\sigma\) are constants. The second term does not involve \(\mathbf{w}\), and \(1/(2\sigma^2)\) is a positive scaling. Dropping both does not move the maximizer:

\[ \max_{\mathbf{w}} \log\mathcal{L}(\mathbf{w}) \;\Longleftrightarrow\; \max_{\mathbf{w}} \sum_{i=1}^n -(\mathbf{w}^\top\mathbf{x}_i - y_i)^2 \;\Longleftrightarrow\; \min_{\mathbf{w}} \sum_{i=1}^n (\mathbf{w}^\top\mathbf{x}_i - y_i)^2 \]

The last problem is exactly linear regression with squared error, whose solution we already know 9:15:

Result \[ \hat{\mathbf{w}}_{\text{ML}} = \mathbf{w}^* = (XX^\top)^{\dagger}X\mathbf{y} \]

A completely different starting point (a generative story for \(y\) given \(\mathbf{x}\), then the most natural estimator) lands on the same answer as least squares.

L33 · 9:49 L33 · 12:08 other noise

Choosing a noise model is choosing a loss

Conclusion

The maximum likelihood estimator under zero-mean Gaussian noise is exactly linear regression with squared error.

Both qualifiers matter. When squared error was introduced in Part 10 it was justified mainly as the intuitive thing to do. Now it has a probabilistic backing. Choosing squared error is implicitly assuming that zero-mean Gaussian noise corrupts the labels.

The converse holds too. If there is reason to believe the noise is not Gaussian, least squares is no longer the maximum likelihood answer. The density assumed for the noise determines the error function.

Example: Laplacian noise L33 · 12:39

The lecture mentions this case only for completeness; here are the steps. The Laplace density has an absolute value in the exponent instead of a square. With a scale \(b > 0\) (the lecture writes it without one):

\[ p(y_i \mid \mathbf{x}_i;\ \mathbf{w}) = \frac{1}{2b}\exp\!\left(-\frac{|\mathbf{w}^\top\mathbf{x}_i - y_i|}{b}\right) \quad\Longrightarrow\quad \log\mathcal{L}(\mathbf{w}) = -\frac1b\sum_{i=1}^n |\mathbf{w}^\top\mathbf{x}_i - y_i| + \text{const} \]

So maximum likelihood under Laplacian noise minimizes \(\sum_i |\mathbf{w}^\top\mathbf{x}_i - y_i|\), the sum of absolute errors, and the two problems are no longer the same.

Noise model in, loss function out. Taking the negative log of the noise density gives the error each residual \(\epsilon = y_i - \mathbf{w}^\top\mathbf{x}_i\) costs. The Gaussian's square in the exponent becomes the squared error; the Laplace density's absolute value becomes the absolute error.

Choosing a noise model implicitly chooses an error function, and choosing an error function implicitly chooses a noise model.

What else does the probabilistic view buy? L33 · 13:45

The link between noise and loss is one gain. The bigger one is that \(\hat{\mathbf{w}}_{\text{ML}}\) is now an estimator. Once learning is viewed as estimation, all the machinery for judging estimators applies, and we can ask what makes an estimator good. Lecture 34 does exactly that.

L34 · 1:08 L34 · 2:48 why it is random

\(\hat{\mathbf{w}}_{\text{ML}}\) is a random variable

The setting is unchanged: there is a true \(\mathbf{w}\) with \(y \mid \mathbf{x} \sim \mathcal{N}(\mathbf{w}^\top\mathbf{x}, \sigma^2)\), and \(\hat{\mathbf{w}}_{\text{ML}} \in \mathbb{R}^d\) is our guess for it. One property of this estimator will suggest a way to extend linear regression. The question is how good a guess it is.

\(\hat{\mathbf{w}}_{\text{ML}}\) is the solution of an optimization problem built from the data \(\mathbf{x}\) and \(\mathbf{y}\). Nothing is assumed about how \(\mathbf{x}\) arises, but \(\mathbf{y}\) has a random component, the noise. So \(\hat{\mathbf{w}}_{\text{ML}}\), a function of \(\mathbf{y}\), is random as well.

That is why one data set cannot settle the question 3:52. On an unlucky day the noise values are large (\(y_1 = \mathbf{w}^\top\mathbf{x}_1 + \epsilon_1\) with \(\epsilon_1\) far from 0, and so on), the labels carry little information about \(\mathbf{w}\), and the estimate is poor. So we ask how \(\hat{\mathbf{w}}_{\text{ML}}\) performs on average.

L34 · 5:02 L34 · 6:13 the thought experiment L34 · 7:54 the value

Average squared distance from the truth

Both \(\hat{\mathbf{w}}_{\text{ML}}\) and \(\mathbf{w}\) are vectors in \(\mathbb{R}^d\), so measure one estimate's error by the squared Euclidean distance \(\lVert\hat{\mathbf{w}}_{\text{ML}} - \mathbf{w}\rVert^2\), and average over the randomness in \(\mathbf{y}\):

\[ \mathbb{E}\Big[\,\lVert \hat{\mathbf{w}}_{\text{ML}} - \mathbf{w} \rVert^2\,\Big] \]

One way to picture the expectation: create a data set whose labels are drawn as \(\mathbf{w}^\top\mathbf{x}_i\) plus fresh zero-mean Gaussian noise, compute \(\hat{\mathbf{w}}_{\text{ML}}\), and record its squared distance to \(\mathbf{w}\). Repeat many times. The average of the recorded distances converges to the expectation: the average deviation over all data sets you could have received.

The lecture does not derive the value. It notes that, given \(\hat{\mathbf{w}}_{\text{ML}} = (XX^\top)^{-1}X\mathbf{y}\) and the properties of the random \(\mathbf{y}\), it is a short algebraic exercise, and states the answer:

Result \[ \mathbb{E}\Big[\,\lVert \hat{\mathbf{w}}_{\text{ML}} - \mathbf{w} \rVert^2\,\Big] = \sigma^2\;\mathrm{tr}\big((XX^\top)^{-1}\big) \]

Use the pseudo-inverse if \(XX^\top\) is singular. This average squared deviation is the mean squared error (MSE) of the estimator, the name L35 uses.

The derivation (the lecture leaves it as an exercise)

Assume \((XX^\top)^{-1}\) exists. Stack the labels: \(\mathbf{y} = X^\top\mathbf{w} + \boldsymbol{\epsilon}\), where \(\boldsymbol{\epsilon} \in \mathbb{R}^n\) has independent \(\mathcal{N}(0, \sigma^2)\) entries, so \(\mathbb{E}[\boldsymbol{\epsilon}] = \mathbf{0}\) and \(\mathbb{E}[\boldsymbol{\epsilon}\boldsymbol{\epsilon}^\top] = \sigma^2 I\). Then

\[ \hat{\mathbf{w}}_{\text{ML}} = (XX^\top)^{-1}X(X^\top\mathbf{w} + \boldsymbol{\epsilon}) = \mathbf{w} + A\boldsymbol{\epsilon}, \qquad A = (XX^\top)^{-1}X \]

\[ \mathbb{E}\lVert A\boldsymbol{\epsilon}\rVert^2 = \mathbb{E}\big[\boldsymbol{\epsilon}^\top A^\top A\,\boldsymbol{\epsilon}\big] = \mathrm{tr}\big(A^\top A\;\mathbb{E}[\boldsymbol{\epsilon}\boldsymbol{\epsilon}^\top]\big) = \sigma^2\,\mathrm{tr}(AA^\top) = \sigma^2\,\mathrm{tr}\big((XX^\top)^{-1}XX^\top(XX^\top)^{-1}\big) = \sigma^2\,\mathrm{tr}\big((XX^\top)^{-1}\big) \]

A by-product of the first line: \(\mathbb{E}[\hat{\mathbf{w}}_{\text{ML}}] = \mathbf{w}\), so on average the estimate is centred on the truth. Standard name: \(\hat{\mathbf{w}}_{\text{ML}}\) is unbiased. Its whole error is spread around \(\mathbf{w}\), which is what the formula measures.

Spread out · σ² tr((XXᵀ)⁻¹)–

Spread out · simulated–

Collinear · σ² tr((XXᵀ)⁻¹)–

Collinear · simulated–

The thought experiment, run 200 times on synthetic data (seeded, generated for this page). True \(\mathbf{w} = (1, 1)\), \(\sigma = 1\), \(n = 20\); the features stay fixed and only the noise is redrawn each time. Each dot is one \(\hat{\mathbf{w}}_{\text{ML}}\). The noise level is the same in both panels, but with nearly collinear features one eigenvalue of \(XX^\top\) is small and the estimates scatter far more. The simulated averages of \(\lVert\hat{\mathbf{w}}_{\text{ML}} - \mathbf{w}\rVert^2\) match the formula.

L34 · 8:24 L34 · 10:34 what we control

Two factors, and only one we can touch

The mean squared error is a product of two terms.

\(\sigma^2\), the noise variance. More noise confuses \(\mathbf{w}^\top\mathbf{x}\) more, so the estimate gets worse. That is expected.
\(\mathrm{tr}\big((XX^\top)^{-1}\big)\), a property of the features alone. \(XX^\top\) behaves like a covariance matrix, so this is like the trace of an inverse covariance. How the features relate to each other also decides how good the estimator is. That is the interesting part.

\(XX^\top = \sum_i \mathbf{x}_i\mathbf{x}_i^\top\) is \(n\) times the matrix \(C\) of Part 2, which is the covariance matrix when the data is centred.

The data-generating process is not ours to change, so \(\sigma^2\) is out of our hands: it is noise we have to suffer. The features, though, are given to us, and they enter through the trace term. Could we tweak something there and get an estimator with a smaller mean squared error? It is not obvious that this is possible. If it is, it has to work on the trace term 11:07.

L35 · 0:13 L35 · 2:31 eigenvalues of the inverse

The trace as a sum of eigenvalues

The trace of a square matrix is the sum of its diagonal entries. It is also the sum of its eigenvalues, a standard fact (it follows from the characteristic polynomial) that the lecture uses without proof:

\[ \mathrm{tr}(A) = \sum_{i=1}^d a_{ii} = \sum_{i=1}^d \lambda_i, \qquad \lambda_i = \text{the } i\text{th eigenvalue of } A \]

Let \(\lambda_1, \dots, \lambda_d\) be the eigenvalues of \(XX^\top\). The lecture pauses here with a question: what are the eigenvalues of \((XX^\top)^{-1}\)? They are the reciprocals \(1/\lambda_1, \dots, 1/\lambda_d\). (If \(XX^\top\mathbf{v} = \lambda\mathbf{v}\), multiply both sides by \((XX^\top)^{-1}\) and divide by \(\lambda\) to get \((XX^\top)^{-1}\mathbf{v} = \tfrac{1}{\lambda}\mathbf{v}\).) Therefore

Mean squared error of the ML estimator \[ \mathrm{MSE}(\hat{\mathbf{w}}_{\text{ML}}) = \mathbb{E}\lVert\hat{\mathbf{w}}_{\text{ML}} - \mathbf{w}\rVert^2 = \sigma^2 \sum_{i=1}^d \frac{1}{\lambda_i}, \qquad \lambda_i = \text{eigenvalues of } XX^\top \]

A single small eigenvalue has a huge reciprocal and drags the whole error up with it. That is what happened in the collinear panel above.

L35 · 3:38 L35 · 5:27 shifting eigenvalues L35 · 8:18 the intuition

A new estimator: add \(\lambda I\)

The lecturer says he will "pull a new estimator out of the hat" and postpones where it comes from. It is not the maximum likelihood estimator:

\[ \hat{\mathbf{w}}_{\text{new}} = (XX^\top + \lambda I)^{-1}X\mathbf{y}, \qquad \lambda > 0, \quad I = d\times d \text{ identity} \]

Compared with \(\hat{\mathbf{w}}_{\text{ML}} = (XX^\top)^{-1}X\mathbf{y}\), the only change is the extra \(\lambda I\) inside the inverse. The lecture uses \(\lambda\) both for this number and for eigenvalues; the eigenvalues always carry a subscript.

Adding \(\lambda I\) shifts every eigenvalue by \(\lambda\)

Take any matrix \(A\) with eigenvalues \(\lambda_1, \dots, \lambda_d\), and an eigenvector \(\mathbf{v}_1\) with \(A\mathbf{v}_1 = \lambda_1\mathbf{v}_1\). Then

\[ (A + \lambda I)\,\mathbf{v}_1 = A\mathbf{v}_1 + \lambda\mathbf{v}_1 = \lambda_1\mathbf{v}_1 + \lambda\mathbf{v}_1 = (\lambda_1 + \lambda)\,\mathbf{v}_1 \]

So \(\mathbf{v}_1\) is still an eigenvector, now with eigenvalue \(\lambda_1 + \lambda\). The eigenvalues of \(A + \lambda I\) are \(\lambda_1 + \lambda, \dots, \lambda_d + \lambda\). With \(A = XX^\top\):

\[ \mathrm{tr}\big((XX^\top + \lambda I)^{-1}\big) = \sum_{i=1}^d \frac{1}{\lambda_i + \lambda} \;<\; \sum_{i=1}^d \frac{1}{\lambda_i} = \mathrm{tr}\big((XX^\top)^{-1}\big) \]

Why this might help

The MSE of \(\hat{\mathbf{w}}_{\text{ML}}\) can be large for two reasons: \(\sigma^2\) is large, which we cannot fix, or \(\mathrm{tr}\big((XX^\top)^{-1}\big)\) is large, which means some eigenvalues \(\lambda_i\) are small and their reciprocals large. The new estimator should be governed by \(\mathrm{tr}\big((XX^\top + \lambda I)^{-1}\big)\) instead. Each denominator has grown by \(\lambda\), so each term has shrunk, and the small eigenvalues shrink the most. The hope is that the overall MSE goes down.

Eigenvalue of \(XX^\top\)	\(1/\lambda_i\)	\(1/(\lambda_i + 1)\)
\(\lambda_1 = 10\)	0.100	0.091
\(\lambda_2 = 2\)	0.500	0.333
\(\lambda_3 = 0.05\)	20.000	0.952
Trace (sum)	20.600	1.377

Illustrative eigenvalues, with \(\lambda = 1\). Almost all of the original trace comes from the one small eigenvalue, and that is the term the shift tames.

The lecture stresses that this is only an intuitive argument for why \(\lambda I\) might lower the MSE. (Standard context, not covered in the lecture: \(\hat{\mathbf{w}}_{\text{new}}\) is shrunk toward \(\mathbf{0}\), so on average it no longer equals \(\mathbf{w}\). Its MSE gains a bias part that grows with \(\lambda\), which is why \(\lambda\) cannot simply be made huge.)

L35 · 9:56 L35 · 11:19 but which λ?

Some \(\lambda\) is always better, but which?

Does the hope always come true? The lecture states an existence theorem, informally and without proof:

Existence theorem (informal) \[ \exists\,\lambda > 0 \;\text{ such that }\; \hat{\mathbf{w}}_{\text{new}} = (XX^\top + \lambda I)^{-1}X\mathbf{y} \;\text{ satisfies }\; \mathbb{E}\lVert\hat{\mathbf{w}}_{\text{new}} - \mathbf{w}\rVert^2 < \mathbb{E}\lVert\hat{\mathbf{w}}_{\text{ML}} - \mathbf{w}\rVert^2 \]

For a suitable \(\lambda\), the new estimator is better on average than the maximum likelihood estimator.

The catch is in the word "exists". The theorem says a good \(\lambda\) is out there; it does not say how to find it. The right \(\lambda\) depends on the data, so in general there may be no automatic formula for it. In practice \(\lambda\) is found by a procedure called cross-validation, which comes up again and again in the rest of the course.

L35 · 11:50 L35 · 14:08 why "validation"

Validation: hold data back to choose \(\lambda\)

We do not know which \(\lambda\) performs well, so we measure it. Say the training data has 10,000 points. Split it into a training set (say 80%) and a validation set (the other 20%).

Choosing λ with a validation set
Split the data into a training set and a validation set.
For each candidate \(\lambda\), compute \(\hat{\mathbf{w}}_{\text{new}} = (XX^\top + \lambda I)^{-1}X\mathbf{y}\) using the training set only.
Predict the label of every validation point and compare with its actual \(y\) to get the validation error. (With finitely many possible labels one would speak of accuracy; for regression it is an error.)
Pick the \(\lambda\) with the least validation error.

The name comes from the role of the held-out part: it is a proxy for test data. The goal is to do well on test data we will meet later but do not have now, so we mimic it with a part of the training data that we never train on.

L35 · 14:41 L35 · 15:23 the folds L35 · 17:11 choosing K

K-fold cross-validation

Which 80% should go to training and which 20% to validation? A random split is possible, but it can be unlucky: the validation 20% may not represent what was learnt from the training 80%, and then the error measured for a \(\lambda\) does not reflect its true error. To avoid this, the validation is repeated over several splits.

K-fold cross-validation
Divide the data into \(K\) equal parts, the folds \(F_1, \dots, F_K\). With 10,000 points and \(K = 10\), each fold has 1,000 points.
For every candidate \(\lambda\) and every \(i = 1, \dots, K\): train on all folds except \(F_i\), that is on \(F_1, \dots, F_{i-1}, F_{i+1}, \dots, F_K\), and compute the error on \(F_i\).
Average the \(K\) errors for each \(\lambda\).
Pick the \(\lambda\) with the least average error.

The lecture's example with \(K = 3\) 15:59: train on \(F_1, F_2\) and validate on \(F_3\); train on \(F_1, F_3\) and validate on \(F_2\); train on \(F_2, F_3\) and validate on \(F_1\). (This \(K\), the number of folds, has nothing to do with the kernel matrix \(K\) of L32.)

One split versus K splits. A single validation set gives one, possibly unrepresentative, error per \(\lambda\). K-fold cross-validation trains \(K\) times per \(\lambda\), each time leaving out a different fold, and averages the \(K\) validation errors before choosing.

How large should K be?

In general, the larger \(K\) the better, because the average runs over more possibilities. The extreme is leave-one-out cross-validation: with \(n\) points, train on \(n - 1\) of them and check the model on the one left out, for each of the \(n\) points in turn. With a billion data points that means a billion train/validation splits, which is computationally a nightmare. In practice \(K\) is chosen according to the time and resources available 17:48. The hope is that the \(\lambda\) found this way also does well on the test data with the new estimator.

λ (log scale) λ = 1.0×10⁻⁵ Folds K = 5

This fold's validation error–

K-fold CV error (average)–

Training error (average)–

Best λ for this K–

Choosing \(\lambda\) for \(\hat{\mathbf{w}}_{\text{new}}\) by K-fold cross-validation, on synthetic data (20 points generated for this page from \(y = \sin(\pi x)\) plus noise). Each point's features are \(1, x, x^2, \dots, x^9\) (\(d = 10\)), so with tiny \(\lambda\) the estimator chases the noise. Left: the current held-out fold (highlighted) and the fit trained on the other folds. Right: both errors averaged over the \(K\) splits. Training error only falls as \(\lambda\) shrinks; validation error is lowest in between, and that is the \(\lambda\) cross-validation picks. \(K = 20\) is leave-one-out.

L35 · 18:19 L35 · 19:31 the hint

Is there a principled route to \(\hat{\mathbf{w}}_{\text{new}}\)?

So far \(\hat{\mathbf{w}}_{\text{new}}\) has been justified by the eigenvalue argument and the existence theorem. One question lingers: how could we have started thinking about the problem so that this estimator came out naturally?

The lecture says the answer is yes, and that the tool is one we have already used to get alternative estimators, though only in an unsupervised setting. It leaves the guess to the student: which alternative to maximum likelihood did we use for probabilistic models? In Part 8 that alternative was Bayesian estimation. Part 12 follows this up with Bayesian modelling for linear regression, which leads to \(\hat{\mathbf{w}}_{\text{new}}\), known there as ridge regression.

Key terms

Span of the data: All linear combinations of the data points. A least-squares minimizer can always be taken there: \(\mathbf{w}^* = X\boldsymbol{\alpha}^*\).
Coefficients \(\boldsymbol{\alpha}^*\): One weight per training point; \(\alpha^*_i\) says how much \(\mathbf{x}_i\) (or \(\phi(\mathbf{x}_i)\)) contributes to \(\mathbf{w}^*\).
Kernel matrix \(K\): The \(n\times n\) matrix \(K_{ij} = k(\mathbf{x}_i, \mathbf{x}_j)\); \(X^\top X\) for the linear kernel.
Kernel regression: \(\boldsymbol{\alpha}^* = K^{-1}\mathbf{y}\), then predict \(\sum_i\alpha^*_i k(\mathbf{x}_i, \mathbf{x}_{\text{test}})\): linear regression in feature space without computing \(\phi\).
Probabilistic model of labels: \(y \mid \mathbf{x} = \mathbf{w}^\top\mathbf{x} + \epsilon\): an unknown fixed \(\mathbf{w}\) (structure) plus noise.
Gaussian noise: \(\epsilon \sim \mathcal{N}(0, \sigma^2)\) with \(\sigma^2\) known; then \(y \mid \mathbf{x} \sim \mathcal{N}(\mathbf{w}^\top\mathbf{x}, \sigma^2)\).
\(\hat{\mathbf{w}}_{\text{ML}}\): The maximum likelihood estimate of \(\mathbf{w}\); under Gaussian noise it equals the least-squares \(\mathbf{w}^* = (XX^\top)^{\dagger}X\mathbf{y}\).
Laplacian noise: Noise with density \(\propto e^{-|\epsilon|/b}\); maximum likelihood then minimizes the sum of absolute errors.
Mean squared error of an estimator: \(\mathbb{E}\lVert\hat{\mathbf{w}} - \mathbf{w}\rVert^2\), averaged over the noise; \(\sigma^2\,\mathrm{tr}((XX^\top)^{-1})\) for \(\hat{\mathbf{w}}_{\text{ML}}\).
Trace: Sum of the diagonal entries of a square matrix, which equals the sum of its eigenvalues.
\(\hat{\mathbf{w}}_{\text{new}}\): \((XX^\top + \lambda I)^{-1}X\mathbf{y}\): adds \(\lambda\) to every eigenvalue of \(XX^\top\). Called ridge regression in Part 12.
Validation set: Training data held out from fitting and used as a stand-in for test data to compare choices such as \(\lambda\).
K-fold cross-validation: Split into \(K\) folds, validate on each fold after training on the rest, average the \(K\) errors, pick the best \(\lambda\).
Leave-one-out cross-validation: K-fold with \(K = n\): each point is held out alone. Accurate but needs \(n\) training runs per \(\lambda\).

Test yourself

Q1Why can we assume the least-squares solution \(\mathbf{w}^*\) is a linear combination of the data points?

Split any minimizer into its projection \(\tilde{\mathbf{w}}\) onto the span of the data and a part \(\tilde{\mathbf{w}}_\perp\) orthogonal to every \(\mathbf{x}_i\). Then \(\mathbf{w}^{*\top}\mathbf{x}_i = \tilde{\mathbf{w}}^\top\mathbf{x}_i\) for all \(i\): same predictions, same error. So \(\tilde{\mathbf{w}} = X\boldsymbol{\alpha}\) is also a minimizer.

Q2You have \(n = 500\) points with \(d = 10{,}000\) features. What are the sizes of \(XX^\top\) and of the kernel matrix \(K\)?

\(XX^\top\) is \(10{,}000 \times 10{,}000\). \(K = X^\top X\) is \(500 \times 500\), and \(\boldsymbol{\alpha}^* \in \mathbb{R}^{500}\). Kernel regression works with the \(n \times n\) matrix, whatever the dimension of the feature space.

Q3Kernel \(k(x, x') = (1 + xx')^2\), training points \(x_1 = 1\), \(x_2 = -1\), and \(\boldsymbol{\alpha}^* = (2, -1)\). Predict at \(x_{\text{test}} = 0.5\).

\(k(1, 0.5) = 1.5^2 = 2.25\) and \(k(-1, 0.5) = 0.5^2 = 0.25\). So \(\hat y = 2(2.25) + (-1)(0.25) = 4.25\).

Q4Under the Gaussian noise model, why does the value of \(\sigma^2\) not change \(\hat{\mathbf{w}}_{\text{ML}}\)?

In the log-likelihood \(\sigma\) appears only as the positive factor \(1/(2\sigma^2)\) in front of the squared errors and inside the constant \(\log\frac{1}{\sqrt{2\pi}\sigma}\). Neither changes which \(\mathbf{w}\) minimizes \(\sum_i(\mathbf{w}^\top\mathbf{x}_i - y_i)^2\).

Q5You believe the label noise is Laplacian, not Gaussian. Which error should the regression minimize?

The sum of absolute errors \(\sum_i |\mathbf{w}^\top\mathbf{x}_i - y_i|\). The log of the Laplace density is \(-|\epsilon|/b\) plus a constant, so maximizing the likelihood minimizes the absolute errors.

Q6\(\sigma^2 = 2\) and \(XX^\top\) has eigenvalues 8, 2 and 0.1. What is the MSE of \(\hat{\mathbf{w}}_{\text{ML}}\)? What does \(\sigma^2\,\mathrm{tr}\big((XX^\top + I)^{-1}\big)\) come to?

\(2\,(1/8 + 1/2 + 1/0.1) = 2\,(0.125 + 0.5 + 10) = 21.25\). With \(\lambda = 1\): \(2\,(1/9 + 1/3 + 1/1.1) = 2\,(0.111 + 0.333 + 0.909) \approx 2.71\). Almost all of the drop comes from the eigenvalue 0.1. (This is the lecture's intuitive comparison, not the exact MSE of \(\hat{\mathbf{w}}_{\text{new}}\).)

Q7A symmetric matrix \(A\) has eigenvalues 3 and 0.5. What are the eigenvalues of \(A + 2I\) and of \((A + 2I)^{-1}\)?

Adding \(2I\) adds 2 to each eigenvalue: 5 and 2.5. The inverse has the reciprocals: 0.2 and 0.4.

Q810,000 points, 10-fold cross-validation, 20 candidate values of \(\lambda\). How many times is \(\hat{\mathbf{w}}_{\text{new}}\) fitted, and on how many points each time? What changes with leave-one-out?

\(10 \times 20 = 200\) fits, each on 9,000 points and validated on the remaining 1,000. Leave-one-out needs \(10{,}000 \times 20 = 200{,}000\) fits on 9,999 points each, which is why it is reserved for small data sets.