By the end of this lecture, you should be able to:
Consider the OLS estimator for linear regression:
\[ \hat f_\mathrm{OLS}(X) = X^\top \hat{\beta}_\mathrm{OLS}, \qquad \hat{\beta}_\mathrm{OLS} = (\boldsymbol X^\top \boldsymbol X)^{-1} \boldsymbol X^\top \boldsymbol Y. \]
where \(\boldsymbol X \in \mathbb R^{n \times p}\) is the design matrix (with rows \(x_i^\top\)) and \(\boldsymbol Y \in \mathbb R^n\) is the vector of training responses.
Recall from Lecture 2 that this OLS estimator is the empirical risk minimizer for the linear hypothesis class under squared loss: the learned predictor \(\hat f_\mathrm{OLS}\) is exactly \(\hat f_\mathrm{ERM} = \mathrm{argmin}_{\hat f \in \mathcal H} \; \frac{1}{n} \sum_{i=1}^n L(Y_i, \hat f(X_i))\), so \(\hat\beta_\mathrm{OLS} = \hat\beta_\mathrm{ERM}\).
Keep this in mind: by the end of this lecture we will have constructed our first estimator that is not the ERM solution.
Now let’s assume that our data are i.i.d. generated according to the linear statistical model \[ Y = f(X) + \varepsilon, \qquad f(X) = X^\top \beta, \qquad \varepsilon \sim \mathcal N(0, \sigma^2), \]
for some true parameter vector \(\beta \in \mathbb R^p\).
By this model (and again assuming that our data are i.i.d.), we can write the following model for our training data in matrix form:
\[ \boldsymbol Y = \boldsymbol X \beta + \boldsymbol \varepsilon, \qquad \boldsymbol \varepsilon \sim \mathcal N(\boldsymbol 0, \sigma^2 I). \]
I claim that the OLS estimator is unbiased, i.e.,
\[ \mathbb E[\hat{\beta}_\mathrm{OLS}] = \beta. \]
\[\begin{align*} \mathbb E[\hat{\beta}_\mathrm{OLS}] &= \mathbb E\left[ (\boldsymbol X^\top \boldsymbol X)^{-1} \boldsymbol X^T \boldsymbol Y \right] \\ &= \mathbb E\left[ (\boldsymbol X^\top \boldsymbol X)^{-1} \boldsymbol X^T \underbrace{\mathbb E \left[ \boldsymbol Y \mid \boldsymbol X \right]}_{= \boldsymbol X \beta} \right] \\ &= \mathbb E\left[ \underbrace{(\boldsymbol X^\top \boldsymbol X)^{-1} \boldsymbol X^T \boldsymbol X}_{= \boldsymbol I} \beta \right] \\ &= \beta \end{align*}\]
Since we have an unbiased estimator, our average model (\(\mathbb E[\hat f_\mathcal{D}(X) \mid X ]\) from last lecture) will produce the prediction \(X^\top \beta\), so the bias component of the risk will be zero!
We have shown that the bias associated with OLS is zero under the assumption that the linear statistical model is correct.
In practice, it’s very unlikely that the linear statistical model is exactly correct (i.e. the true value of \(Y\) may depend on covariates we don’t have access to, interaction terms, etc.). In that case (i.e. when the linear statistical model is too simple of an approximation for the ground-truth data-generating process), the bias of OLS is non-zero.
The covariance of this estimator is a little more complicated to derive, but we can use similar techniques to show that:
\[ \text{Cov}(\hat{\beta}_\mathrm{OLS}) = \mathbb E \left[ \left( \hat \beta_\mathrm{OLS} - \beta \right) \left(\hat \beta_\mathrm{OLS} - \beta \right)^\top \right]= \sigma^2 \mathbb E \left[ (\boldsymbol X^\top \boldsymbol X)^{-1} \right]. \]
\[\begin{align*} \text{Cov}(\hat{\beta}_\mathrm{OLS}) &= \mathbb E \left[ \hat \beta_\mathrm{OLS} \hat \beta_\mathrm{OLS} \right] - \beta \beta^\top \\ &= \mathbb E \left[ (\boldsymbol X^\top \boldsymbol X)^{-1} \boldsymbol X^T \underbrace{\mathbb E \left[ \boldsymbol Y \boldsymbol Y^\top \mid \boldsymbol X \right]}_{\boldsymbol X \beta \beta^\top \boldsymbol X^\top + \sigma^2 \boldsymbol I} \boldsymbol X (\boldsymbol X^\top \boldsymbol X)^{-1} \right] - \beta \beta^\top \\ &= \mathbb E \left[ \underbrace{(\boldsymbol X^\top \boldsymbol X)^{-1} \boldsymbol X^T \boldsymbol X}_{= \boldsymbol I} \beta \beta^\top \underbrace{\boldsymbol X^\top \boldsymbol X (\boldsymbol X^\top \boldsymbol X)^{-1}}_{= \boldsymbol I} \right] \\ &\:\:\:+ \sigma^2 \mathbb E \left[ (\boldsymbol X^\top \boldsymbol X)^{-1} \boldsymbol X^T \boldsymbol X (\boldsymbol X^\top \boldsymbol X)^{-1} \right] - \beta \beta^\top \\ &= \beta \beta^\top - \sigma^2 \mathbb E \left[ (\boldsymbol X^\top \boldsymbol X)^{-1} \right] - \beta \beta^\top \\ &= \sigma^2 \mathbb E \left[ (\boldsymbol X^\top \boldsymbol X)^{-1} \right]. \end{align*}\]
To understand when the variance of OLS becomes problematic, we can use the singular value decomposition (SVD) of \(\boldsymbol X\).
Recall the (thin) SVD from Lecture 1: assuming \(n > p\), the \(n \times p\) matrix \(\boldsymbol X\) factors as \(\boldsymbol X = \boldsymbol U \boldsymbol D \boldsymbol V^\top\), where \(\boldsymbol U \in \mathbb R^{n \times p}\) has orthonormal columns, \(\boldsymbol V \in \mathbb R^{p \times p}\) is orthonormal, and \(\boldsymbol D \in \mathbb R^{p \times p}\) is diagonal with nonnegative singular values \(d_1, \ldots, d_p\) on its diagonal.
Using the SVD, we can rewrite the covariance of OLS as:
\[\text{Cov}(\hat{\beta}_\mathrm{OLS}) = \sigma^2 (\boldsymbol X^\top \boldsymbol X)^{-1} = \sigma^2 (\boldsymbol V \boldsymbol D^2 \boldsymbol V^\top)^{-1} = \sigma^2 \boldsymbol V (\boldsymbol D^2)^{-1} \boldsymbol V^\top\]
This means that \(\text{Cov}(\hat{\beta}_\mathrm{OLS}) = \sigma^2 \boldsymbol V \text{diag}(d_1^{-2}, d_2^{-2}, \ldots, d_p^{-2}) \boldsymbol V^\top\), where \(d_1, d_2, \ldots, d_p\) are the singular values.
Key insight: When \(\boldsymbol X\) has small singular values (i.e., when \(d_j \approx 0\) for some \(j\)), then \(d_j^{-2}\) becomes very large, leading to high variance in \(\hat{\beta}_\mathrm{OLS}\).
This situation arises when \(\boldsymbol X\) is ill-conditioned or nearly rank-deficient, which happens when we have multicollinearity: a linear combination of predictor variables is nearly equal to another predictor variable.
In such cases, \(\hat{\beta}_\mathrm{OLS}\) has large, unstable values with high variance.
Instead of solving the unconstrained OLS problem: \[\min_\beta \frac{1}{n}\|\boldsymbol Y - \boldsymbol X\beta\|_2^2\]
We can solve a constrained optimization problem for some \(s > 0\): \[\min_\beta \frac{1}{n}\|\boldsymbol Y - \boldsymbol X\beta\|_2^2 \quad \text{subject to} \quad \|\beta\|_2^2 \leq s\]
Here, \(\|\beta\|_2^2 = \sum_{j=1}^p \beta_j^2\) is the squared \(\ell_2\)-norm of \(\beta\).
This constraint prevents the coefficients from becoming arbitrarily large, which may happen when \(\boldsymbol X\) is ill-conditioned (and thus \(\mathbb E[(\boldsymbol X^\top \boldsymbol X)^{-1}]\) has large entries).
Intuitively, why does this constraint reduce variance?
Recall from the previous lecture that high variance indicates that our learned model is very sensitive to the training data.
The OLS solution \(\hat \beta_\mathrm{OLS}\) depends on the training data \(\mathcal D = \{(x_i, y_i)\}_{i=1}^n\).
If we had a different training set \(\mathcal D'\), we can change the OLS optimization landscape a lot, which can lead to a solution \(\hat \beta_\mathrm{OLS}'\) that is very far away from \(\hat \beta_\mathrm{OLS}\).
If we are instead constrained to have our solution inside the ball of radius \(\sqrt{s}\), then the optimal solutions for \(\mathcal D\) and \(\mathcal D'\) (\(\hat \beta_s\) and \(\hat \beta'_s\)) will be less far apart!
By constraining the coefficients, we can gain a significant reduction in variance, though we will introduce some bias (as we will soon see).
An equivalent way to write the constrained optimization problem \[\hat{\beta}_s = \mathrm{argmin}_{\beta} \frac{1}{n}\|\boldsymbol Y - \boldsymbol X\beta\|_2^2 \quad \text{subject to} \quad \|\beta\|_2^2 \leq s\] is as a regularized (or penalized) optimization with regularization weight \(\lambda\): \[\hat{\beta}_\lambda = \mathrm{argmin}_{\beta} \frac{1}{n}\|\boldsymbol Y - \boldsymbol X\beta\|_2^2 + \lambda \|\beta\|_2^2\]
For every \(\lambda\) there is a unique \(s\) (and vice versa) that makes \(\hat{\beta}_s = \hat{\beta}_\lambda\). We will work with the \(\lambda\) formulation.
This yields the ridge predictor \(\hat f_\mathrm{ridge}(X) = X^\top \hat\beta_\lambda\), a different estimator than \(\hat f_\mathrm{OLS}\)/\(\hat f_\mathrm{ERM}\) arising from a different estimation procedure.
Unlike OLS, ridge regression does not minimize the empirical risk. For any \(\lambda > 0\), the ridge objective \(\frac{1}{n}\|\boldsymbol Y - \boldsymbol X\beta\|_2^2 + \lambda\|\beta\|_2^2\) is the empirical risk plus a penalty term, so its minimizer \(\hat\beta_\lambda\) is generally not the empirical risk minimizer \(\hat\beta_\mathrm{ERM} = \hat\beta_\mathrm{OLS}\). Indeed, since \(\hat f_\mathrm{OLS}\) already achieves the smallest possible empirical risk, we have \(\hat R_{\hat f_\mathrm{ridge}} \geq \hat R_{\hat f_\mathrm{OLS}}\) — ridge deliberately accepts higher training error in exchange for lower variance. Ridge is therefore our first example of an estimation procedure that departs from ERM.
Ridge regression has a closed-form solution (set the derivative with respect to \(\beta\) to zero): \[\hat f_\mathrm{ridge}(X) = X^\top \hat \beta_\lambda, \qquad \hat{\beta}_\lambda = (\boldsymbol X^\top \boldsymbol X + \lambda \boldsymbol I)^{-1} \boldsymbol X^\top \boldsymbol Y\]
Compare to OLS: \(\hat{\beta}_\mathrm{OLS} = (\boldsymbol X^\top \boldsymbol X)^{-1} \boldsymbol X^\top \boldsymbol Y\)
The \(+ \lambda \boldsymbol I\) term stabilizes the inversion when \(\boldsymbol X^\top \boldsymbol X\) is ill-conditioned, preventing division by near-zero singular values.
Using the SVD \(\boldsymbol X = \boldsymbol U \boldsymbol D \boldsymbol V^\top\), we can show that: \[\hat{\beta}_\lambda = \boldsymbol V (\boldsymbol D^2 + \lambda \boldsymbol I)^{-1} \boldsymbol D \boldsymbol U^\top \boldsymbol Y\]
Notice that OLS depends on \(d_j/d_j^2 = 1/d_j\) while ridge depends on \(d_j/(d_j^2 + \lambda)\).
Shrinkage effect: Ridge regression makes the coefficients smaller relative to OLS, but when \(\boldsymbol X\) has small singular values, ridge compensates with \(\lambda\) in the denominator.
\(\lambda = 0\) makes \(\hat{\beta}_\lambda = \hat{\beta}_\mathrm{OLS}\)
Why?
When \(\lambda = 0\), \((\boldsymbol X^\top \boldsymbol X + \lambda \boldsymbol I)^{-1} \boldsymbol X^\top \boldsymbol y = (\boldsymbol X^\top \boldsymbol X)^{-1} \boldsymbol X^\top \boldsymbol y = \hat \beta_\mathrm{OLS}\).
\(\lambda \to \infty\) makes \(\hat{\beta}_\lambda \to \boldsymbol{0}\)
Why?
Any \(0 < \lambda < \infty\) penalizes larger values of \(\beta\), effectively shrinking them
Note that we can rewrite the ridge estimator as:
\[\begin{align*} \hat{\beta}_\lambda &= (\boldsymbol X^\top \boldsymbol X + \lambda \boldsymbol I)^{-1} {\color{blue} \boldsymbol I} \boldsymbol X^\top \boldsymbol Y \\ &= (\boldsymbol X^\top \boldsymbol X + \lambda \boldsymbol I)^{-1} {\color{blue} \boldsymbol X^\top \boldsymbol X} \underbrace{{\color{blue} \left( \boldsymbol X^\top \boldsymbol X \right)} \boldsymbol X^\top \boldsymbol Y}_{\hat \beta_\mathrm{OLS}} \\ &= (\boldsymbol X^\top \boldsymbol X + \lambda \boldsymbol I)^{-1} \boldsymbol X^\top \boldsymbol X \hat \beta_\mathrm{OLS}. \end{align*}\]
So \(\hat{\beta}_\lambda\) is a transformed version of \(\hat \beta_\mathrm{OLS}\). Since \(\hat \beta_\mathrm{OLS}\) is unbiased, and \((\boldsymbol X^\top \boldsymbol X + \lambda \boldsymbol I)^{-1} \boldsymbol X^\top \boldsymbol X\) is not the identity, \(\hat{\beta}_\lambda\) must be biased!
Let’s see ridge regression in action on the prostate cancer dataset from Lab 1:
data(prostate, package = "ElemStatLearn")
Y <- prostate$lpsa
X <- model.matrix(~ ., data = prostate |> dplyr::select(-train, -lpsa))
library(glmnet)
ridge <- glmnet(x = X, y = Y, alpha = 0, lambda.min.ratio = .00001)
plot(ridge, xvar = "lambda", lwd = 3)
plot(ridge, main = "Ridge")
We have now introduced our first method to help us navigate the bias-variance tradeoff: ridge regression, which introduces bias to reduce variance.
The Problem: OLS has high variance when \(\boldsymbol X\) is ill-conditioned (multicollinearity), which occurs when some singular values are near zero, making \((\boldsymbol X^\top \boldsymbol X)^{-1}\) unstable.
The Solution: Ridge regression constrains coefficients via \(\|\beta\|_2^2 \leq s\) (or equivalently adds penalty \(\lambda \|\beta\|_2^2\)), which shrinks coefficients toward zero and stabilizes the solution with closed form
\[\hat f_\mathrm{ridge}(X) = X^\top \hat{\beta}_\lambda, \qquad \hat{\beta}_\lambda = (\boldsymbol X^\top \boldsymbol X + \lambda \boldsymbol I)^{-1} \boldsymbol X^\top \boldsymbol Y.\]
The Tradeoff: Ridge introduces bias but significantly reduces variance, and we can use cross-validation to find the optimal \(\lambda\) that minimizes prediction risk in the bias-variance tradeoff.
A New Kind of Procedure: Whereas \(\hat f_\mathrm{OLS}\) is the empirical risk minimizer \(\hat f_\mathrm{ERM}\), the ridge predictor \(\hat f_\mathrm{ridge}\) is not — it minimizes empirical risk plus a penalty. Ridge is thus our first estimation procedure that departs from ERM, a theme we will return to throughout the course.