Linear regression: squared error, a closed form, and gradient descent

Parts 2 to 9 covered unsupervised learning: representation learning, clustering and estimation. L27 begins the second half of the course with a different paradigm, supervised learning, which the lecture calls the more popular one. From here on most algorithms are supervised, and the course also asks why some algorithms are better than others. The paradigms themselves were introduced in Part 1.

The data looks the same as before: points \(\mathbf{x}_1, \dots, \mathbf{x}_n \in \mathbb{R}^d\). Their coordinates are the features. Other names are attributes, and statisticians say predictors or predictor variables; the course says features. What is new is a label \(y_i\) attached to each point 1:37. The labels are the supervision, which is where the name comes from.

The values a label may take split supervised learning into different problems. Each kind needs its own type of algorithm, and the course will justify why as it goes.

• Binary classification. Two possible labels, hence "binary". The two encodings \(\pm1\) and \(\{0,1\}\) are almost equivalent for the course. For a new feature vector not seen in the data, predict which of the two buckets it falls in. For spam, each email is first turned into numbers somehow (how to do this well comes later); someone marks some emails as spam or not spam; the classifier must then label new emails automatically.
• Multi-class classification. A finite set of more than two classes. In digit recognition the image is the feature vector and the label is the digit written in it. People write the same digit in different ways, and the system should still recognize it.
• Regression. The label is any real number. To predict rainfall, the features might be temperature, precipitation and atmospheric pressure, and the label is an amount such as 100.23 cm. Regression is perhaps the earliest supervised problem and has strong roots in statistics.

The plan 6:39: regression for a couple of weeks, then binary classification for a long time, then a brief look at how binary-classification ideas give multi-class algorithms. Multi-class is not a major focus of the course.

The training data is \(\mathbf{x}_1, \dots, \mathbf{x}_n \in \mathbb{R}^d\) with labels \(y_1, \dots, y_n \in \mathbb{R}\). The goal is to learn a function \(h : \mathbb{R}^d \to \mathbb{R}\), from the feature space to the label space, that predicts well on new feature vectors.

There are endlessly many such functions. To choose among them we need a score for any given \(h\) on the data set we have: one number saying how good \(h\) is. Then we can search for the \(h\) with the best score.

Start with one point. Its label is \(y_i\) and \(h\) predicts \(h(\mathbf{x}_i)\). A natural comparison is the squared difference. Adding this over all points gives the error of \(h\) on the data set:

Why square? The error should be non-negative 3:17. With the plain difference \(h(\mathbf{x}_i) - y_i\), a function that outputs \(-\infty\) for every point would make the "error" as small as possible while being useless for prediction. Squared error is not the only reasonable choice; it is where the course starts, and it gets a proper justification later from a probabilistic view (Part 11).

The goal is now to find the \(h\) with the smallest error. How small can it be? It is a sum of squares, so it cannot go below 0, and it equals 0 exactly when

Such an \(h\) is easy to build: when given a training point \(\mathbf{x}_i\), output its label \(y_i\). Training error is then fully minimized. But this \(h\) has only memorized the mapping from \(\mathbf{x}_i\) to \(y_i\). It is the exam analogy from Part 1: memorizing the exercise problems scores well on the exercises and says little about the exam.

The lecture names the worry: we may be overfitting to the training data, so that doing well on it does not carry over to test data. Our goal is never the training data itself 8:15. A rainfall predictor that is perfect on the last three years and wrong from tomorrow onward is useless. What counts is performance on test data.

What went wrong? The squared error is not the problem. The problem is that the search ran over all functions from \(\mathbb{R}^d\) to \(\mathbb{R}\). A space that large always contains a function that memorizes the data, and that function fits the noise as well.

The remedy is to search a smaller space. That amounts to assuming that the true mapping from features to labels has some structure. The simplest structure is linear:

