Ridge and LASSO: priors, penalties and sparse weights

At the end of Part 11 (L35), a new estimator \(\hat{\mathbf{w}}_{\text{new}} = (XX^\top + \lambda I)^{-1}X\mathbf{y}\) was proposed: for a suitable \(\lambda\) it has a smaller mean squared error than the maximum likelihood estimator, and \(\lambda\) is chosen by cross-validation. It looked like a lucky guess. L36 shows where it comes from, using the Bayesian approach of Part 8.

The Bayesian recipe has three pieces: a prior over the parameter, the usual likelihood of the data, and the posterior obtained by combining them. The parameter here is the weight vector \(\mathbf{w} \in \mathbb{R}^d\). It ranges over a continuous space, so the prior must be a density over all of \(\mathbb{R}^d\).

The likelihood is the model from Part 11: each label is the linear prediction plus Gaussian noise.

The lecture sets the noise variance to 1 only to keep the algebra less messy. Everything below goes through with a general \(\sigma^2\).

Which prior? A suitable one makes the posterior easy to compute. In Part 8, coin tosses had a Bernoulli likelihood, we used a beta prior, and the posterior was again a beta distribution. Because prior and posterior had the same functional form, updating the belief after seeing data was quick.

A prior with this property has a name. Given a likelihood, a prior is conjugate to it if prior × likelihood gives a posterior of the same form as the prior. The beta prior is conjugate to the Bernoulli likelihood.

Our likelihood is Gaussian, which is \(e\) raised to some power. Multiply it by a prior that is also \(e\) to some power and the exponents add, so the posterior is again \(e\) to some power. The natural choice is a Gaussian prior:

The mean \(\mathbf{0}\) is the zero vector in \(\mathbb{R}^d\), and the covariance has \(\gamma^2\) on the diagonal and zeros elsewhere. So each component \(w_j\) is Gaussian with mean 0 and variance \(\gamma^2\), and the components are independent. To draw a \(\mathbf{w}\) from this prior, sample each of \(w_1, \dots, w_d\) separately from \(\mathcal{N}(0, \gamma^2)\).

Given data \((\mathbf{x}_1, y_1), \dots, (\mathbf{x}_n, y_n)\), Bayes' rule gives the posterior up to a constant:

The points are generated independently, so the likelihood is a product over \(i\). With noise variance 1, the chance of seeing \(y_i\) at \(\mathbf{x}_i\) is proportional to \(\exp\!\big(-(y_i - \mathbf{w}^\top\mathbf{x}_i)^2/2\big)\); the factors \(1/\sqrt{2\pi}\) are dropped into the proportionality sign. The prior is a product over the \(d\) independent components, and the sum of squares in its exponent is the squared length of \(\mathbf{w}\):

The maximum a posteriori (MAP) estimate is the \(\mathbf{w}\) that maximizes the posterior. The missing constants do not move the maximizer, and log is increasing, so we can maximize the log instead. The log turns the product into a sum. Flipping the sign turns the max into a min:

Maximum likelihood had only the first term. The prior adds a second one: \(\lVert\mathbf{w}\rVert^2\) scaled by \(1/(2\gamma^2)\).

To minimize \(f\), take the gradient and set it to zero. The lecture states the gradient and asks you to verify it. Write the data as \(X \in \mathbb{R}^{d\times n}\) with the points \(\mathbf{x}_i\) as columns and \(\mathbf{y} \in \mathbb{R}^n\), so the first term is \(\tfrac12\lVert X^\top\mathbf{w} - \mathbf{y}\rVert^2\). Its gradient is \(X(X^\top\mathbf{w} - \mathbf{y})\), as in Part 10, and the gradient of \(\tfrac{1}{2\gamma^2}\mathbf{w}^\top\mathbf{w}\) is \(\mathbf{w}/\gamma^2\):

Standard result, not derived in the lecture: if the noise variance \(\sigma^2\) is kept, the same steps give \((XX^\top + \tfrac{\sigma^2}{\gamma^2} I)^{-1}X\mathbf{y}\).

