Bias and Variance in Practice: Ridge Regression on Polynomial Features

The bias–variance decomposition is usually drawn as a cartoon: two smooth curves crossing, with a U-shaped total on top. It is more convincing when you measure it. This post sets up a regression problem small enough that every term in the decomposition can be estimated by resampling, fits ridge regression on polynomial features for forty values of the penalty $λ$ , and compares the measured test error with the sum of bias², variance and noise. They agree to three decimals, which is a reassuring check that the formula is not just a story.

Setup

The target is $f (x) = \sin (2 π x)$ on $[0, 1]$ , observed with Gaussian noise of standard deviation $σ = 0.3$ . Each training set has $n = 30$ points. The model is a degree-12 polynomial, which is deliberately too flexible: thirteen parameters for thirty noisy points means an unregularized fit will chase the noise.

Two implementation details matter more than they look:

Standardize each feature column. $x^{12}$ on $[0, 1]$ has a tiny variance compared to $x$ , so without standardization the penalty would treat the columns completely differently.
Do not penalize the intercept. Centering $y$ and fitting the slope-only problem is equivalent and keeps the maths clean.

import numpy as np

SIGMA, DEG = 0.3, 12
f = lambda x: np.sin(2 * np.pi * x)

def make(n, rng):
    x = rng.uniform(0, 1, n)
    return x, f(x) + SIGMA * rng.normal(size=n)

def features(x, mu=None, sd=None):
    X = np.vander(x, DEG + 1, increasing=True)[:, 1:]   # x, x^2, ..., x^12
    if mu is None:
        mu, sd = X.mean(0), X.std(0)
    return (X - mu) / sd, mu, sd

def ridge_fit(X, y, lam):
    ym = y.mean()                                      # unpenalized intercept via centering
    w = np.linalg.solve(X.T @ X + lam * np.eye(X.shape[1]), X.T @ (y - ym))
    return w, ym

The closed form is $\hat{w} = (X^{⊤} X + λ I)^{- 1} X^{⊤} y$ . Writing $X = U D V^{⊤}$ for the SVD shows what $λ$ does: along each principal direction $v_{j}$ the least-squares coefficient is multiplied by the shrinkage factor $d_{j}^{2} / (d_{j}^{2} + λ)$ , so directions with little data support (small $d_{j}$ ) are shrunk first. The sum of these factors is the effective number of parameters, $df (λ) = \sum_{j} d_{j}^{2} / (d_{j}^{2} + λ)$ , which slides from 12 down to 0 as $λ$ grows.

The decomposition

For a fixed input $x$ , let ${\hat{f}}_{D} (x)$ be the prediction of the model trained on a random training set $D$ and $\bar{f} (x) = E_{D} [{\hat{f}}_{D} (x)]$ its average over training sets. For a fresh noisy label $y = f (x) + ε$ ,

E_{D, ε} [(y - {\hat{f}}_{D} (x))^{2}] = \underset{{bias}^{2}}{\underset{⏟}{(f (x) - \bar{f} (x))^{2}}} + \underset{variance}{\underset{⏟}{E_{D} [({\hat{f}}_{D} (x) - \bar{f} (x))^{2}]}} + \underset{noise}{\underset{⏟}{σ^{2}}} .

Derivation of the three terms

Add and subtract $\bar{f} (x)$ and $f (x)$ inside the square:

y - {\hat{f}}_{D} = ε + (f - \bar{f}) + (\bar{f} - {\hat{f}}_{D}) .

Square and take expectations. The three squared terms give noise, bias² and variance. The cross terms vanish: $ε$ is independent of everything else and has mean zero, and $E_{D} [\bar{f} - {\hat{f}}_{D}] = 0$ by definition of $\bar{f}$ , while $f - \bar{f}$ is a constant with respect to $D$ .

Every term can be estimated by brute force: draw 200 training sets, fit each one for every $λ$ , and look at the spread of predictions on a fixed grid of 200 test inputs.

lams = np.logspace(-6, 2, 41)
x_grid = np.linspace(0, 1, 200)
preds = np.zeros((200, len(lams), len(x_grid)))        # (training set, lambda, test point)
for r in range(200):
    x_tr, y_tr = make(30, rng)
    X_tr, mu, sd = features(x_tr)
    for j, lam in enumerate(lams):
        w, ym = ridge_fit(X_tr, y_tr, lam)
        preds[r, j] = features(x_grid, mu, sd)[0] @ w + ym

