Ridge (L2) and Lasso (L1) Regression

$$

Goals

• Introduce ridge (L2) and Lasso (l1) regression
  • As a complexity penalty
  • As a tuneable hierarchy of models to be selected by cross-validation
  • As a Bayesian posterior estimate

Reading

These notes are a supplement to Section 6.2 of Gareth et al. (2021).

Approximately colinear regressors

Sometime it makes sense to penalize very large values of $\hat{\beta }$.

Consider the following contrived example. Let

$$\begin{array}{rl}{y}_n=& z_n+{\epsilon }_n\text{where}{z}_n\sim \mathcal{N}(0,1)\text{and}\epsilon \sim \mathcal{N}(0,1).\end{array}$$

If we regress $y\sim \beta z$, then $\hat{\beta }\sim \mathcal{N}(1,1/N)$ approximately, with no problems. But suppose we actually also observe $x_n=z_n+{ϵ}_n$ with ${ϵ}_n\sim \mathcal{N}(0,\delta )$ for very small $\delta \ll 1$. Then $x_n$ and $z_n$ are highly correlated:

$${\mathit{M}}_{\mathit{X}}={\mathbb{E}}_{}\left[\left(\begin{array}{c}{x}_n\\ z_n\end{array}\right)\left(\begin{array}{cc}{x}_n& z_n\end{array}\right)\right]=\left(\begin{array}{cc}1+\delta & 1\\ 1& 1\end{array}\right)\Rightarrow {\mathit{M}}_{\mathit{X}}^{-1}=\frac{1}{\delta }\left(\begin{array}{cc}1& -1\\ -1& 1+\delta \end{array}\right).$$

Therefore, if we try to regress $y\sim {\beta }_{x}x+{\beta }_{z}z$, we get

$${\mathrm{Cov}}_{}(\left(\begin{array}{c}{\hat{\beta }}_x\\ {\hat{\beta }}_z\end{array}\right)|\mathit{X},\mathit{Z})=\frac{1}{N}{\mathit{M}}_{\mathit{X}}^{-1}=\frac{1}{N}\frac{1}{\delta }\left(\begin{array}{cc}1& -1\\ -1& 1+\delta \end{array}\right).$$

Note that ${\mathrm{Var}}_{}\left({\beta }_x\right)=\frac{N}{\delta }$, which is very large. With high probability, we will estimate $\hat{\beta }$ that has a very large magnitude, ${\Vert \hat{\beta }\Vert }_2^2$, although its two components should nearly the negative of one another. This could be problematic in practice. For example, in our test set or application, if $x_n$ and $z_n$ are not as well–correlated as in our training set, this will lead to crazy and highly variable predicted values.

Does it make sense to permit such a large variance? Wouldn’t it be better to choose slightly smaller $\hat{\beta }$, which are in turn somewhat more “stable”?

Penalizing large regressors with ridge regression

Recall that one perspective on regression is that we choose $\hat{\beta }$ to minimize the loss

$$\hat{\beta }:=\underset{\beta }{\mathrm{argmin}}\sum _{n=1}^N(y_n-{\beta }^{⊺}{x}_n{)}^2=:RSS(\beta ).$$

We motivated this as an approximation to the expected predicted loss, $L(\beta )={\mathbb{E}}_{}[y_{\mathrm{new}},x_{\mathrm{new}}](y_{\mathrm{new}}-{\beta }^{⊺}{x}_{\mathrm{new}}{)}^2$. But that made sense when we had a fixed set of regressors, and have shown that the correspondence breaks down when we are searching the space of regressors. In particular, $RSS(\beta )$ always decreases as we add more regressors, but $L(\beta )$ may not.

Instead, let us consider minimizing $RSS(\beta )$, but with an additional penalty for large $\hat{\beta }$. There are many ways to do this! But one convenient one is as follows. For now, pick a $\lambda$, and choose $\hat{\beta }$ to minimize:

$$\hat{\beta }(\lambda ):=\underset{\beta }{\mathrm{argmin}}{L}_{ridge}(\beta ,\lambda ):=RSS(\beta )+\lambda {\Vert \beta \Vert }_2^2=RSS(\beta )+\lambda \sum _{p=1}^P{\beta }_p^2.$$