The MAP estimate is exactly \(\hat{\mathbf{w}}_{\text{new}}\) from L35, with \(1/\gamma^2\) in place of \(\lambda\). There it was pulled out of a hat and justified by its mean squared error (Part 11). Here it comes out of a completely different line of thought: impose a prior on \(\mathbf{w}\) and ask for the most probable \(\mathbf{w}\) after seeing the data.

Only the name of the constant changes. In L35 we cross-validated \(\lambda\); here we would cross-validate \(\gamma^2\) in the same way. Either way it is a single number, and we pick the value that best fits held-out data.

MAP estimation for linear regression with the Gaussian prior \(\mathcal{N}(\mathbf{0}, \gamma^2 I)\) is the same as the estimator \((XX^\top + \lambda I)^{-1}X\mathbf{y}\) with \(\lambda = 1/\gamma^2\).

The route so far: solve linear regression as an optimization problem (Part 10); read it as maximum likelihood under Gaussian noise (Part 11); take the Bayesian version of that estimation problem; arrive at a new estimator that is better motivated and may well do better in practice.

L37 puts the two estimators side by side. Plain linear regression, which is also the maximum likelihood estimator, and the MAP estimator, written \(\hat{\mathbf{w}}_{\text{R}}\):

The R stands for ridge, a name with roots in classical statistics that the lecture does not go into. The problem on the right is called ridge regression. (Multiply the MAP objective by 2 to get exactly this form, with \(\lambda = 1/\gamma^2\).)

The first term is the original loss. The second is a regularizer: a penalty, added to the loss, for choosing some kinds of \(\mathbf{w}\) over others. Doing Bayesian modelling with a prior on \(\mathbf{w}\) is one way of saying that we prefer some \(\mathbf{w}\)'s; the regularizer is the same preference written as a penalty. The prior's variance turns into \(\lambda\). Such a number, which we fix by cross-validation instead of learning it from the loss, is a hyperparameter.

The prior has mean zero, so its mode, the single most probable \(\mathbf{w}\), is the zero vector. In the objective this shows up as the regularizer pulling the answer towards zero: it penalizes the length \(\lVert\mathbf{w}\rVert^2\), so it wants every component small. The loss term pulls the other way, towards a good fit. If several \(\mathbf{w}\)'s have the same loss, ridge prefers the one with the smallest length.

The regularizer cannot win outright, since \(\mathbf{w} = \mathbf{0}\) predicts nothing. Why push towards zero at all? Real data sets can have 10,000 or even a million features, and many of them may be redundant. When the label is a linear combination of the features, redundant features mean that many different combinations explain it equally well.

The height and weight example L37 · 5:39

Take three features: \(f_1\) = height, \(f_2\) = weight, and \(f_3 = 2\times\text{height} + 3\times\text{weight}\), as if some sensor recorded that combination. The label is a noisy version of \(3\times\text{height} + 4\times\text{weight}\). Several weight vectors reproduce it exactly:

The lecture leaves \(c_1, c_2\) unspecified; \(c_1 = -\tfrac12,\ c_2 = \tfrac32\) is the solution of \(2c_2 = 3,\ c_1 + 3c_2 = 4\).

The lecture's argument: we are probably better off not picking up redundant features, as \((3, 4, 0)\) does by ignoring \(f_3\). So push as many weights to zero as possible, assuming most features are redundant, while keeping the ones needed to keep the loss small. That is what the zero-mean prior is trying to encourage.

A worked check on these numbers (not in the lecture). Every exact fit has the form \((3 - 2t,\ 4 - 3t,\ t)\); \(t = 0, 1, \tfrac32\) give the three rows. Their squared lengths are 25, 3 and 2.5, and \(\lVert\mathbf{w}\rVert^2 = 25 - 36t + 14t^2\) is smallest at \(t = 9/7\): \(\mathbf{w} = (\tfrac37, \tfrac17, \tfrac97)\), with no zero entry. (This is also the least-norm solution that the pseudo-inverse of Part 10 would return.) So among equally good fits, the squared-length penalty spreads weight over all three features instead of dropping one. This gap between "towards zero" and "exactly zero" is what L38 makes precise and what LASSO (L39) addresses. For comparison, the smallest \(\lVert\mathbf{w}\rVert_1\) among these fits is at \(t = 4/3\): \((\tfrac13, 0, \tfrac43)\), which drops a feature.

