The Bias-Variance Decomposition

there is exactly one theorem in machine learning that every practitioner rederives on a whiteboard at least once a year, and this is it. 𐃏 the squared-error risk of any learned predictor splits into three non-negative pieces — irreducible noise, squared bias, and variance — and every design decision you make (model class, regularisation strength, \(k\), ensemble size, early stopping) is secretly a transaction between the last two.

motivation#

think of a learning algorithm as a dart thrower. the bullseye is the true regression function at a query point; each dart is the prediction of a model trained on a different random sample of data. two independent failure modes emerge:

• bias: the thrower’s aim is systematically off — the darts cluster around the wrong spot. no amount of extra throws fixes this; the model class simply cannot represent the target.
• variance: the aim wobbles — the darts scatter widely around their own centre. any single throw (a single trained model) may land far from where the thrower is aiming, even if the aim is true.

the decomposition#

setup#

fix a query point \(x_0\). data are generated as

\begin{equation} y = f(x) + \varepsilon, \qquad \mathbb{E}[\varepsilon] = 0, \quad \operatorname{Var}(\varepsilon) = \sigma^2, \end{equation}

with \(\varepsilon\) independent of everything else. a learning algorithm consumes a training set \(\mathcal{D} = \{(x_i, y_i)\}_{i=1}^{n}\) (itself random) and emits a predictor \(\hat{f}_{\mathcal{D}}\). write the average prediction over training sets as

\begin{equation} \bar{f}(x_0) = \mathbb{E}_{\mathcal{D}}\!\left[\hat{f}_{\mathcal{D}}(x_0)\right]. \end{equation}

we care about the expected squared error on a fresh observation \(y_0 = f(x_0) + \varepsilon_0\):

\begin{equation} \operatorname{Err}(x_0) = \mathbb{E}_{\mathcal{D},\,\varepsilon_0}\!\left[\big(y_0 - \hat{f}_{\mathcal{D}}(x_0)\big)^2\right]. \end{equation}

derivation#

two applications of the same trick — add and subtract a mean, watch the cross term die.¹

step 1 — peel off the noise. substitute \(y_0 = f(x_0) + \varepsilon_0\):

\begin{align*} \operatorname{Err}(x_0) &= \mathbb{E}\!\left[\big(\varepsilon_0 + f(x_0) - \hat{f}_{\mathcal{D}}(x_0)\big)^2\right] \\ &= \mathbb{E}[\varepsilon_0^2] + 2\,\mathbb{E}[\varepsilon_0]\,\mathbb{E}\!\left[f(x_0) - \hat{f}_{\mathcal{D}}(x_0)\right] + \mathbb{E}_{\mathcal{D}}\!\left[\big(f(x_0) - \hat{f}_{\mathcal{D}}(x_0)\big)^2\right] \\ &= \sigma^2 + \mathbb{E}_{\mathcal{D}}\!\left[\big(f(x_0) - \hat{f}_{\mathcal{D}}(x_0)\big)^2\right]. \end{align*}

step 2 — split the estimation error. add and subtract \(\bar{f}(x_0)\):

\begin{align*} \mathbb{E}_{\mathcal{D}}\!\left[\big(f(x_0) - \hat{f}_{\mathcal{D}}(x_0)\big)^2\right] &= \mathbb{E}_{\mathcal{D}}\!\left[\big(f(x_0) - \bar{f}(x_0) + \bar{f}(x_0) - \hat{f}_{\mathcal{D}}(x_0)\big)^2\right] \\ &= \big(f(x_0) - \bar{f}(x_0)\big)^2 + \mathbb{E}_{\mathcal{D}}\!\left[\big(\hat{f}_{\mathcal{D}}(x_0) - \bar{f}(x_0)\big)^2\right], \end{align*}

where the cross term is \(2\big(f(x_0)-\bar f(x_0)\big)\,\mathbb{E}_{\mathcal{D}}\big[\bar f(x_0) - \hat f_{\mathcal{D}}(x_0)\big] = 0\) by definition of \(\bar f\). hence

\begin{equation} \boxed{\;\operatorname{Err}(x_0) = \underbrace{\sigma^2}_{\text{noise}} + \underbrace{\big(f(x_0) - \bar{f}(x_0)\big)^2}_{\text{bias}^2} + \underbrace{\mathbb{E}_{\mathcal{D}}\big[\big(\hat{f}_{\mathcal{D}}(x_0) - \bar{f}(x_0)\big)^2\big]}_{\text{variance}}\;} \end{equation}

this is eq. (7.9) of (Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome, 2009). three remarks:

• \(\sigma^2\) is the bayes error at \(x_0\): no learner beats it, however clever.
• bias measures the gap between the truth and the average model — a property of the model class and algorithm, not of any single fit.
• variance measures how much a single fit wobbles around that average — a property of the data-sensitivity of the algorithm.

worked closed forms#

k-nearest neighbours. for a knn regressor with neighbours \(x_{(1)}, \dots, x_{(k)}\) of \(x_0\) (inputs held fixed),

\begin{equation} \operatorname{Err}(x_0) = \sigma^2 + \left(f(x_0) - \frac{1}{k}\sum_{\ell=1}^{k} f(x_{(\ell)})\right)^{2} + \frac{\sigma^2}{k}, \end{equation}

so \(k\) is a literal bias–variance dial: growing \(k\) averages away variance as \(\sigma^2/k\) while dragging in ever-farther neighbours whose \(f\) values pollute the mean — bias creeps up. small \(k\) reverses the trade (Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome, 2009).