Each function in this set is fixed by its parameter vector \(\mathbf{w}\). In one dimension it is the line \(y = mx\). The line \(y = mx + c\) corresponds to \(\mathbf{w}^\top\mathbf{x} + b\); the lecture uses plain \(\mathbf{w}^\top\mathbf{x}\) for simplicity and says the shifted version is a straightforward adaptation, shown later 12:44.

The assumption being made 13:17: the real relation between features and labels is linear, and the labels may not look exactly linear only because of noise in the output. Searching all functions would fit that noise too, which is what overfitting means here. Searching only linear functions fits just the structure.

Keep the squared error and minimize it over \(h_{\mathbf{w}} \in \mathcal{H}_{\text{linear}}\) only. Since each \(h_{\mathbf{w}}\) is fixed by \(\mathbf{w}\), this is a search over \(\mathbf{w} \in \mathbb{R}^d\):

Can this objective reach 0? Yes: it is still a sum of non-negative numbers, and it is 0 when \(\mathbf{w}^\top\mathbf{x}_i = y_i\) for every \(i\). This time that is not overfitting. It means the labels follow the hypothesized linear structure exactly, with no noise, so the best structured fit happens to be perfect. Overfitting is what happens when the search space is bigger than the actual relation from inputs to outputs, so the fit starts to follow noise. Noise gets a probabilistic treatment later (Part 11); for now it is whatever the labels contain beyond the structure, and it is useless for prediction.

Matrices make the problem easier to handle. Stack the data points as the columns of \(X\), so \(X^\top\) has them as rows. Collect the labels in a vector \(\mathbf{y}\):

Then \(X^\top\mathbf{w} = [\mathbf{x}_1^\top\mathbf{w}, \dots, \mathbf{x}_n^\top\mathbf{w}]^\top\) is the vector of all \(n\) predictions. Subtracting \(\mathbf{y}\) leaves the \(n\) errors, and the squared L2 norm of that vector is the sum of their squares. A squared norm is also the vector transposed times itself:

This is an unconstrained quadratic optimization problem. Unconstrained, because any \(\mathbf{w} \in \mathbb{R}^d\) is allowed. Quadratic, because of the squared error. The simplest approach is to take the derivative with respect to \(\mathbf{w}\), which in higher dimensions is the gradient \(\nabla f(\mathbf{w})\), and set it to zero. The gradient is the vector of partial derivatives: for each component \(w_j\), differentiate while treating all other components as constants, and collect the results. Sometimes this gives a closed form, sometimes not.

The lecture leaves the gradient as an exercise. Expanding the objective:

With \(\nabla_{\mathbf{w}}(\mathbf{w}^\top A\mathbf{w}) = 2A\mathbf{w}\) for a symmetric \(A\) (here \(A = XX^\top\)), \(\nabla_{\mathbf{w}}(\mathbf{w}^\top\mathbf{b}) = \mathbf{b}\), and the last term constant:

Setting this to zero gives the equation the solution \(\mathbf{w}^*\) must satisfy. Standard name: the normal equations.

When is \(XX^\top\) not invertible 5:24? When it is not full rank, that is, it has a non-trivial null space. Then it maps many vectors to the same vector, a whole subspace landing on one point, so the normal equations have many solutions. The pseudo-inverse is a matrix with the good properties of an inverse that is available even when no inverse exists. Among the many solutions it returns the one with the smallest length, which is why it is called the least-norm solution 6:33. For now the lecture assumes \(XX^\top\) is invertible; nothing conceptual changes, and the pseudo-inverse is discussed later.

With a single feature (\(d = 1\)), \(XX^\top = \sum_i x_i^2\) is a number and \(X\mathbf{y} = \sum_i x_i y_i\), so \(w^* = \sum_i x_i y_i \big/ \sum_i x_i^2\). The gradient-descent widget in Lecture 31 uses this case.