fbar = preds.mean(0)                                   # average prediction per lambda
bias2 = ((fbar - f(x_grid)) ** 2).mean(1)
var = preds.var(0).mean(1)

Training and validation error versus lambda, and the bias-variance decomposition — Left: training error (30 points) and validation error (2000 points) as a function of the ridge penalty, averaged over 200 training sets. Right: the measured bias², variance and their sum with the noise floor, compared with the directly measured expected test error. The two coincide.

λ	Train MSE	Validation MSE	Bias²	Variance	Bias² + Var + $σ^{2}$
$10^{- 6}$	0.0638	5.715	0.016	6.156	6.262
$10^{- 4}$	0.0678	1.679	0.007	1.724	1.820
$10^{- 3}$	0.0704	0.649	0.009	0.597	0.695
$10^{- 2}$	0.0750	0.361	0.014	0.277	0.381
$10^{- 1}$	0.0952	0.165	0.025	0.050	0.165
$10^{0}$	0.1428	0.253	0.088	0.077	0.255
$10^{1}$	0.1882	0.238	0.120	0.023	0.233
$10^{2}$	0.3456	0.418	0.289	0.033	0.412

Three things stand out. The training error is below the noise floor of $0.09$ for small $λ$ , which is the signature of fitting noise. The validation error is minimised at $λ = 0.1$ , where variance has dropped by two orders of magnitude while bias² is still small. And the variance is not perfectly monotone: around $λ = 1$ it rises slightly before falling again, because the shrunken fits swing in different ways near the interval edges; the decomposition still adds up exactly, so the non-monotonicity is real rather than a measurement artefact.

What the fits look like

Ridge fits to the same 30 points at three values of lambda — The same 30 training points fitted with three penalties. Left: the near-unregularized fit oscillates between points. Middle: close to the validation optimum. Right: so much shrinkage that the fit cannot follow the sine.

The middle panel is what "a good fit" means in this context: not the curve that passes closest to the training points, but the one that would change least if we had drawn a different thirty points.

The regularization path

Plotting each coefficient against $λ$ makes the shrinkage visible. With standardized features and an unregularized fit, the high-degree coefficients take values near $\pm 10$ with alternating signs, cancelling each other almost perfectly to produce a wiggly curve. As $λ$ increases they collapse first, because the high-degree columns are nearly collinear and correspond to small singular values.

Ridge regularization path for twelve polynomial coefficients — The regularization path: mean standardized coefficient for each polynomial degree (darker lines are higher degrees) as the penalty grows. Large cancelling coefficients disappear first; everything is near zero by λ = 10.

TIP

In a real project, do not pick $λ$ by eye from a plot like this. Standardize the features inside the cross-validation loop (fit the scaler on the training fold only) and let something like scikit-learn's RidgeCV search a logarithmic grid; the shape of the validation curve above is typical and a grid of 20–40 values across six decades is plenty.

Learning curves: the other axis

Bias and variance also trade off against the amount of data. Holding $λ$ fixed and growing the training set from 15 to 400 points gives learning curves. The validation error falls towards the noise floor from above and the training error rises towards it from below; the gap between them is, roughly, the variance term.

Learning curves for weak and strong ridge regularization — Learning curves for a weak (left) and a strong (right) penalty, averaged over 100 draws per training-set size. With 15 points and 13 parameters the weakly regularized model is close to interpolation and its validation error explodes; with 400 points both models approach the noise floor, but the strongly regularized one is held back by bias.

With $λ = 10^{- 4}$ the validation MSE goes from 37.0 at $n = 15$ to 0.106 at $n = 50$ and 0.091 at $n = 400$ , essentially the noise floor. With $λ = 1$ the small- $n$ disaster is milder (23.4 at $n = 15$ , 0.263 at $n = 30$ ) but the curve flattens at 0.111, a bias floor that no amount of data removes. The same regularization that saves you with little data costs you with a lot of it, which is why $λ$ should be re-tuned whenever the dataset size changes.

Takeaways

Bias² and variance are not metaphors; they can be estimated by resampling, and their sum plus noise reproduces the measured test error.
Ridge regression trades variance for bias one direction at a time, ordered by singular value, so it hits the nearly-collinear high-degree columns first.
The optimal penalty depends on $n$ . Learning curves show both the overfitting regime and the bias floor on a single plot.