linear regression with \(p\) features. the fit at \(x_0\) is \(\hat f(x_0) = x_0^\top(\mathbf{X}^\top \mathbf{X})^{-1}\mathbf{X}^\top \mathbf{y}\), a fixed linear functional \(h(x_0)^\top \mathbf{y}\) of the responses, so its variance is \(\lVert h(x_0)\rVert^2 \sigma^2\); averaged over the training inputs this is \(\tfrac{p}{n}\sigma^2\). model complexity here is the parameter count — variance grows linearly in \(p\), and bias is whatever the best linear approximant to \(f\) leaves behind.

an estimator, not a curve. the decomposition is really a statement about estimators. for the plug-in variance estimator \(\hat\sigma^2 = \tfrac1n \sum_i (X_i - \bar X)^2\), one computes \(\mathbb{E}[\hat\sigma^2] = \tfrac{n-1}{n}\sigma^2\): bias \(-\sigma^2/n\), vanishing at rate \(1/n\) — the mse of any estimator likewise splits as \(\operatorname{mse} = \operatorname{bias}^2 + \operatorname{variance}\) (Wasserman, Larry, 2010).

the trade-off#

sweep model complexity and the three terms move in opposite directions: \(\mathrm{bias}^2\) falls (richer classes approximate better), variance rises (richer classes chase noise), noise stays put. the sum is u-shaped, and the art of supervised learning is to sit at the bottom of the u.

every regularisation device in the standard toolkit is a move along this curve:

bagging deserves its footnote: averaging \(B\) identically-distributed fits with pairwise correlation \(\rho\) leaves variance \(\rho\sigma_{\!f}^2 + \tfrac{1-\rho}{B}\sigma_{\!f}^2\) — the reason random forests decorrelate trees with feature subsampling (Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome, 2009).

a simulation#

monte carlo makes the decomposition tangible: fix \(f(x) = \sin 2\pi x\), draw many training sets, fit polynomials of increasing degree, and measure \(\mathrm{bias}^2\) and variance empirically at a test grid.

import numpy as np

rng = np.random.default_rng(0)
f = lambda x: np.sin(2 * np.pi * x)
n, sigma, trials = 30, 0.4, 400
xg = np.linspace(0.05, 0.95, 50)          # test grid

def experiment(degree):
    preds = np.empty((trials, xg.size))
    for t in range(trials):
        x = rng.uniform(0, 1, n)
        y = f(x) + rng.normal(0, sigma, n)
        coef = np.polyfit(x, y, degree)     # least-squares polynomial fit
        preds[t] = np.polyval(coef, xg)
    fbar = preds.mean(axis=0)               # average model
    bias2 = ((fbar - f(xg)) ** 2).mean()    # squared bias, averaged over grid
    var = preds.var(axis=0).mean()          # variance, averaged over grid
    return bias2, var

print(f"{'deg':>3} {'bias^2':>9} {'var':>9} {'sum+noise':>10}")
for d in (1, 2, 3, 5, 7, 9):
    b2, v = experiment(d)
    print(f"{d:>3} {b2:>9.4f} {v:>9.4f} {b2 + v + sigma**2:>10.4f}")

output — \(\mathrm{bias}^2\) collapses at degree 3 (a cubic can bend like one period of a sine; note degree 2 barely helps, a parabola can’t), variance climbs overall — pausing for a small dip at degree 3, where the fit first stops fighting the data — then blows up violently as high-degree fits on \(n=30\) points start whipping at the grid edges:

deg    bias^2       var  sum+noise
  1    0.1555    0.0257     0.3412
  2    0.1508    0.0473     0.3581
  3    0.0028    0.0243     0.1872
  5    0.0001    0.0394     0.1994
  7    0.0028    1.0079     1.1707
  9    0.0365    9.1342     9.3307

(push the degree past 10 and unscaled np.polyfit additionally becomes numerically ill-conditioned — the vandermonde matrix is nearly singular — which inflates the measured “variance” with floating-point noise. real lesson: centre and scale your features.)

beyond squared loss#

the clean additive split is special to squared error. two caveats worth carrying around:

• 0–1 loss does not decompose additively. for classification the interaction is stranger: on the wrong side of the decision boundary, extra variance can reduce expected error by occasionally flipping a systematically-wrong prediction to the right class. squared-loss intuition transfers qualitatively, not quantitatively (Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome, 2009).
• modern over-parameterised models bend the u. past the interpolation threshold, test error can descend a second time — the “double descent” phenomenon: enormously over-parameterised networks generalise despite fitting noise, because implicit regularisation of the training dynamics selects low-norm interpolants.² the decomposition still holds — it is an identity — but variance no longer explodes monotonically with parameter count.²

the probabilistic view of the same identity — decomposing posterior predictive uncertainty — is developed in ch. 8 of (Deisenroth, Marc Peter and Faisal, A. Aldo and Ong, Cheng Soon, 2020).

see also#

• cross validation — estimating the test error you just decomposed
• regularised regression — buying bias to sell variance, with a dial
• loss functions — why squared error owns this identity
• performance metrics — what to report once you stop decomposing

References#

Deisenroth, Marc Peter and Faisal, A. Aldo and Ong, Cheng Soon (2020). Mathematics for Machine Learning, Cambridge University Press.

Hastie, Trevor and Tibshirani, Robert and Friedman, Jerome (2009). The Elements of Statistical Learning, Springer.

Wasserman, Larry (2010). All of Statistics: A Concise Course in Statistical Inference, Springer.