This is known as ridge regression, L2–penalized regression. The latter is because the penalty ${\Vert \beta \Vert }_2^2$ is the L2 norm of the regressor; next time we will study the L1 version, which is also known as the Lasso.

The term $\lambda {\Vert \beta \Vert }_2^2$ is known as a “regularizer,” since it imposes some “regularity” to the estimate $\hat{\beta }(\lambda )$. Note that

• As $\lambda \to \mathrm{\infty }$, then $\hat{\beta }(\lambda )\to \mathbf{0}$
• When $\lambda =0$, then $\hat{\beta }(\lambda )=\hat{\beta }$, the OLS estimator.

So the inclusion of $\lambda$ is to “shrink” the estimate $\hat{\beta }(\lambda )$ towards zero. Note that since the ridge loss has an extra penalty for the norm, it is impossible for the OLS solution to have a smaller norm than the ridge solution.

The ridge regression regularizer has the considerable advantage that the optimum is available in closed form, since

$$\begin{array}{rl}{L}_{ridge}(\beta ,\lambda )=& (\mathit{Y}-\mathit{X}\beta {)}^{⊺}(\mathit{Y}-\mathit{X}\beta )+\lambda {\beta }^{⊺}\beta \\ =& {\mathit{Y}}^{⊺}\mathit{Y}-2{\mathit{Y}}^{⊺}\mathit{X}\beta +{\beta }^{⊺}{\mathit{X}}^{⊺}\mathit{X}\beta +\lambda {\beta }^{⊺}\beta \\ =& {\mathit{Y}}^{⊺}\mathit{Y}-2{\mathit{Y}}^{⊺}\mathit{X}\beta +{\beta }^{⊺}({\mathit{X}}^{⊺}\mathit{X}+\lambda {\mathit{I}}_{})\beta \Rightarrow \\ \frac{\partial L_{ridge}(\beta ,\lambda )}{\partial \beta }=& -2{\mathit{X}}^{⊺}\mathit{Y}+2({\mathit{X}}^{⊺}\mathit{X}+\lambda {\mathit{I}}_{})\beta \Rightarrow \\ \hat{\beta }(\lambda )=& {({\mathit{X}}^{⊺}\mathit{X}+\lambda {\mathit{I}}_{})}^{-1}{\mathit{X}}^{⊺}\mathit{Y}.\end{array}$$

Note that ${\mathit{X}}^{⊺}\mathit{X}+\lambda {\mathit{I}}_{}$ is always invertible if $\lambda >0$, even if $\mathit{X}$ is not full–column rank. In this sense, the ridge regression can deal safely with colinearity.

Exercise

Prove that ${\mathit{X}}^{⊺}\mathit{X}+\lambda {\mathit{I}}_{}$ is invertible if $\lambda >0$. Hint: using the fact that ${\mathit{X}}^{⊺}\mathit{X}$ positive semi–definite because it’s symmetric, find a lower bound on the smallest eigenvalue of ${\mathit{X}}^{⊺}\mathit{X}+\lambda {\mathit{I}}_{}$.

Standardizing regressors

Suppose we re-scale one of the regressors $x_{np}$ by some value $\alpha$ for a very small $\alpha \ll 1$, i.e., regressing on $x_{np}^{\prime }=\alpha x_{np}$ instead of $x_{np}$. As we know, the OLS minimizer ${\hat{\beta }}_p^{\prime }={\hat{\beta }}_p/\alpha$ and the fitted value $\hat{\mathit{Y}}$ is unchanged at $\lambda =0$. But for a particular $\lambda >0$, how does this affect the ridge solution? We can write

$$\lambda {\hat{\beta }}_p^{\prime 2}=\frac{\lambda }{{\alpha }^2}{\hat{\beta }}_p^2.$$

That is, we will effectively “punish” large values of ${\hat{\beta }}_p^{\prime }$ much more than we would “punish” the corresponding values of $\hat{\beta }$. In turn, for a particular $\lambda$, we will tend to set ${\hat{\beta }}_p^{\prime }<{\hat{\beta }}_p$ (although this is not necessarily the case).

The point is that re-scaling the regressors affects the meaning of $\lambda$. Correspondingly, if different regressors have very different typical scales, such as age versus income, then ridge regression will drive their coefficients to zero to very different degrees.

