The Bias-Variance Decomposition, Derived and Measured

The bias–variance decomposition separates expected squared prediction error into squared bias, variance, and irreducible noise. This article derives the identity, estimates its terms in simulations, and examines how averaging and model size affect them — including a case where the familiar U-shaped curve does not appear.

Suppose that, across repeated training sets, the fitted model predicts 8 half the time and 12 half the time at one fixed input, where the average outcome is 9. The average prediction is 10, so the bias is 1 and its square is 1. The predictions spread \(\pm 2\) around their own average, so the variance is 4. If the observation noise has variance 1, the expected squared prediction error is \(1 + 4 + 1 = 6\). The variance here is the spread of the prediction at this one input across training sets; it says nothing about how predictions differ from one input to another.

The derivation

Write \(f(x) = \mathbb{E}[Y \mid X = x]\) for the best possible prediction at \(x\) under squared-error loss, and \(Y = f(x) + \varepsilon\) with \(\mathbb{E}[\varepsilon \mid X = x] = 0\). Let \(\hat{f}\) be the model fitted on a random training set, so \(\hat{f}\) is itself random, and let the test observation be independent of that training set. Assume the second moments involved are finite. Then the expected squared error at a fixed point \(x\), taken over both the test noise and the choice of training set (and any independent randomness used in fitting), decomposes as:

\(\mathbb{E}\big[(Y-\hat{f}(x))^2 \mid X=x\big] = \underbrace{\big(\mathbb{E}[\hat{f}(x)]-f(x)\big)^2}_{\text{bias}^2} + \underbrace{\mathbb{E}\big[(\hat{f}(x)-\mathbb{E}[\hat{f}(x)])^2\big]}_{\text{variance}} + \underbrace{\sigma^2}_{\text{irreducible}}\)

Here \(\sigma^2 = \operatorname{Var}(Y \mid X = x)\); writing it as a single constant assumes the noise variance is the same at every \(x\). If it is not, the last term is simply \(\operatorname{Var}(Y \mid X = x)\) and the rest of the identity is unchanged.

Write \(\mu(x) = \mathbb{E}[\hat{f}(x)]\) for the average prediction over training sets. The error splits into three pieces:

\(Y – \hat{f}(x) = \varepsilon + \big[f(x) – \mu(x)\big] + \big[\mu(x) – \hat{f}(x)\big]\)

Squaring and taking expectations, the three squared terms give \(\sigma^2\), \(\text{bias}^2\), and variance. The cross terms vanish: \(\varepsilon\) is independent of the training set and has mean zero, so its products with the other two terms in the expansion have mean zero; and \(\mu(x) – \hat{f}(x)\) has mean zero by the definition of \(\mu\), while \(f(x) – \mu(x)\) is a constant, so their product averages to zero too.

Bias describes the average prediction produced by a learning procedure. It depends on the model class, the fitting method, the sample size, and the data distribution — two procedures that search the same class can have different bias, as when one takes the sample mean and the other shrinks it toward zero. Variance is the spread of the fitted model around its own average, a property of how sensitive fitting is to the sample. The third term is the noise, which no model can reduce given the features and data-generating process at hand; different or more accurately measured features can change it.

Both bias and variance are defined by an expectation over training sets. A single fitted model does not reveal these quantities. Estimating them from observed data requires repeated fitting and additional assumptions, particularly because the true regression function is unknown.

Measuring it

With simulated data the expectations can be approximated by brute force: fit the same model on many independently drawn training sets and look at the spread of the predictions. In the code below, preds is an array whose rows are repeated fits and whose columns are fixed input positions on a grid of 200 points. Each column is one instance of the fixed-\(x\) decomposition above; the reported numbers average those columns over the grid, a convenient summary, though not an exact integral over the input distribution.

import numpy as np

rng = np.random.default_rng(0)
NOISE = 0.3
x_grid = np.linspace(-1, 1, 200)
truth = np.sin(3 * x_grid)

def decompose(degree, n=25, reps=800):
    preds = np.empty((reps, len(x_grid)))
    for r in range(reps):
        xs = rng.uniform(-1, 1, n)
        ys = np.sin(3 * xs) + NOISE * rng.normal(size=n)
        preds[r] = np.polyval(np.polyfit(xs, ys, degree), x_grid)
    bias2 = np.mean((preds.mean(0) - truth) ** 2)
    variance = np.mean(preds.var(0))
    return bias2, variance, NOISE ** 2, variance / reps   # last: var/B, see text

print(f"{'deg':>4} {'bias^2':>14} {'variance':>18} {'noise':>8} {'total':>18} {'var/B':>14}")
for d in (0, 1, 2, 3, 5, 9, 15):
    b, v, nz, se = decompose(d)
    print(f"{d:4d} {b:14.4f} {v:18.4f} {nz:8.4f} {b + v + nz:18.4f} {se:14.4f}")