What does \(\mathbf{w}^*\) depend on? \(X\) is \(d \times n\) and \(X^\top\) is \(n \times d\), so \(XX^\top\) is \(d \times d\). It is a scaled covariance matrix: \(XX^\top = \sum_i \mathbf{x}_i\mathbf{x}_i^\top = n\,C\), where \(C = \frac1n\sum_i\mathbf{x}_i\mathbf{x}_i^\top\) is the matrix of Part 2. So, as in PCA (Part 3), a covariance-like matrix shows up, here through its inverse. Its exact role is explained later.

The key difference is the labels. PCA was completely unsupervised. Regression tries to match the labels as well as possible, so the labels must shape the answer, and they enter through \(X\mathbf{y}\). Given a data set, \(\mathbf{w}^*\) can be computed directly in code from the data and the labels, though possibly at a high computational cost (Lecture 31).

Having a closed form is a great start, but \((XX^\top)^{\dagger}X\mathbf{y}\) still needs an interpretation. L29 lists four questions. Answering them gives a solid understanding of linear regression:

1. Geometry. Can the solution be described geometrically? Answered in Lecture 30 below.
2. Computation. Are there computational issues, as there were for PCA (Part 4)? Answered in Lecture 31 below.
3. Non-linear relationships. The features were assumed to be linearly related to the labels. What if the relation is quadratic or cubic? PCA faced a similar question (Part 5), but the answer here has a different flavour. This becomes kernel regression in Part 11.
4. Noise. What exactly is the noise mentioned in L28? The right language for noise is probability, so: is there a probabilistic view of linear regression? Also in Part 11.

The answers lead to two algorithmic changes to the standard problem, covered in Part 12.

To be able to draw it, take \(d = 2\) features and \(n = 3\) data points: the height and weight of three people. Everything said here carries over to higher dimensions.

The usual picture would plot the three data points, each in \(\mathbb{R}^2\). This picture is different. It draws one vector per feature, in \(\mathbb{R}^n = \mathbb{R}^3\):

• \(\mathbf{f}_1 \in \mathbb{R}^3\): the heights of all three people (in general an \(n\)-dimensional vector);
• \(\mathbf{f}_2 \in \mathbb{R}^3\): their weights;
• \(\mathbf{y} \in \mathbb{R}^3\): their labels, one per person 3:22.