Similarly, it often doesn’t make sense to penalize the constant, so we might take ${\beta }_1$ to be the constant ($x_{n1}=1$), and write

$$\hat{\beta }(\lambda ):=RSS(\beta )+\lambda {\Vert \beta \Vert }_2^2=RSS(\beta )+\lambda \sum _{p=2}^P{\beta }_p^2.$$

But this gets tedious, and assumes we have included a constant.
Instead, we might invoke the FWL theorem, center the response and regressors at their mean values, and then do penalized regression.

For both these reasons, before performing ridge regression (or any other similar penalized regression), we typically standardize the regressors, defining

$$x_n^{\prime }:=\frac{x_n-{\overline{x}}_n}{\sqrt{\frac{1}{N}\sum _{n=1}^N(x_n-\overline{x}{)}^2}}.$$

We then run ridge regression on ${\mathit{x}}^{\prime }$ rather than $x$, so that we

• Don’t penalize the constant term and
• Penalize every coefficient the same regardless of its regressor’s typical scale.

A complexity penalty

Suppose that $y_n={\beta }^{⊺}{\mathit{x}}_n+{\epsilon }_n$ for some $\beta$. Note that, for a fixed ${\mathit{x}}_{\mathrm{new}}$, $y_{\mathrm{new}}$, and fixed $\mathit{X}$,

$$\begin{array}{rl}{\mathbb{E}}_{}[y_{\mathrm{new}}-{\mathit{x}}_{\mathrm{new}}^{⊺}\hat{\beta }(\lambda )]=& {\mathit{x}}_{\mathrm{new}}^{⊺}(\beta -{\mathbb{E}}_{}[\hat{\beta }(\lambda )])\\ =& {\mathit{x}}_{\mathrm{new}}^{⊺}({\mathit{I}}_{}-{({\mathit{X}}^{⊺}\mathit{X}+\lambda {\mathit{I}}_{})}^{-1}({\mathit{X}}^{⊺}\mathit{X}))\beta \end{array}$$

That is, as $\lambda$ grows, $\hat{\beta }(\lambda )$ becomes biased. However, by the same reasoning as in the standard case, under homoskedasticity,

$${\mathrm{Cov}}_{}(\hat{\beta }(\lambda ))={\sigma }^2{({\mathit{X}}^{⊺}\mathit{X}+\lambda {\mathit{I}}_{})}^{-1}({\mathit{X}}^{⊺}\mathit{X}){({\mathit{X}}^{⊺}\mathit{X}+\lambda {\mathit{I}}_{})}^{-1},$$

which is smaller that ${\mathrm{Cov}}_{}\left(\hat{\beta }\right)$ in the sense that ${\mathrm{Cov}}_{}\left(\hat{\beta }\right)-{\mathrm{Cov}}_{}(\hat{\beta }(\lambda ))$ is a positive definite matrix for $\lambda >0$ and full-rank $\mathit{X}$. In this sense, the family $\hat{\beta }(\lambda )$ is a one-dimensional family that trades off bias and variance. We can thus use cross validation to choose $\lambda$ that minimizes estimated MSE.

Constrained optimization and the minimum norm interpolator

We can rewrite the ridge loss in a suggestive way. Fix $\lambda =\hat{\lambda }$, write $b={\Vert \hat{\mathit{\beta }}(\hat{\lambda })\Vert }_2^2$, and write

$${\mathcal{L}}_{ridge}^{\prime }(\mathit{\beta },\lambda )=RSS(\mathit{\beta })+\lambda ({\Vert \mathit{\beta }\Vert }_2^2-b).$$

Since $b$ is fixed, $\hat{\mathit{\beta }}(\hat{\lambda })$ still satisfies

$${\frac{\partial {\mathcal{L}}_{ridge}(\mathit{\beta },\lambda )}{\partial \mathit{\beta }}|}_{\hat{\mathit{\beta }},\hat{\lambda }}={\frac{\partial RSS(\mathit{\beta })}{\partial \mathit{\beta }}|}_{\hat{\mathit{\beta }}}+\hat{\lambda }{\frac{\partial {\Vert \mathit{\beta }\Vert }_2^2}{\partial \mathit{\beta }}|}_{\hat{\mathit{\beta }}}=\mathbf{0}.$$