\(\lambda\) plays the role of \(1/\gamma^2\). Suppose we strongly believe that most features are redundant, so most weights should be near zero. Then the prior should have a small variance \(\gamma^2\), with its mass packed around zero. Small \(\gamma^2\) means large \(\lambda\), and the minimization pays a large penalty for any \(\mathbf{w}\) with a large length. A weak belief means large \(\gamma^2\), small \(\lambda\), and the loss term dominates.

So the two terms of the objective balance how much loss we accept against how much redundancy we want to avoid.

L37 closes with questions. Why not directly try to make as many weights as possible exactly zero? Is ridge the only way to regularize? Would understanding ridge geometrically lead to better algorithms, or to a different formulation of linear regression? The answer to all of them is yes, and it leads to a modified version of linear regression (LASSO, in L39). The first step is the geometry of ridge.

L38 asks how the solutions of linear regression and ridge regression relate, and whether seeing that relation clearly suggests another way to regularize. Recall that \(\lVert\mathbf{w}\rVert^2 = \sum_{j} w_j^2\) is the squared L2 norm, the squared length of the vector \(\mathbf{w}\).

To draw pictures, take two features, say height and weight. Then \(\mathbf{w} = (w_1, w_2)\) and the parameter space is the whole plane, negative values included. Solve plain linear regression (equivalently maximum likelihood with Gaussian noise) and mark \(\hat{\mathbf{w}}_{\text{ML}}\) in the plane. Where is \(\hat{\mathbf{w}}_{\text{R}}\)?

The lecture uses an equivalence that it states without proof:

In (A), \(\mathbf{w}\) is unconstrained and pays a penalty. In (B), the loss is the same with no penalty, but \(\mathbf{w}\) may only come from a restricted set. In two dimensions that set is \(w_1^2 + w_2^2 \le \theta\): with equality it is a circle centred at the origin of radius \(\sqrt\theta\), so the set is the circle together with its inside. Ridge searches for the lowest-loss \(\mathbf{w}\) in this disc. \(\hat{\mathbf{w}}_{\text{ML}}\) was found with no restriction on length, so it may well lie outside it.

Standard context: the general tool behind such equivalences is constrained optimization with Lagrange multipliers, which the course develops for support vector machines in Part 16.

To find where in the disc ridge lands, look more closely at \(\hat{\mathbf{w}}_{\text{ML}}\). Call the objective \(f(\mathbf{w}) = \sum_i(\mathbf{w}^\top\mathbf{x}_i - y_i)^2\). Its value at \(\hat{\mathbf{w}}_{\text{ML}}\) need not be zero: the points may not have an exact linear relationship, so even the best \(\mathbf{w}\) suffers some loss. But by design \(f(\hat{\mathbf{w}}_{\text{ML}})\) is the smallest possible value; no \(\mathbf{w}\) has strictly less loss.

Now fix a constant \(c > 0\) and collect every \(\mathbf{w}\) whose loss is exactly \(c\) more than the minimum. If \(\hat{\mathbf{w}}_{\text{ML}}\) has loss 100, then \(S_{10}\) is the set of \(\mathbf{w}\)'s with loss 110.

In vector form, with \(X\) holding the data points as columns, \(X^\top\mathbf{w}\) is the vector of predictions \(\mathbf{w}^\top\mathbf{x}_i\). Every \(\mathbf{w} \in S_c\) satisfies

The lecture leaves the simplification as an exercise (expand both sides into four terms each; some cancel). The result is

