A Bayesian perspective on ridge or ‘L2’ regression
Reading
Bayesian inference
Frequentist inference proceeds by constructing a sensible estimator (\(\hat{\boldsymbol{\beta}}\) in this course), making modeling assumptions, and estimating how our estimator varies randomly around the truth under repeated sampling of new datasets.
There’s a conceptually very different way of solving the inference problem known generally as “Bayesian inference,” though a better term would probably be “inverse inference.” A Bayesian approach attempts to directly answer the question: which parameters \(\boldsymbol{\beta}\) are most consistent with the data I saw, given assumptions about the randomness in the data?
Let’s begin with a simple example of a sample mean. Suppose we observe data \(x\), which is a scalar-valued tests of mathematical ability for a single student. Suppose we know that each student has an ability \(\mu\), and each test has some error of the form \[ x\sim \mathcal{N}\left(\mu, \sigma^2\right). \] Suppose furthermore we know that students have a range of abilities in the whole population, and that abilities are distributed as \[ \mu \sim \mathcal{N}\left(0, s^2\right). \] Given that a student got the score \(x\), what’s our best guess about their ability?
Note that among students that get high scores relative to the true distribution of \(\mu\) (i.e. higher than \(2 s\)), then many such scores are due to error in \(x\) rather than an acutal high ability. So it makes sense to “shrink” our guess of \(\mu\) towards what we know the true average ability is.
This can be done formally with the posterior distribution. We can ask: among all students who got a score equal to \(x\), what kinds of \(\mu\) did they have? This is formally the conditional distribution \(p(\mu \vert x)\), and it happens that \[ p\left( \mu \vert x\right) = % \gauss{ % \frac{\sigma^{-2} \x}{\sigma^{-2} + \s^{-2}}, % \frac{1}{\sigma^{-2} + \s^{-2}}} = \mathcal{N}\left( \frac{x}{1 + \frac{\sigma^2}{s^2}}, \frac{1}{1 + \frac{\sigma^2}{s^2}}\right). \]
Note that, as \(\sigma \rightarrow \infty\), our guess becomes \(0\), since our test tells us nothing. And as \(s\rightarrow \infty\), our guess becomes \(x\), because the distribution of students tells us nothing.
Bayesain posterior for regression
Under the following model, we can compute the posterior for Bayesian regression:
\[ \boldsymbol{\beta}\sim \mathcal{N}\left(\boldsymbol{0}, \sigma_\beta^2 \boldsymbol{I}\right) \quad\textrm{and}\quad \boldsymbol{Y}\vert \boldsymbol{\beta}, \boldsymbol{X}\sim \mathcal{N}\left(\boldsymbol{X}\boldsymbol{\beta}, \sigma^2 \boldsymbol{I}\right), \]
then
\[ \boldsymbol{\beta}\vert \boldsymbol{Y}, \boldsymbol{X}\sim \mathcal{N}\left( \left(\boldsymbol{X}^\intercal\boldsymbol{X}+ \frac{\sigma^{2}}{\sigma_\beta^{2}} \boldsymbol{I}\right)^{-1} \boldsymbol{X}^\intercal\boldsymbol{Y}, \sigma^{2}\left( \boldsymbol{X}^\intercal\boldsymbol{X}+ \frac{\sigma^{2}}{\sigma_\beta^{2}} \boldsymbol{I}\right)^{-1} \right). \]
One way to derive this is to recognize that, if \(\boldsymbol{\beta}\sim \mathcal{N}\left(\mu, \Sigma\right)\), then
\[ \log \mathbb{P}_{\,}\left(\boldsymbol{\beta}| \mu, \Sigma\right) = -\frac{1}{2} \beta^\intercal\Sigma^{-1} \beta + \beta^\intercal\Sigma^{-1} \mu + \textrm{Terms that do not depend on }\beta. \]
We can write out the distribution of \(\mathbb{P}_{\,}\left(\beta | \boldsymbol{Y}\right) = \mathbb{P}_{\,}\left(\beta, \boldsymbol{Y}\right) / \mathbb{P}_{\,}\left(\boldsymbol{Y}\right)\), gather terms that depend on \(\beta\), and read off the mean and covariance:
\[ \begin{aligned} \log \mathbb{P}_{\,}\left(\beta, \boldsymbol{Y}\right) ={}& -\frac{1}{2} \sigma^{-2} (\boldsymbol{Y}- \boldsymbol{X}\beta)^\intercal(\boldsymbol{Y}- \boldsymbol{X}\beta) -\frac{1}{2} \sigma_\beta^{-2} \beta^\intercal\beta + \textrm{Terms that do not depend on }\beta \\={}& -\frac{1}{2} \sigma^{-2} \beta^\intercal\boldsymbol{X}^\intercal\boldsymbol{X}\beta + \sigma^{-2} \beta^\intercal\boldsymbol{X}^\intercal\boldsymbol{Y} -\frac{1}{2} \sigma_\beta^{-2} \beta^\intercal\beta + \textrm{Terms that do not depend on }\beta \\={}& -\frac{1}{2} \beta^\intercal\left(\sigma^{-2} \boldsymbol{X}^\intercal\boldsymbol{X}+ \sigma_\beta^{-2} \boldsymbol{I}\right) \beta + \sigma^{-2} \beta^\intercal\boldsymbol{X}^\intercal\boldsymbol{Y}+ \textrm{Terms that do not depend on }\beta. \end{aligned} \]
From this, we can read off \(\Sigma\) and \(\mu\), and get the above expression.
If we take \(\lambda = \sigma^2 / \sigma_\beta^2\), then we can see that
\[ \mathbb{E}_{\vphantom{}}\left[\beta | \boldsymbol{Y}\right] = \left(\boldsymbol{X}^\intercal\boldsymbol{X}+ \lambda \boldsymbol{I}\right)^{-1} \boldsymbol{X}^\intercal\boldsymbol{Y}= \hat{\beta}(\lambda). \]
This posterior mean is known as ridge regression, or L2–penalized regression. The reason for the term L2 is discussed below. Note that \(\boldsymbol{X}^\intercal\boldsymbol{X}+ \lambda \boldsymbol{I}\) is always invertible if \(\lambda > 0\), even if \(\boldsymbol{X}\) is not full–column rank. In this sense, ridge regression can deal safely with colinearity.
The prior corresponds to a belief that the “true” coefficients are not too large. Specifically, 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.
A penalty on epmirical loss
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^\intercal x_n)^2 =: RSS(\beta). \]
We motivated this as an approximation to the expected predicted loss, \(L(\beta) = \mathbb{E}_{\vphantom{}}\left[y_\mathrm{new},x_\mathrm{new}\right]{(y_\mathrm{new}- \beta^\intercal 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 \left\Vert\beta\right\Vert_2^2 = RSS(\beta) + \lambda \sum_{p=1}^P \beta_p^2. \]
This is called “ridge” or L2 regression. The term “L2” is because the penalty \(\left\Vert\beta\right\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 \left\Vert\beta\right\Vert_2^2\) is known as a “regularizer,” since it imposes some “regularity” to the estimate \(\hat{\beta}(\lambda)\). Note that
- As \(\lambda \rightarrow \infty\), then \(\hat{\beta}(\lambda) \rightarrow \boldsymbol{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{aligned} L_{ridge}(\beta, \lambda) ={}& (\boldsymbol{Y}- \boldsymbol{X}\beta)^\intercal(\boldsymbol{Y}- \boldsymbol{X}\beta) + \lambda \beta^\intercal\beta \\={}& \boldsymbol{Y}^\intercal\boldsymbol{Y}-2 \boldsymbol{Y}^\intercal\boldsymbol{X}\beta + \beta^\intercal\boldsymbol{X}^\intercal\boldsymbol{X}\beta+ \lambda \beta^\intercal\beta \\={}& \boldsymbol{Y}^\intercal\boldsymbol{Y}-2 \boldsymbol{Y}^\intercal\boldsymbol{X}\beta + \beta^\intercal\left(\boldsymbol{X}^\intercal\boldsymbol{X}+ \lambda \boldsymbol{I}\right) \beta \quad \Rightarrow \\ \frac{\partial L_{ridge}(\beta, \lambda)}{\partial \beta} ={}& -2 \boldsymbol{X}^\intercal\boldsymbol{Y}+ 2 \left(\boldsymbol{X}^\intercal\boldsymbol{X}+ \lambda \boldsymbol{I}\right) \beta \quad\Rightarrow \\ \hat{\beta}(\lambda) ={}& \left(\boldsymbol{X}^\intercal\boldsymbol{X}+ \lambda \boldsymbol{I}\right)^{-1} \boldsymbol{X}^\intercal\boldsymbol{Y}. \end{aligned} \]
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} = \alpha x_{np}\) instead of \(x_{np}\). As we know, the OLS minimizer \(\hat{\beta}_p' = \hat{\beta}_p / \alpha\) and the fitted value \(\hat{\boldsymbol{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'^2 = \frac{\lambda}{\alpha^2} \hat{\beta}_p^2. \]
That is, we will effectively “punish” large values of \(\hat{\beta}_p'\) 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 < \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 \left\Vert\beta\right\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 := \frac{x_n - \bar{x}_n}{\sqrt{\frac{1}{N} \sum_{n=1}^N(x_n - \bar{x})^2}}. \]
We then run ridge regression on \(\boldsymbol{x}'\) 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.