#  deg         bias^2           variance    noise              total          var/B
#    0         0.5208             0.0252   0.0900             0.6360         0.0000
#    1         0.1686             0.0229   0.0900             0.2816         0.0000
#    2         0.1706             0.0611   0.0900             0.3217         0.0001
#    3         0.0033             0.0230   0.0900             0.1163         0.0000
#    5         0.0000             0.0963   0.0900             0.1864         0.0001
#    9         0.3736           338.2195   0.0900           338.6831         0.4228
#   15  25368426.6909   20056441252.7004   0.0900   20081809679.4813  25070551.5659

Estimated squared bias falls from 0.52 at degree 0 to a value that rounds to zero at degree 5. A degree-5 polynomial can approximate \(\sin(3x)\) closely on this interval. Estimated variance moves up and down at low degrees (0.025, 0.023, 0.061, 0.023) and then grows very large: 338 at degree 9 and \(2 \times 10^{10}\) at degree 15 with only 25 training points.

The bias^2 column is estimated from 800 fits. The var/B column estimates the sampling error contribution to that squared-bias estimate. At degree 9 the two columns are about 0.37 and 0.42, and at degree 15 both are about 25 million. These runs therefore do not support a precise estimate of the true squared bias in those rows.

Among the degrees tested, degree 3 has the lowest estimated total, 0.1163, of which 0.09 is noise. Degree 5 has lower estimated bias but a larger variance term and a worse total. That is the trade in this experiment: past a point, the extra flexibility is spent fitting noise, and the fitted noise differs from sample to sample.

Under the symmetric uniform input distribution, the best population least-squares quadratic approximation to the odd function \(\sin(3x)\) has zero constant and quadratic coefficients, so it is also a linear approximation. This describes the best function in each class. The bias in our table instead concerns the average of models fitted on only 25 observations; that average need not equal the population optimum. In each finite sample, the fitted quadratic coefficient can be nonzero. The degree-2 fits here have higher estimated variance than the degree-1 fits, 0.061 against 0.023, with no improvement in estimated total error. The two squared-bias estimates, 0.1706 and 0.1686, alone do not establish a precise difference between their true biases.

Which one do you have?

An aggregate loss gives a small subgroup relatively little weight, so a model can score well overall while doing poorly there; poor subgroup performance may also reflect limited data, missing features, or greater noise within that group. Aggregate averages can hide subgroup errors — but the decomposition itself holds at each fixed \(x\), so for a subgroup defined by the observed inputs it can be examined by averaging over that subgroup’s inputs. Estimating its terms there requires enough data from that subgroup, or additional assumptions.

Where the U-shaped curve stops applying (further reading)

The familiar picture has bias falling and variance rising as model size grows, with test error rising steadily once the model has enough parameters to fit the training data exactly. The decomposition itself does not predict that shape; it is an identity that holds under its assumptions whatever the curve looks like. What need not hold is the pattern. The experiment below shows one run in which it does not.

This section is harder than the rest of the article and can be skipped on a first pass. To read the code: W is a random feature transform, drawn once for each \(p\) and never trained; the only thing fitted is beta. F_tr has shape (n, p) — one row per training example, one column per random feature. np.linalg.pinv returns the pseudo-inverse of that matrix, which exists even when an ordinary inverse does not; multiplying it by y_tr gives a least-squares solution, and when several coefficient vectors fit the training data exactly, the one with the smallest norm.

import numpy as np

rng = np.random.default_rng(0)
n, d = 40, 20
w = rng.normal(size=d)

def gen(m):
    X = rng.normal(size=(m, d))
    return X, X @ w + 0.3 * rng.normal(size=m)

X_tr, y_tr = gen(n)
X_te, y_te = gen(4000)

for p in (5, 20, 39, 40, 41, 60, 120, 400, 2000):
    W = rng.normal(size=(d, p)) / np.sqrt(d)
    F_tr, F_te = np.tanh(X_tr @ W), np.tanh(X_te @ W)
    beta = np.linalg.pinv(F_tr) @ y_tr         # minimum-norm least-squares solution
    print(f"p={p:5d}  train {np.mean((F_tr @ beta - y_tr) ** 2):8.4f}"
          f"  test {np.mean((F_te @ beta - y_te) ** 2):10.4f}"
          f"  |beta| {np.linalg.norm(beta):8.2f}")
