Bias-Variance Decomposition

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.

How repeated-fit uncertainty affects the decomposition estimates

The last column is there because the bias^2 column is itself an estimate. At each input it squares the difference between the mean of \(B = 800\) fitted predictions and the true mean outcome. That average prediction has its own sampling error: with independent repetitions and finite variance, \(\mathbb{E}\big[(\bar{f}_B(x) - f(x))^2\big] = \text{Bias}(x)^2 + \operatorname{Var}(\hat{f}(x))/B\). The var/B column estimates this contribution using the variance across the 800 fits; because that variance uses ddof=0, the estimate is slightly downward-biased.

At degrees 0 to 3 this estimated contribution is much smaller than the reported squared bias. That comparison helps assess how much estimating the mean prediction adds to squared bias, but it is not a confidence interval for either estimate. At degree 5 both are near zero at this precision, consistent with small bias. At the higher degrees, the comparable sizes of the two columns limit what we can infer about squared bias, as noted above.

The bias^2 and variance columns, before display rounding, still add up exactly to the mean squared distance from the truth, because preds.var(0) uses ddof=0; that identity holds for any finite array and says nothing about whether each term is a good estimate of its population counterpart.

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?

Symptom	Possible explanations	What to check or try
High training error, test error close to it	limited capacity; also label noise, missing features, or an optimizer that did not converge	more capacity, better features, less regularization — after checking the optimization and the noise floor
Low training error, much higher test error	high variance; also a distribution shift or a small, unrepresentative test set	more data, more regularization, less capacity — after confirming the two splits come from the same distribution
Both errors near a reliably estimated noise floor	little room left under this squared-error criterion	check subgroup performance and whether better measurements or additional features could lower the noise floor; a different loss would change the evaluation goal
Test error unstable across seeds or folds	high variance in fitting; also variation in how the folds are composed	compare repeated fits on different training samples against a common evaluation set; varying the seed on a fixed sample separately checks algorithmic randomness; then averaging, ensembling, more data
Errors concentrated in one subgroup	large bias, variance, or noise within that subgroup, hidden by the aggregate	examine the decomposition within the subgroup; this needs enough data from it

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.

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.