where \(c'\) is a constant that depends on \(c\), the data and \(\hat{\mathbf{w}}_{\text{ML}}\), but not on \(\mathbf{w}\). We read this as an equation in \(\mathbf{w}\).

Completing the algebra

One clean way to do the exercise: split the residual vector in two and expand the squared length.

The middle term vanishes, because \(\hat{\mathbf{w}}_{\text{ML}}\) satisfies the normal equations \(XX^\top\hat{\mathbf{w}}_{\text{ML}} = X\mathbf{y}\) from Part 10. So

and the constant \(c'\) is in fact \(c\) itself.

What shape is this?

If \(XX^\top\) were the identity, the equation would read \(\lVert\mathbf{w} - \hat{\mathbf{w}}_{\text{ML}}\rVert^2 = c'\): a circle centred at \(\hat{\mathbf{w}}_{\text{ML}}\). In general \(XX^\top\) is some other matrix. It linearly transforms \(\mathbf{w} - \hat{\mathbf{w}}_{\text{ML}}\) before the dot product, which stretches the circle into an ellipse, still centred at \(\hat{\mathbf{w}}_{\text{ML}}\).

So the \(\mathbf{w}\)'s that make, say, 10 units more error than \(\hat{\mathbf{w}}_{\text{ML}}\) lie on an elliptical contour around it. Asking for 20 units more gives a larger concentric ellipse, and so on: the loss grows as we move outward from \(\hat{\mathbf{w}}_{\text{ML}}\) in any direction.

Put the two pictures together. \(\hat{\mathbf{w}}_{\text{ML}}\) sits outside the disc \(\lVert\mathbf{w}\rVert^2 \le \theta\), surrounded by its ellipses of extra loss. Grow the ellipse. Small ones miss the disc entirely; at some size the ellipse first touches it. That point of first contact is \(\hat{\mathbf{w}}_{\text{R}}\).

The reason: any \(\mathbf{w}\) other than \(\hat{\mathbf{w}}_{\text{ML}}\) costs extra loss, so ridge must pay something. It also has to stay inside the disc. Moving outward from \(\hat{\mathbf{w}}_{\text{ML}}\), the first point of the disc we reach is the feasible point with the least extra loss, which is exactly what problem (B) asks for.

What the picture says

• Ridge pushes the weights towards zero. It searches inside a disc, so \(\hat{\mathbf{w}}_{\text{R}}\) is a shorter vector than \(\hat{\mathbf{w}}_{\text{ML}}\).
• It does not necessarily make any weight exactly zero. The ellipse can touch the round boundary anywhere, and in a typical picture both coordinates of the contact point are nonzero.

Why exact zeros matter 25:40: if a feature's weight is exactly 0, the feature is not useful and can be thrown away. If the weight is merely small, we are stuck deciding whether to keep it or not.

The motivating question for L39 26:10: the regularizer decides the region in which we search for \(\mathbf{w}\). Can we change that region so that the growing ellipses are more likely to first touch it at a point where some coordinates are exactly zero?

The goal is a sparse solution: one with many entries exactly zero. The proposal is to regularize with the L1 norm instead of the squared L2 norm. The L1 norm of a vector just adds up the absolute values of its components:

The same kind of equivalence holds: for a given \(\lambda\) there is a \(\theta\) for which the two problems agree. It is generally a different \(\theta\) from the ridge case; the lecture reuses the letter on purpose. What matters is that L1 regularization is also a constrained problem, now with a constraint on the L1 norm of \(\mathbf{w}\).

In two dimensions \(|w_1| + |w_2| \le \theta\) is a diamond: a square turned by 45°, with its four corners on the axes at \((\pm\theta, 0)\) and \((0, \pm\theta)\). It replaces the disc as the region where we search for \(\mathbf{w}\).

Keep \(\hat{\mathbf{w}}_{\text{ML}}\) and its ellipses, and grow them as before. The circle bulges outward between the axes, so an ellipse tends to meet it on the bulge; that point was \(\hat{\mathbf{w}}_{\text{R}}\). The diamond is flat there. The ellipse grows past where it would have met the circle and is more likely to first touch at a corner. Call the contact point \(\hat{\mathbf{w}}_{\text{L}}\), the L1-regularized solution.

At a corner, one weight is nonzero and the other is exactly 0: the solution uses only one feature. In two dimensions that is the only kind of sparsity available, since setting both weights to zero is meaningless. The only points of the diamond with a zero coordinate are its four corners.

No. In two dimensions, a different position of \(\hat{\mathbf{w}}_{\text{ML}}\) or differently shaped ellipses can make the first contact happen on a flat edge, with both weights nonzero. (In the widget above, turn \(\hat{\mathbf{w}}_{\text{ML}}\) towards 45° and the diamond solution leaves the corner.) What people typically observe in high dimensions is that an L1 regularizer is more likely to give sparse solutions than an L2 regularizer: the L2 region has a bulge, while the L1 region is flat, so the contours move past the bulge and hit a flat point where many components can be zero.

This is an intuitive argument, not a proof. When L1 regularization is guaranteed to give a sparse solution is proved in advanced courses; this course does not do it.

Linear regression with an L1 regularizer is called LASSO, an acronym for Least Absolute Shrinkage and Selection Operator. Each word describes something we have already seen.

• Least: we minimize a loss.
• Absolute: the penalty is the L1 norm, a sum of absolute values.
• Shrinkage: we only search among \(\mathbf{w}\)'s of smaller length, so \(\mathbf{w}\) is shrunk. Ridge is a shrinkage method too; anything that restricts the length of the parameter is. LASSO shrinks using absolute values.
• Selection: the aim is more than smaller weights. Pushing many components to exactly 0 leaves nonzero weights only on the features that matter for the loss, so those features are selected and the rest left out.
• Operator: a fancy word for the regularizer.

LASSO is very popular for linear regression 10:16, especially when there are a lot of features and we expect most of them to be useless or redundant: LASSO pushes their weights to exactly 0.

If LASSO pushes weights to exactly zero, why would anyone still use ridge? The first point the lecture makes: ridge has a closed-form solution and LASSO does not.

For LASSO we cannot take the gradient, set it to zero and solve, because the penalty \(|w_j|\) is not differentiable at \(w_j = 0\). The problem has to be solved some other way.

Gradient methods follow the negative gradient direction (Part 10). When a function is not differentiable everywhere, an iterative method still works if a subgradient exists at every point. Where the gradient exists, it is the subgradient. Where it does not, a subgradient is the slope of a line through the point that stays entirely below the function: the function is lower-bounded by it everywhere.

Take a piecewise linear function. At a point where two pieces meet, there is no gradient, but many lines through that point stay below the whole function; the slope of each one is a subgradient, so there are many. At a point in the middle of a piece, only one line works, the piece itself extended. When there is only one subgradient, it is the gradient. In one dimension, think of a subgradient as a slope.

This matters because the L1 penalty in one dimension is the absolute value \(|x|\): two linear pieces joined at 0, not differentiable there. Any line through the origin with slope between −1 and 1 stays below \(|x|\). So at 0 every value in \([-1, 1]\) is a subgradient of \(|x|\); slope −1 and slope 1 are the two pieces themselves.

Why subgradients are useful: if the function to minimize is convex, subgradient descent can be used instead of gradient descent, and it also converges. The function does not need to be differentiable. Move along the negative of a subgradient instead of the negative gradient. Where several subgradients exist, it does not matter which one is picked; an arbitrary one will do.

1. Subgradient descent
2. Start from some \(\mathbf{w}^{0}\).
3. At step \(t\), pick any subgradient \(\mathbf{g}^{t}\) of the objective at \(\mathbf{w}^{t}\) (the gradient, wherever it exists).
4. Update \(\mathbf{w}^{t+1} = \mathbf{w}^{t} - \eta_t\,\mathbf{g}^{t}\) with a step size \(\eta_t\). The lecture does not go into how to choose it.
5. Repeat. For a convex objective the iterates approach the minimum.

The LASSO problem is convex: the squared loss is a convex quadratic, the L1 norm is a convex function, and a sum of convex functions is convex (convexity was introduced in Part 9). So LASSO is a convex optimization problem, and subgradient descent can find its minimum.

Standard, not stated in the lecture: one subgradient of the LASSO objective at \(\mathbf{w}\) is \(2X(X^\top\mathbf{w} - \mathbf{y}) + \lambda\mathbf{s}\), where \(s_j = \operatorname{sign}(w_j)\) when \(w_j \ne 0\) and \(s_j\) is any number in \([-1, 1]\) when \(w_j = 0\).

Subgradient descent is by no means the only way. LASSO has special structure, a squared loss plus an L1 penalty, and special-purpose methods exploit it. One example is IRLS, iteratively reweighted least squares. The idea is that we can solve least squares easily, in closed form, so IRLS uses least squares as a black box, repeatedly, to solve the LASSO problem. The course does not cover IRLS.

LASSO has no closed form. It gives sparse solutions, yet it needs an optimizer such as subgradient descent or a special method like IRLS. Ridge has a closed form; use it for small \(n\) and \(d\), and stochastic gradient descent when \(n\) or \(d\) is large.

The interactive above solves LASSO by coordinate descent, another special-purpose method; it is used there only to draw the picture.

L40 closes the regression unit with a high-level summary.

1. Least squares. Minimize the squared loss; it has a closed-form solution and a geometric interpretation (Part 10).
2. Computation. Concerns about the cost of the closed form led to gradient descent and stochastic gradient descent (Part 10).
3. Probabilistic view. With Gaussian noise, least squares is the maximum likelihood estimator (Part 11).
4. Bayesian counterpart. A zero-mean Gaussian prior on \(\mathbf{w}\) gives the L2-regularized version, ridge regression (L36, L37).
5. Geometry. The L2 regularizer pushes weights close to zero but not exactly to zero. The geometry of the problem suggested an L1 regularizer instead: LASSO, squared loss plus an L1 penalty (L38, L39).
6. Solving LASSO. Subgradients and subgradient descent, or special techniques such as IRLS (L40).

People use many other regularizers. The elastic net penalty mixes the LASSO and ridge penalties (in its usual form, \(\lambda_1\lVert\mathbf{w}\rVert_1 + \lambda_2\lVert\mathbf{w}\rVert^2\)). There is a whole cottage industry of regularizers for solutions with a particular structure. One example is group structure: features are collected from many sources, and you know some of them should be selected together or not at all. Such domain-specific constraints can be turned into regularizers, and the problem solved in regularized form.

Ridge and LASSO are the basic, most commonly used regularizers. If your domain knowledge says something about the structure of the solution you expect, consider turning it into a more meaningful regularizer of your own.

This completes regression. The course now turns to supervised learning from the classification point of view, starting in Part 13.

Conjugate prior

A prior that, multiplied by a given likelihood, gives a posterior of the same form: beta for Bernoulli, Gaussian for Gaussian.

Gaussian prior \(\mathcal{N}(\mathbf{0}, \gamma^2 I)\)

Each weight independent with mean 0 and variance \(\gamma^2\); prefers short \(\mathbf{w}\).

MAP estimate

The \(\mathbf{w}\) that maximizes the posterior. Here \((XX^\top + \tfrac{1}{\gamma^2}I)^{-1}X\mathbf{y}\).

Ridge regression

Squared loss plus \(\lambda\lVert\mathbf{w}\rVert^2\); solution \((XX^\top + \lambda I)^{-1}X\mathbf{y}\).

Regularizer

A penalty added to the loss that makes some kinds of \(\mathbf{w}\) preferable to others.

Hyperparameter \(\lambda\)

The strength of the penalty, playing the role of \(1/\gamma^2\); chosen by cross-validation.

Constrained form

Minimize the loss subject to \(\lVert\mathbf{w}\rVert^2 \le \theta\) (or \(\lVert\mathbf{w}\rVert_1 \le \theta\)); same solution as the penalty form for a matching \(\theta\).

Loss contour

The set of \(\mathbf{w}\) with loss \(f(\hat{\mathbf{w}}_{\text{ML}}) + c\): the ellipse \((\mathbf{w}-\hat{\mathbf{w}}_{\text{ML}})^\top XX^\top(\mathbf{w}-\hat{\mathbf{w}}_{\text{ML}}) = c\).

L1 norm

\(\lVert\mathbf{w}\rVert_1 = \sum_j|w_j|\). Its region in 2-D is a diamond with corners on the axes.

Sparse solution

A weight vector with many entries exactly 0, so those features can be dropped.

LASSO

Least Absolute Shrinkage and Selection Operator: squared loss plus \(\lambda\lVert\mathbf{w}\rVert_1\).

Subgradient

A \(\mathbf{g}\) with \(f(\mathbf{z}) \ge f(\mathbf{x}) + \mathbf{g}^\top(\mathbf{z}-\mathbf{x})\) for all \(\mathbf{z}\); the gradient where \(f\) is differentiable.

Subgradient descent

Step along the negative of any subgradient. Converges for convex functions, so it can solve LASSO.

Elastic net

A regularizer that mixes the LASSO and ridge penalties.