So the optimum is actually the same. The loss ${\mathcal{L}}_{ridge}^{\prime }(\mathit{\beta },\lambda )$ is the Lagrange multiplier version of the constrained optimization problem

$$\hat{\mathit{\beta }}(b):=\underset{\mathit{\beta }:{\Vert \mathit{\beta }\Vert }_2^2\le b}{\mathrm{argmin}}RSS(\mathit{\beta }).$$

Taking $b={\Vert \hat{\mathit{\beta }}(\lambda )\Vert }_2^2$, we see that, for every $\lambda$, there is a $b$, and vice–versa. The ridge regression is thus equivalent to minimizing the squared error loss subject to the constraint that it lies within an L2 ball. Making the ball larger allows better fit to the data, but by using larger $\mathit{\beta }$. This intuition is particularly useful when contrasting ridge regression with the lasso.

This intuition is also useful when understanding what happens as $\lambda \to 0$. Write $r=ESS(\hat{\mathit{\beta }}(\hat{\lambda }))$, and note that we could also have written

$${\mathcal{L}}_{ridge}^{″}(\mathit{\beta },\lambda )=\frac{1}{\lambda }(RSS(\mathit{\beta })-r)+{\Vert \mathit{\beta }\Vert }_2^2.$$

As before, for fixed $\lambda$ (and so fixed $r$), ${\mathcal{L}}_{ridge}^{″}(\mathit{\beta },\lambda )$ has the same ridge optimum. However, we can write

$$\hat{\mathit{\beta }}(r):=\underset{\mathit{\beta }:RSS(\mathit{\beta })<r}{\mathrm{argmin}}{\Vert \mathit{\beta }\Vert }_2^2.$$

We can thus equivalently interpret the ridge estimator as the one that produces the smallest $\mathit{\beta }$ in the L2 norm, subject to the $RSS$ being no larger than $r$. As $\lambda \to 0$, infinite weight gets put on the $RSS(\mathit{\beta })$ term, and $r\to 0$ if $\mathit{X}$ is rank $N$ or higher. So we see that

$$\underset{r\to 0}{lim}\hat{\mathit{\beta }}(r)=\underset{\mathit{\beta }:RSS(\mathit{\beta })=0}{\mathrm{argmin}}{\Vert \mathit{\beta }\Vert }_2^2.$$

That is, when there are many $\mathit{\beta }$ that have $RSS(\mathit{\beta })=0$ — i.e., which “interpolate” the data — the limiting ridge estimator chooses the one with minimum norm. This is call the “ridgeless interpolator.” (See the “ridgeless” lecture in Tibshirani (2023).)

A Bayesian posterior

Another way to interpret the ridge penalty is as a Bayesian posterior mean. If

$$\mathit{\beta }\sim \mathcal{N}(\mathbf{0},{\sigma }_{\beta }^2{\mathit{I}}_{})\text{and}\mathit{Y}|\mathit{\beta },\mathit{X}\sim \mathcal{N}(\mathit{X}\mathit{\beta },{\sigma }^2{\mathit{I}}_{}),$$

then

$$\mathit{\beta }|\mathit{Y},\mathit{X}\sim \mathcal{N}({({\mathit{X}}^{⊺}\mathit{X}+\frac{{\sigma }^2}{{\sigma }_{\beta }^2}{\mathit{I}}_{})}^{-1}{\mathit{X}}^{⊺}\mathit{Y},{\sigma }^2{({\mathit{X}}^{⊺}\mathit{X}+\frac{{\sigma }^2}{{\sigma }_{\beta }^2}{\mathit{I}}_{})}^{-1}).$$

One way to derive this is to recognize that, if $\mathit{\beta }\sim \mathcal{N}(\mu ,\mathrm{\Sigma })$, then

$$\log {\mathbb{P}}_{}(\mathit{\beta }|\mu ,\mathrm{\Sigma })=-\frac{1}{2}{\beta }^{⊺}{\mathrm{\Sigma }}^{-1}\beta +{\beta }^{⊺}{\mathrm{\Sigma }}^{-1}\mu +\text{Terms that do not depend on }\beta .$$