# p=    5  train  11.8870  test    14.5171  |beta|     2.83
# p=   20  train   0.6633  test     9.4984  |beta|    12.44
# p=   39  train   0.0022  test    18.1100  |beta|    23.12
# p=   40  train   0.0000  test   101.1249  |beta|    53.75
# p=   41  train   0.0000  test     9.2991  |beta|    15.61
# p=   60  train   0.0000  test     2.0096  |beta|     6.38
# p=  120  train   0.0000  test     1.1771  |beta|     3.26
# p=  400  train   0.0000  test     0.6155  |beta|     1.59
# p= 2000  train   0.0000  test     0.6516  |beta|     0.67

Test error climbs from 9.50 at \(p = 20\) to 18.11 at \(p = 39\) and peaks at 101.12 at \(p = 40\), where the number of features equals the number of training rows. Past that point it falls, to 9.30 at \(p = 41\) and 0.6155 at \(p = 400\), better by more than an order of magnitude than anything below the threshold, and then rises slightly again to 0.6516 at \(p = 2000\). This run shows the double-descent shape. The design limits what it can show. It uses one training set and one draw of W per \(p\), so it measures test error along a single path and says nothing about variance across repeated training sets. And each \(p\) draws a fresh W, so the models are not nested: a larger \(p\) is a different model, not the smaller one with features added.

The coefficient norm is the proposed explanation. At \(p = 40\), when the feature matrix is invertible, there is exactly one coefficient vector that fits the training data, and here it is extreme: \(\|\beta\| = 53.75\) against 12.44 at \(p = 20\). Above the threshold, provided the feature matrix has full row rank, infinitely many vectors fit the data exactly, and the pseudo-inverse selects the one with the smallest norm, which falls to 0.67 at \(p = 2000\). Below the threshold the same call returns the ordinary least-squares fit, which does not interpolate in these runs. Numerically, pinv discards singular values below a tolerance, so these solution descriptions hold up to its numerical rank decision and floating-point error. As a mechanism the norm is suggestive, and no more than that: a raw coefficient norm depends on how many features there are and how they are scaled (duplicating every feature lets the same function be written with a smaller norm), so a falling norm does not by itself establish why the test error falls.

The decomposition is not contradicted here; it remains an identity under its assumptions. What this run shows is that test error does not follow a single U-shaped curve in model size. The usual account of why is that above the threshold the minimum-norm choice acts as an implicit regularizer and the variance falls as capacity grows — a claim about repeated fits, which the single-path design above leaves untested.

Stopping a search after the first rise in validation error can miss a later decrease. If the search budget allows it, compare a wider range of model sizes on validation data, reserving the test set for the final evaluation. An overparameterized model is not automatically overfitted: what matters is the effective complexity that the fitting procedure selects, and the parameter count alone does not measure it.

Exercises

1. Bagging mainly reduces variance in this example. Average decision trees fitted on bootstrap samples and decompose the error at 1, 5, 25, and 100 members, alongside a single tree fitted on the original sample. Report which of the two terms moves.

You should get: a variance term that falls and a bias term that stays small, with the reduction far short of the ideal factor.

import numpy as np
from sklearn.tree import DecisionTreeRegressor

rng = np.random.default_rng(0)
x_grid = np.linspace(-1, 1, 200).reshape(-1, 1)
truth = np.sin(3 * x_grid.ravel())

NOISE = 0.3

def decompose(bag, bootstrap=True, n=60, reps=300):
    preds = np.empty((reps, len(x_grid)))
    for r in range(reps):
        xs = rng.uniform(-1, 1, (n, 1))
        ys = np.sin(3 * xs.ravel()) + NOISE * rng.normal(size=n)
        members = []
        for _ in range(bag):
            idx = rng.integers(0, n, n) if bootstrap else np.arange(n)
            t = DecisionTreeRegressor(random_state=0).fit(xs[idx], ys[idx])
            members.append(t.predict(x_grid))
        preds[r] = np.mean(members, axis=0)
    return np.mean((preds.mean(0) - truth) ** 2), np.mean(preds.var(0))

b, v = decompose(1, bootstrap=False)
print(f"single tree, original data  bias^2 {b:.4f}  variance {v:.4f}"
      f"  bias^2+var {b + v:.4f}  +noise {b + v + NOISE**2:.4f}")
for bag in (1, 5, 25, 100):
    b, v = decompose(bag)
    print(f"bootstrap members {bag:4d}       bias^2 {b:.4f}  variance {v:.4f}"
          f"  bias^2+var {b + v:.4f}  +noise {b + v + NOISE**2:.4f}")
# single tree, original data  bias^2 0.0004  variance 0.0939  bias^2+var 0.0943  +noise 0.1843
# bootstrap members    1       bias^2 0.0008  variance 0.0945  bias^2+var 0.0953  +noise 0.1853
# bootstrap members    5       bias^2 0.0005  variance 0.0542  bias^2+var 0.0546  +noise 0.1446
# bootstrap members   25       bias^2 0.0007  variance 0.0457  bias^2+var 0.0464  +noise 0.1364
# bootstrap members  100       bias^2 0.0008  variance 0.0442  bias^2+var 0.0450  +noise 0.1350

