Bias-Variance Tradeoff

In the previous page, we developed intuition for why Ridge regression's shrinkage can improve predictions despite introducing bias. Now we make this precise: we'll derive exact expressions for bias and variance of Ridge estimators, analyze how they depend on $\lambda$, and understand when the tradeoff is most favorable.

This analysis answers the fundamental question: How do we know that trading bias for variance reduction is worthwhile?

To derive bias and variance precisely, we need a clear probabilistic framework.

The linear model:

$$\mathbf{y} = \mathbf{X}\boldsymbol{\beta}^* + \boldsymbol{\epsilon}$$

where:

• $\boldsymbol{\beta}^* \in \mathbb{R}^p$ is the true (unknown) coefficient vector
• $\mathbf{X} \in \mathbb{R}^{n \times p}$ is the design matrix (treated as fixed/non-random)
• $\boldsymbol{\epsilon} \sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I}_n)$ is i.i.d. Gaussian noise

Key quantities:

Let $\mathbf{X}^T\mathbf{X}$ have eigenvalue decomposition: $$\mathbf{X}^T\mathbf{X} = \mathbf{V}\mathbf{D}\mathbf{V}^T = \mathbf{V} \cdot \text{diag}(d_1, \ldots, d_p) \cdot \mathbf{V}^T$$

Define:

• Shrinkage factors: $s_j(\lambda) = \frac{d_j}{d_j + \lambda}$
• Rotated true coefficients: $\boldsymbol{\gamma}^* = \mathbf{V}^T \boldsymbol{\beta}^*$ (coefficients in principal component basis)
• Ridge estimator: $\hat{\boldsymbol{\beta}}_{\lambda} = (\mathbf{X}^T\mathbf{X} + \lambda\mathbf{I})^{-1}\mathbf{X}^T\mathbf{y}$

Definition of bias:

The bias of an estimator is the difference between its expected value and the true parameter:

$$\text{Bias}(\hat{\boldsymbol{\beta}}\lambda) = \mathbb{E}[\hat{\boldsymbol{\beta}}\lambda] - \boldsymbol{\beta}^*$$

Computing the expected value:

\begin{align} \mathbb{E}[\hat{\boldsymbol{\beta}}_\lambda] &= \mathbb{E}[(\mathbf{X}^T\mathbf{X} + \lambda\mathbf{I})^{-1}\mathbf{X}^T\mathbf{y}] &= (\mathbf{X}^T\mathbf{X} + \lambda\mathbf{I})^{-1}\mathbf{X}^T\mathbb{E}[\mathbf{y}] &= (\mathbf{X}^T\mathbf{X} + \lambda\mathbf{I})^{-1}\mathbf{X}^T\mathbf{X}\boldsymbol{\beta}^* \end{align}

The last step uses $\mathbb{E}[\mathbf{y}] = \mathbf{X}\boldsymbol{\beta}^*$ (noise has zero mean).

The bias formula:

\begin{align} \text{Bias}(\hat{\boldsymbol{\beta}}_\lambda) &= (\mathbf{X}^T\mathbf{X} + \lambda\mathbf{I})^{-1}\mathbf{X}^T\mathbf{X}\boldsymbol{\beta}^* - \boldsymbol{\beta}^* &= \left[(\mathbf{X}^T\mathbf{X} + \lambda\mathbf{I})^{-1}\mathbf{X}^T\mathbf{X} - \mathbf{I}\right]\boldsymbol{\beta}^* &= -\lambda(\mathbf{X}^T\mathbf{X} + \lambda\mathbf{I})^{-1}\boldsymbol{\beta}^* \end{align}

The last equality uses the identity: $$(\mathbf{A} + \lambda\mathbf{I})^{-1}\mathbf{A} = \mathbf{I} - \lambda(\mathbf{A} + \lambda\mathbf{I})^{-1}$$

Boxed result:

$$\boxed{\text{Bias}(\hat{\boldsymbol{\beta}}_\lambda) = -\lambda(\mathbf{X}^T\mathbf{X} + \lambda\mathbf{I})^{-1}\boldsymbol{\beta}^*}$$

Eigenvalue decomposition form:

Using $\mathbf{X}^T\mathbf{X} = \mathbf{V}\mathbf{D}\mathbf{V}^T$ and $\boldsymbol{\gamma}^* = \mathbf{V}^T\boldsymbol{\beta}^*$:

$$\text{Bias}(\hat{\boldsymbol{\beta}}_\lambda) = -\mathbf{V}\cdot\text{diag}\left(\frac{\lambda}{d_1 + \lambda}, \ldots, \frac{\lambda}{d_p + \lambda}\right)\cdot\boldsymbol{\gamma}^*$$

Squared bias (summed over coordinates):

$$|\text{Bias}|2^2 = \sum{j=1}^{p} \left(\frac{\lambda}{d_j + \lambda}\right)^2 (\gamma_j^)^2 = \lambda^2 \sum_{j=1}^{p} \frac{(\gamma_j^)^2}{(d_j + \lambda)^2}$$

This is the total squared bias across all coefficient components.

Definition of variance:

The variance (covariance matrix) of the estimator measures how it fluctuates across different realizations of the noise:

$$\text{Var}(\hat{\boldsymbol{\beta}}\lambda) = \mathbb{E}\left[(\hat{\boldsymbol{\beta}}\lambda - \mathbb{E}[\hat{\boldsymbol{\beta}}\lambda])(\hat{\boldsymbol{\beta}}\lambda - \mathbb{E}[\hat{\boldsymbol{\beta}}_\lambda])^T\right]$$

Computing the variance:

Since $\hat{\boldsymbol{\beta}}_\lambda$ is linear in $\mathbf{y}$:

$$\hat{\boldsymbol{\beta}}_\lambda = \mathbf{W}\mathbf{y}, \quad \text{where } \mathbf{W} = (\mathbf{X}^T\mathbf{X} + \lambda\mathbf{I})^{-1}\mathbf{X}^T$$

Using $\text{Var}(\mathbf{W}\mathbf{y}) = \mathbf{W}\text{Var}(\mathbf{y})\mathbf{W}^T$ and $\text{Var}(\mathbf{y}) = \sigma^2\mathbf{I}$:

$$\text{Var}(\hat{\boldsymbol{\beta}}_\lambda) = \sigma^2 \mathbf{W}\mathbf{W}^T = \sigma^2 (\mathbf{X}^T\mathbf{X} + \lambda\mathbf{I})^{-1}\mathbf{X}^T\mathbf{X}(\mathbf{X}^T\mathbf{X} + \lambda\mathbf{I})^{-1}$$

Eigenvalue decomposition form:

Substituting $\mathbf{X}^T\mathbf{X} = \mathbf{V}\mathbf{D}\mathbf{V}^T$:

$$\text{Var}(\hat{\boldsymbol{\beta}}_\lambda) = \sigma^2 \mathbf{V} \cdot \text{diag}\left(\frac{d_1}{(d_1 + \lambda)^2}, \ldots, \frac{d_p}{(d_p + \lambda)^2}\right) \cdot \mathbf{V}^T$$

Total variance (trace):

$$\text{tr}(\text{Var}(\hat{\boldsymbol{\beta}}\lambda)) = \sigma^2 \sum{j=1}^{p} \frac{d_j}{(d_j + \lambda)^2}$$

Boxed result:

$$\boxed{\text{Total Variance} = \sigma^2 \sum_{j=1}^{p} \frac{d_j}{(d_j + \lambda)^2}}$$

Comparison with OLS variance:

For OLS ($\lambda = 0$):

$$\text{Var}(\hat{\boldsymbol{\beta}}_{\text{OLS}}) = \sigma^2 (\mathbf{X}^T\mathbf{X})^{-1} = \sigma^2 \mathbf{V} \cdot \text{diag}\left(\frac{1}{d_1}, \ldots, \frac{1}{d_p}\right) \cdot \mathbf{V}^T$$