These feature vectors are the rows of \(X\), so \(X^\top = [\,\mathbf{f}_1 \ \ \mathbf{f}_2\,]\). Combining them gives other vectors, such as \(3\mathbf{f}_1 + 2.5\mathbf{f}_2\) or \(-5\mathbf{f}_1 + 3\mathbf{f}_2\) (the lecture's arbitrary examples). All combinations \(a\,\mathbf{f}_1 + b\,\mathbf{f}_2\) together form a 2-D plane sitting inside \(\mathbb{R}^3\), the plane spanned by \(\mathbf{f}_1\) and \(\mathbf{f}_2\) 4:30.

In matrix form a combination is \([\,\mathbf{f}_1 \ \ \mathbf{f}_2\,]\begin{bmatrix} a \\ b \end{bmatrix} = X^\top\boldsymbol\alpha\). So the prediction vector \(X^\top\mathbf{w}\) of any linear model is a point of this plane, and every point of the plane is the prediction vector of some \(\mathbf{w}\).

The label vector \(\mathbf{y}\) need not lie in the plane. We want to explain the labels as a combination of height and weight, so we look for the point of the plane closest to \(\mathbf{y}\). That point is the projection of \(\mathbf{y}\) onto the subspace spanned by the feature vectors 5:33.

Being in the plane, the closest point is \(\alpha_1^*\mathbf{f}_1 + \alpha_2^*\mathbf{f}_2 = X^\top\boldsymbol\alpha^*\) for some unknown \(\boldsymbol\alpha^*\). Being the closest, the gap \(\mathbf{y} - X^\top\boldsymbol\alpha^*\), the vector to add to reach \(\mathbf{y}\), meets the plane at a right angle. In particular it is orthogonal to \(X^\top\boldsymbol\alpha^*\) itself:

Does the closed-form w* satisfy it? L30 · 8:48

Substitute \(\boldsymbol\alpha^* = \mathbf{w}^* = (XX^\top)^{-1}X\mathbf{y}\) from L29 and check that the left side is 0. Write \(A = XX^\top\) and use the inverse for readability. The lecture leaves the algebra as an exercise and points to the key fact: \(A\) is symmetric, so \(A^{-1}\) is symmetric too, and \(\mathbf{w}^{*\top} = \mathbf{y}^\top X^\top A^{-1}\). The two terms become

They are equal, so their difference is 0. With the pseudo-inverse the same steps work, because \(A^{\dagger}AA^{\dagger} = A^{\dagger}\) is one of its defining properties.

Standard fact, not in the lecture: the gap \(\mathbf{y} - X^\top\mathbf{w}^*\) is perpendicular to every vector in the plane, in particular to each feature vector \(\mathbf{f}_j\). Stacking those \(d\) conditions gives \(X(\mathbf{y} - X^\top\mathbf{w}^*) = \mathbf{0}\), which is exactly the normal equations \(XX^\top\mathbf{w}^* = X\mathbf{y}\) of L29. Projection onto a single line in Part 2 is the one-feature version of this picture.

An exercise from the lecture 12:07: what if \(\mathbf{y}\) already lies in the plane? Then its projection is \(\mathbf{y}\) itself. The lecture leaves the meaning as an exercise; in terms of the error it means \(X^\top\mathbf{w}^* = \mathbf{y}\), so the squared error is 0 and every label is exactly \(\mathbf{w}^{*\top}\mathbf{x}_i\). That is the noise-free case from L28.

L31 takes up the second question, computation. The closed form \(\mathbf{w}^* = (XX^\top)^{\dagger}X\mathbf{y}\) is convenient: given data, compute \(\mathbf{w}^*\) once and use it to predict on test data. The question is how expensive that computation is.

The obvious culprit is the inverse (or pseudo-inverse) of the \(d \times d\) matrix \(XX^\top\). Inverting it takes on the order of \(d^3\) operations. That is the same cost as the eigenvector computation for PCA (Part 4), and the lecture notes the two are equally hard: once the eigenvectors are known, the inverse follows by inverting the eigenvalues, \(XX^\top = \sum_j \lambda_j\mathbf{u}_j\mathbf{u}_j^\top \Rightarrow (XX^\top)^{-1} = \sum_j \frac{1}{\lambda_j}\mathbf{u}_j\mathbf{u}_j^\top\) 2:46. Special-purpose algorithms can do slightly better, but \(d^3\) is the general picture.

With a billion features, inverting a billion-by-billion matrix takes far too much computation. Can we reach the same \(\mathbf{w}^*\) without computing an inverse? The theory is the foundation, but what goes wrong in practice matters too.

\(\mathbf{w}^*\) is the minimizer of an unconstrained problem: the sum (or average) of squared errors, over all of \(\mathbb{R}^d\). So general-purpose tools for unconstrained minimization apply, and the most popular of them is gradient descent.

Gradient descent is an iterative way to find minimizers that uses only first-order information, the gradient. Start at some guess \(\mathbf{w}_0\), an arbitrary point in \(\mathbb{R}^d\). The gradient at the current point gives the direction in which the function increases fastest. We want to decrease it, so we move the opposite way: to climb a mountain follow the gradient, to reach the bottom follow its negative.

At the new point the gradient may point a different way, so the step is repeated. If the step size is chosen carefully, repeating it eventually reaches a local minimum, a point where the gradient is 0. There the update stops changing \(\mathbf{w}\).

Local and global minima L31 · 7:48

A local minimum need not be the global minimum. On a wiggly function with several valleys, gradient descent started on one side slides into the nearest valley and stops there, even if a deeper valley exists elsewhere.

On a "nice" function without wiggles, every local minimum is a global minimum, and gradient descent converges to it wherever it starts 8:53. The squared-error objective is of this kind. It is quadratic in \(\mathbf{w}\), and the lecture places it among the convex functions (see Part 9), for which all local minima are global. So gradient descent converges to the optimal \(\mathbf{w}^*\). (Standard detail: it is convex because its matrix of second derivatives, \(2XX^\top\), is positive semidefinite.)

Does gradient descent actually remove the inverse here? The objective and its gradient, as in L29:

The lecture again suggests verifying the gradient component by component. Each partial derivative is \(\partial f/\partial w_j = \sum_i 2(\mathbf{w}^\top\mathbf{x}_i - y_i)\,x_{ij}\), and stacking them over \(j\) gives \(2X(X^\top\mathbf{w} - \mathbf{y})\), the same vector.

There is no inverse anywhere 12:09. Each step is the current \(\mathbf{w}\) minus a scalar times a quantity computed directly. Gradient descent mimics \((XX^\top)^{-1}X\mathbf{y}\) in many steps instead of one, and converges to the same \(\mathbf{w}^*\). So it still runs when inverting is the bottleneck, even with millions of features.

Why "chosen carefully" matters (standard analysis, not derived in the lecture). With one feature and the average error, \(f(w) = a w^2 - 2bw + c\) with \(a = \frac1n\sum_i x_i^2\), \(b = \frac1n\sum_i x_i y_i\) and \(w^* = b/a\). One step gives \(w_{t+1} - w^* = (1 - 2\eta a)(w_t - w^*)\). The gap shrinks only if \(|1 - 2\eta a| < 1\), that is \(0 < \eta < 1/a\). Between \(1/(2a)\) and \(1/a\) the factor is negative, so every step lands on the other side of \(w^*\).

Gradient descent fixed the problem of large \(d\). What if \(n\), the number of data points, is large, say hundreds of millions? Gradient descent is iterative and needs the gradient in every round, and the gradient needs \(XX^\top\) 13:47. That matrix is only \(d \times d\), but it comes from multiplying a \(d \times n\) matrix by an \(n \times d\) one, so a dependence on \(n\) hides inside it. With billions of points even computing \(XX^\top\) may be too hard. (Writing the gradient as \(2X(X^\top\mathbf{w} - \mathbf{y})\) avoids forming \(XX^\top\), but each step still touches all \(n\) points.)

Can gradient descent be adapted to this situation? Surprisingly, yes. A variant called stochastic gradient descent (SGD) removes the \(n\) problem in the way gradient descent removed the \(d\) problem, and it is used in practice all the time. There are several variants; this is one simple version.

1. Stochastic gradient descent for linear regression
2. Start from some \(\mathbf{w}_0 \in \mathbb{R}^d\). Fix a budget of \(T\) rounds (or use a convergence criterion, as for gradient descent).
3. In each round \(t = 0, 1, \dots, T-1\), sample \(k\) data points uniformly at random from all \(n\). Put them in \(\tilde X \in \mathbb{R}^{d \times k}\) and their labels in \(\tilde{\mathbf{y}} \in \mathbb{R}^k\).
4. Pretend the sample is the entire data set and take one gradient step for it: \(\mathbf{w}_{t+1} = \mathbf{w}_t - \eta_t \cdot 2\big(\tilde X\tilde X^\top\mathbf{w}_t - \tilde X\tilde{\mathbf{y}}\big)\).
5. After \(T\) rounds output the average iterate \(\bar{\mathbf{w}}_T = \frac1T\sum_{t=1}^{T}\mathbf{w}_t\), not the last one.

Sampling is cheap 16:34. With a billion points and \(k = 100\), draw 100 indices uniformly from 1 to \(10^9\), which is fast, then fetch those 100 points, which is easy if the data is stored sensibly. The matrix \(\tilde X\) is only \(d \times k\), so the step is manageable 19:11. The most important detail is that every round draws a fresh sample; the \(k\) points are not the same from round to round.

Noisy steps and the averaged output L31 · 20:12

The gradient computed from a sample is a noisy version of the full gradient, so SGD moves in noisy directions. Because the sampling is uniform, the noise cancels out over time and SGD still makes progress toward the same \(\mathbf{w}^*\) that gradient descent approaches. (Standard reasoning: averaged over the random choice of sample, the sampled gradient points the same way as the full gradient, scaled by \(k/n\).)

One caveat 20:46. Gradient descent outputs its last iterate \(\mathbf{w}_T\). In SGD the last round may have drawn an unlucky sample and moved in a direction that leads away from the minimizer. So SGD outputs the average of all the iterates it has seen.

\(k\) is the batch size 22:56. Pick it to match the computation you can afford, larger or smaller. The guarantee holds whatever the batch size; a smaller batch typically needs somewhat more rounds than standard gradient descent to reach the minimizer. That extra work is still far more manageable than computing \(XX^\top\) when \(n\) is very large. In short, switching to gradient descent removed the problem with large \(d\), and switching to its stochastic version removed the problem with large \(n\). In the widget above, tick Stochastic to watch \(w_t\) jitter around \(w^*\) while the running average \(\bar w\) settles near it.

L31 ends with the next question 24:01: why linear? If the relation between features and labels is quadratic, or otherwise non-linear, what then? Just as PCA was kernelized into kernel PCA (Part 5), linear regression gets a kernelized version, kernel regression, in Part 11.

Supervised learning

Learning from data points \(\mathbf{x}_i\) that each come with a label \(y_i\), the supervision.

Features and labels

Features (attributes, predictors) are the coordinates of \(\mathbf{x}_i \in \mathbb{R}^d\); the label \(y_i\) is what we want to predict.

Binary / multi-class classification

Labels from two values (\(\pm1\) or \(\{0,1\}\)) / from a finite set such as \(\{0, \dots, 9\}\).

Regression

Supervised learning with real-valued labels, such as a rainfall amount.

Squared error

\(\sum_i (h(\mathbf{x}_i) - y_i)^2\): a non-negative score of a function \(h\) on the data set.

Overfitting

Fitting the training data, noise included, so closely (for example by memorizing) that performance on test data suffers.

\(\mathcal{H}_{\text{linear}}\)

The restricted search space of linear functions \(h_{\mathbf{w}}(\mathbf{x}) = \mathbf{w}^\top\mathbf{x}\), \(\mathbf{w} \in \mathbb{R}^d\).

Linear regression

\(\min_{\mathbf{w}} \sum_i(\mathbf{w}^\top\mathbf{x}_i - y_i)^2 = \min_{\mathbf{w}}\lVert X^\top\mathbf{w} - \mathbf{y}\rVert^2\), with \(X \in \mathbb{R}^{d\times n}\).

Normal equations

\(XX^\top\mathbf{w}^* = X\mathbf{y}\), from setting the gradient to zero; solved by \(\mathbf{w}^* = (XX^\top)^{\dagger}X\mathbf{y}\).

Pseudo-inverse \(A^{\dagger}\)

Stands in for the inverse when \(A\) is not invertible; it returns the least-norm solution.

Projection view

\(X^\top\mathbf{w}^*\) is the projection of the label vector onto the subspace spanned by the feature vectors.

Gradient descent, step size

\(\mathbf{w}_{t+1} = \mathbf{w}_t - \eta_t\nabla f(\mathbf{w}_t)\); \(\eta_t\) sets how far each step goes.

Convex function

A function whose local minima are all global; the squared-error objective is one.

SGD, batch size

Gradient steps on \(k\) randomly sampled points per round (\(k\) is the batch size), with the average iterate as output.