We can write out the distribution of ${\mathbb{P}}_{}(\beta |\mathit{Y})={\mathbb{P}}_{}(\beta ,\mathit{Y})/{\mathbb{P}}_{}\left(\mathit{Y}\right)$, gather terms that depend on $\beta$, and read off the mean and covariance:

$$\begin{array}{rl}\log {\mathbb{P}}_{}(\beta ,\mathit{Y})=& -\frac{1}{2}{\sigma }^{-2}(\mathit{Y}-\mathit{X}\beta {)}^{⊺}(\mathit{Y}-\mathit{X}\beta )-\frac{1}{2}{\sigma }_{\beta }^{-2}{\beta }^{⊺}\beta +\text{Terms that do not depend on }\beta \\ =& -\frac{1}{2}{\sigma }^{-2}{\beta }^{⊺}{\mathit{X}}^{⊺}\mathit{X}\beta +{\sigma }^{-2}{\beta }^{⊺}{\mathit{X}}^{⊺}\mathit{Y}-\frac{1}{2}{\sigma }_{\beta }^{-2}{\beta }^{⊺}\beta +\text{Terms that do not depend on }\beta \\ =& -\frac{1}{2}{\beta }^{⊺}({\sigma }^{-2}{\mathit{X}}^{⊺}\mathit{X}+{\sigma }_{\beta }^{-2}{\mathit{I}}_{})\beta +{\sigma }^{-2}{\beta }^{⊺}{\mathit{X}}^{⊺}\mathit{Y}+\text{Terms that do not depend on }\beta .\end{array}$$

From this, we can read off $\mathrm{\Sigma }$ and $\mu$, and get the above expression.

If we take $\lambda ={\sigma }^2/{\sigma }_{\beta }^2$, then we can see that

$${\mathbb{E}}_{}[\beta |\mathit{Y}]={({\mathit{X}}^{⊺}\mathit{X}+\lambda {\mathit{I}}_{})}^{-1}{\mathit{X}}^{⊺}\mathit{Y}=\hat{\beta }(\lambda ).$$

This gives the ridge procedure some interpretability. First of all, the use of the ridge penalty corresponds to a prior belief that $\beta$ is not too large.

Second, the ridge penalty you use reflects the relative scale of the noise variance and prior variance in a way that makes sense:

• If $\sigma \gg {\sigma }_{\beta }$, then the data is noisy (relative to our prior beliefs). We should not take fitting the data too seriously, and so should estimate a smaller $\beta$ than $\hat{\beta }$. And indeed, in this case $\lambda$ is large, a large $\lambda$ shrinks the estimated coefficients.
• If ${\sigma }_{\beta }\gg \sigma$, then we find it plausible that $\beta$ is very large (relative to the variability in our data). We should not take our prior beliefs too seriously, and estimate a coefficient that matches $\hat{\beta }$. And indeed, in this case, $\lambda$ is small, and we do not shrink the coefficients much.

Sparse regression with the L1 norm (lasso)

One problem with the L2 solution might be that the solution ${\hat{\mathit{\beta }}}_{L2}(\lambda )$ is still “dense”, meaning that, in general, every entry of it is nonzero, and we still have to invert a $P\times P$ matrix.

For example, consider our highly correlated regressor example from the previous lecture. The ridge regression will still include both regressors, and their coefficient estimates will still be highly negatively correlated, but both will be shrunk towards zero. Maybe it would make more sense to select only one variable to include. Let us try to think of how we can change the penalty term to achieve this.

A “sparse” solution is an estimator $\hat{\mathit{\beta }}$ in which many of the entries are zero — that is, an estimated regression line that does not use many of the available regressors.

In a word — ridge regression estimates are not sparse. Let’s try to derive one that is by changing the penalty.

A very intuitive way to produce a sparse estimate is as follows: $${\hat{\mathit{\beta }}}_{L0}(\lambda ):=\underset{\mathit{\beta }}{\mathrm{argmin}}(RSS(\mathit{\beta })+\lambda \sum _{p}1({\beta }_p\ne 0))\text{(practically difficult)}$$

This finds a tradeoff between the best fit to the data, but with a penalty for using more regressors. This makes sense, but is very difficult to compute. In particular, this objective is very non-convex. Bayesian statisticians do attempt to estimate models with a similar kind of penalty (they are called “spike and slab” models), but they are extremeley computationally intensive and beyond the scope of this course.