2. The decomposition depends on the loss. For a binary classification, take voters that are each correct with the same probability \(p\), independently of one another on every example, and report the accuracy of their majority vote at 1, 11, and 101 members, for \(p\) below and above one half.

You should get: averaging that helps rapidly on one side of 0.5 and hurts rapidly on the other.

import numpy as np

rng = np.random.default_rng(0)
n_examples = 100_000

for p in (0.45, 0.51, 0.55, 0.60):
    for m in (1, 11, 101):
        votes = rng.random((n_examples, m)) < p     # each voter correct w.p. p
        print(f"p={p:.2f}  m={m:4d}  ensemble accuracy {(votes.sum(1) > m / 2).mean():.4f}")
    print()
# p=0.45  m=   1  ensemble accuracy 0.4502
# p=0.45  m=  11  ensemble accuracy 0.3656
# p=0.45  m= 101  ensemble accuracy 0.1565
#
# p=0.51  m=   1  ensemble accuracy 0.5082
# p=0.51  m=  11  ensemble accuracy 0.5275
# p=0.51  m= 101  ensemble accuracy 0.5797
#
# p=0.55  m=   1  ensemble accuracy 0.5506
# p=0.55  m=  11  ensemble accuracy 0.6347
# p=0.55  m= 101  ensemble accuracy 0.8444
#
# p=0.60  m=   1  ensemble accuracy 0.6009
# p=0.60  m=  11  ensemble accuracy 0.7563
# p=0.60  m= 101  ensemble accuracy 0.9790

3. More data changes the answer. Repeat the degree sweep at \(n = 25\) and \(n = 400\) and report the degree that minimizes the error against \(f(x)\), \(\text{bias}^2 + \text{variance}\), in each case.

You should get: a best-performing tested degree that moves upward in this simulation as the training set grows, and, at low degrees, a bias term that changes little.

import numpy as np

rng = np.random.default_rng(0)
x_grid = np.linspace(-1, 1, 200)
truth = np.sin(3 * x_grid)

def decompose(degree, n, reps=400):
    preds = np.empty((reps, len(x_grid)))
    for r in range(reps):
        xs = rng.uniform(-1, 1, n)
        ys = np.sin(3 * xs) + 0.3 * rng.normal(size=n)
        preds[r] = np.polyval(np.polyfit(xs, ys, degree), x_grid)
    return np.mean((preds.mean(0) - truth) ** 2), np.mean(preds.var(0))

for n in (25, 400):
    print(f"n = {n}")
    for d in (1, 3, 5, 7, 9, 11):
        b, v = decompose(d, n)
        print(f"  degree {d:2d}  bias^2 {b:10.4f}  variance {v:14.4f}"
              f"  bias^2+var {b + v:14.4f}  +noise {b + v + 0.09:14.4f}")
# n = 25
#   degree  1  bias^2     0.1683  variance         0.0245  bias^2+var         0.1928  +noise         0.2828
#   degree  3  bias^2     0.0035  variance         0.0239  bias^2+var         0.0273  +noise         0.1173
#   degree  5  bias^2     0.0008  variance         0.1002  bias^2+var         0.1010  +noise         0.1910
#   degree  7  bias^2     0.0072  variance         6.0503  bias^2+var         6.0575  +noise         6.1475
#   degree  9  bias^2   127.8970  variance     52636.8896  bias^2+var     52764.7866  +noise     52764.8766
#   degree 11  bias^2    47.0277  variance     21032.4639  bias^2+var     21079.4916  +noise     21079.5816
# n = 400
#   degree  1  bias^2     0.1681  variance         0.0015  bias^2+var         0.1696  +noise         0.2596
#   degree  3  bias^2     0.0031  variance         0.0010  bias^2+var         0.0041  +noise         0.0941
#   degree  5  bias^2     0.0000  variance         0.0014  bias^2+var         0.0014  +noise         0.0914
#   degree  7  bias^2     0.0000  variance         0.0019  bias^2+var         0.0019  +noise         0.0919
#   degree  9  bias^2     0.0000  variance         0.0025  bias^2+var         0.0025  +noise         0.0925
#   degree 11  bias^2     0.0000  variance         0.0030  bias^2+var         0.0031  +noise         0.0931

References

• Geman, S., Bienenstock, E., & Doursat, R. (1992). Neural Networks and the Bias/Variance Dilemma. Neural Computation.
• Belkin, M., Hsu, D., Ma, S., & Mandal, S. (2019). Reconciling modern machine-learning practice and the classical bias-variance trade-off. PNAS.