Total OLS variance: $\sigma^2 \sum_{j=1}^{p} \frac{1}{d_j}$

Variance reduction factor:

For each principal direction $j$:

$$\frac{\text{Var}{\text{Ridge}, j}}{\text{Var}{\text{OLS}, j}} = \frac{d_j/(d_j + \lambda)^2}{1/d_j} = \frac{d_j^2}{(d_j + \lambda)^2} = s_j^2$$

Since $s_j < 1$, Ridge variance is strictly less than OLS variance in every direction when $\lambda > 0$.

MSE decomposition:

The Mean Squared Error of the Ridge estimator is:

$$\text{MSE}(\hat{\boldsymbol{\beta}}\lambda) = \mathbb{E}|\hat{\boldsymbol{\beta}}\lambda - \boldsymbol{\beta}^*|_2^2 = |\text{Bias}|_2^2 + \text{tr}(\text{Var})$$

Substituting our results:

$$\text{MSE}(\lambda) = \underbrace{\lambda^2 \sum_{j=1}^{p} \frac{(\gamma_j^*)^2}{(d_j + \lambda)^2}}{\text{Squared Bias}} + \underbrace{\sigma^2 \sum{j=1}^{p} \frac{d_j}{(d_j + \lambda)^2}}_{\text{Variance}}$$

This can be written per-component in principal direction $j$:

$$\text{MSE}_j(\lambda) = \frac{\lambda^2 (\gamma_j^*)^2 + \sigma^2 d_j}{(d_j + \lambda)^2}$$

Behavior of MSE with λ:

• At $\lambda = 0$ (OLS): Bias = 0, Variance = $\sigma^2 \sum_j 1/d_j$ (can be huge)
• As $\lambda \to \infty$: Bias² → $|\boldsymbol{\beta}^*|_2^2$ (total signal), Variance → 0
• Optimal $\lambda^*$: Somewhere in between, minimizing total MSE

The tradeoff visualized:

As $\lambda$ increases from 0:

1. Variance decreases monotonically (good)
2. Squared bias increases monotonically (bad)
3. Total MSE initially decreases (variance reduction dominates)
4. Eventually MSE increases (bias becomes too large)
5. Optimal MSE occurs at some $\lambda^* > 0$

Oracle optimal λ:

In principle, the optimal $\lambda^*$ minimizes the MSE:

$$\lambda^* = \arg\min_\lambda \text{MSE}(\lambda)$$

Taking the derivative of MSE and setting to zero yields a complex expression that depends on the true parameters $\boldsymbol{\beta}^*$ and noise variance $\sigma^2$—quantities we don't know in practice.

Simplified case (orthogonal design):

When $\mathbf{X}^T\mathbf{X} = \mathbf{I}$ (orthonormal features), the optimal λ for each coefficient $\beta_j^*$ is:

$$\lambda_j^* = \frac{\sigma^2}{(\beta_j^*)^2}$$

This is the ratio of noise variance to squared signal in that direction.

Per-component analysis:

For a single component with eigenvalue $d$ and true coefficient $\gamma^*$:

$$\text{MSE}(\lambda) = \frac{\lambda^2 (\gamma^*)^2 + \sigma^2 d}{(d + \lambda)^2}$$

Taking derivative and setting to zero:

$$\frac{d}{d\lambda}\text{MSE} = \frac{2\lambda (\gamma^)^2 (d + \lambda)^2 - 2(d + \lambda)[\lambda^2(\gamma^)^2 + \sigma^2 d]}{(d + \lambda)^4} = 0$$

Solving:

$$\lambda^_j = \frac{\sigma^2 d_j}{(\gamma_j^)^2}$$

The challenge: we don't know β:*

The oracle optimal λ depends on $\boldsymbol{\beta}^*$, which is exactly what we're trying to estimate. This creates a circular problem:

• To find optimal λ, we need to know $\boldsymbol{\beta}^*$
• To estimate $\boldsymbol{\beta}^*$, we need to choose λ

Practical solutions include:

1. Cross-validation: Use held-out data to estimate prediction error for different λ values
2. Generalized cross-validation (GCV): An efficient approximation to leave-one-out CV
3. AIC/BIC-type criteria: Use effective degrees of freedom to penalize complexity
4. Empirical Bayes: Estimate λ by maximizing marginal likelihood

We'll explore these methods in detail in the next page.

Existence of improvement:

A fundamental theorem guarantees that Ridge regression (for some $\lambda > 0$) always achieves lower MSE than OLS, except in degenerate cases.

Theorem (Hoerl and Kennard, 1970):

For any true parameter $\boldsymbol{\beta}^ eq \mathbf{0}$ and any design matrix $\mathbf{X}$, there exists $\lambda > 0$ such that:*

$$\text{MSE}(\hat{\boldsymbol{\beta}}\lambda) < \text{MSE}(\hat{\boldsymbol{\beta}}{\text{OLS}})$$

This is because the MSE curve has negative slope at $\lambda = 0$: the marginal variance reduction from a small λ always exceeds the marginal bias increase.

Magnitude of improvement:

The amount of MSE reduction depends on:

1. Condition number of $\mathbf{X}^T\mathbf{X}$: Higher condition number → larger potential improvement
2. Signal-to-noise ratio: Lower SNR → larger potential improvement
3. Dimensionality: More parameters relative to observations → larger improvement
4. Distribution of true coefficients: If true coefficients are small, shrinkage toward zero helps more

Quantitative bound:

For the optimal Ridge estimator:

$$\frac{\text{MSE}(\hat{\boldsymbol{\beta}}{\lambda^*})}{\text{MSE}(\hat{\boldsymbol{\beta}}{\text{OLS}})} \leq \frac{\text{effective df}}{p}$$

where effective df = $\sum_j d_j/(d_j + \lambda^*)$. Since effective df < p, this ratio is always less than 1.

We've focused on estimation MSE: how well we estimate $\boldsymbol{\beta}^*$. But often we care more about prediction: how well we predict new responses.

Prediction risk:

For a new observation $\mathbf{x}{\text{new}}$ with true response $y{\text{new}} = \mathbf{x}{\text{new}}^T\boldsymbol{\beta}^* + \epsilon{\text{new}}$, the prediction risk is:

$$\text{Prediction MSE} = \mathbb{E}[(\hat{y}{\text{new}} - y{\text{new}})^2]$$

where $\hat{y}{\text{new}} = \mathbf{x}{\text{new}}^T \hat{\boldsymbol{\beta}}_\lambda$.

Decomposition:

$$\text{Prediction MSE} = \underbrace{\sigma^2}{\text{Irreducible}} + \underbrace{\mathbf{x}{\text{new}}^T \cdot \text{MSE}(\hat{\boldsymbol{\beta}}\lambda) \cdot \mathbf{x}{\text{new}}}_{\text{Estimation error contribution}}$$

Minimizing estimation MSE also minimizes prediction MSE (for any fixed test point).

In-sample vs. out-of-sample:

• Training error: Typically increases with λ (less fit to training data)
• Test error: Typically decreases then increases with λ (optimal λ minimizes)
• Generalization gap = Test error - Training error: Typically decreases with λ

Ridge closes the generalization gap by trading some training fit for better generalization.

Expected out-of-sample error:

For random test points $\mathbf{x}_{\text{new}}$ with the same distribution as training:

$$\mathbb{E}[\text{Prediction MSE}] = \sigma^2 + \text{tr}(\boldsymbol{\Sigma}_x \cdot \text{MSE Matrix})$$

where $\boldsymbol{\Sigma}_x$ is the covariance of the features.

This shows that prediction performance integrates estimation MSE over the feature distribution—directionally weighting by where test points are likely to fall.