A convex approximation to the preceding loss is the L1 or Lasso loss, leading to Lasso or L1 regression:

$${\hat{\mathit{\beta }}}_{L1}(\lambda ):=\underset{\mathit{\beta }}{\mathrm{argmin}}(RSS(\mathit{\beta })+\lambda \sum _p\left|{\beta }_b\right|)=\underset{\mathit{\beta }}{\mathrm{argmin}}(RSS(\mathit{\beta })+\lambda {\Vert \mathit{\beta }\Vert }_1).$$

This loss is convex (beacuse it is the sum of two convex functions), and so is much easier to minimize. Furthermore, as $\lambda$ grows, it does produce sparser and sparser solutions — though it may not be obvious at first.

Just as with the ridge regression, you should standardize variables before appyling the Lasso.

The Lasso produces sparse solutions

One way to see that the Lasso produces sparse solutions is to start with a very large $\lambda$ and see what happens as it is slowly decreased.

Start at $\lambda$ very large, so that ${\hat{\mathit{\beta }}}_{L1}(\lambda )=\mathbf{0}$. If we take small step of size $\delta$ in a particular direction away from zero in entry ${\beta }_p$, then $\lambda {\Vert \hat{\beta }\Vert }_1$ increases by $\delta \lambda$, and the RSS changes by the gradient of the squared error,

$$\delta \sum _{n=1}^N(y_n-\hat{\mathit{\beta }}(\lambda {)}^{⊺}{\mathit{x}}_n){\mathit{x}}_{np}=\delta \sum _{n=1}^N{\hat{\epsilon }}_n{\mathit{x}}_{np}=\delta \sum _{n=1}^N{y}_n{\mathit{x}}_{np}\text{ (because }\hat{\mathit{\beta }}(\lambda )=\mathbf{0}\text{)}.$$

As long as $\left|\sum _{n=1}^N{y}_n{\mathit{x}}_{np}\right|<\lambda$ for all $p$, we cannot improve the loss by moving away from $\mathbf{0}$. Since the loss is convex, that means $\mathbf{0}$ is the minimum.

Eventually, we decrease $\lambda$ until $\sum _{n=1}^N{y}_n{\mathit{x}}_{np}=\lambda$ for some $p$. At that point, ${\beta }_p$ moves away from zero as $\lambda$ decreases, and the ${\hat{\epsilon }}_n$ also change. However, until $\sum _{n=1}^N{\hat{\epsilon }}_n{\mathit{x}}_{nq}=\lambda$ for some other $q$, only ${\beta }_p$ will be nonzero. As $\lambda$ decreases more and more variables tend to get added to the model, until $\lambda =0$, when of course ${\hat{\mathit{\beta }}}_{L1}(0)=\hat{\mathit{\beta }}$, the OLS solution. Along the path, variables may come in and out of the regression in complex ways.

The Lasso as a constrained optimization problem

See figure 6.7 from Gareth et al. (2021) for an interpretation of the Lasso as a constrained optimization problem. The shape of the L1 ball provides an intuitive way to understand the sparsity of the solution compared to ridge.

The Bayesian Lasso is not sparse

In contrast to the ridge case, it is not hard to show that if you take a prior

$${\mathbb{P}}_{}\left(\mathit{\beta }\right)\propto \exp \left(\lambda {\Vert \mathit{\beta }\Vert }_1\right)$$

you do not recover a sparse solution for the posterior mean! The difference is that, for the ridge prior, the posterior remained normal, so that the “maximum a–posteriori” (MAP) estimator was equal to the mean. In the Lasso case, the MAP is sparse, but the mean is not, and the two do not coincide because the posterior is not normal.

A more Bayesian way to produce sparse posteriors is by setting a non–zero probability that ${\beta }_p=0$. These are called “spike and slab priors,” and are beyond the scope of this course.

References

Gareth, J., W. Daniela, H. Trevor, and T. Robert. 2021. An Introduction to Statistical Learning: With Applications in Python. Spinger.

Tibshirani, R. 2023. “Advanced Topics in Statistical Learning: Spring 2023.” https://www.stat.berkeley.edu/~ryantibs/statlearn-s23/.