Bias–variance decomposition

The decision-tree note ended with two observations. A fully grown tree fits its training set perfectly and still loses to a tree of depth 3 on new data; and two trees grown on samples from the same distribution disagree almost everywhere. Both are statements about a learning method, the procedure that turns a training set into a model, rather than about one fitted model, because the training set is random. This note makes them precise. Over the randomness of the training set, the expected squared error of a method splits into three terms: the noise in the targets, the squared bias of the method and its variance. The split explains how every complexity knob of the course works, from $k$ in nearest neighbours to the depth of a tree, and it shows exactly what averaging can and cannot fix, which is the subject of the next two notes.

Notation. The data are pairs $(X, Y)$ from a fixed distribution. $m(x) = E[Y \mid X = x]$ is the regression function and $\sigma^2(x) = \operatorname{Var}(Y \mid X = x)$ the noise at $x$. A training set $D = \{(x_i, y_i)\}_{i=1}^n$ consists of $n$ independent draws, and a learning method maps it to a function $\hat f_D$. Expectations $E_D$ and variances $\operatorname{Var}_D$ are taken over training sets of size $n$. $\bar f(x) = E_D[\hat f_D(x)]$ is the average model: the prediction at $x$ averaged over all training sets the method could have seen.

The decomposition over training sets#

Fix a test point $x$ and draw a new target $Y$ at it, independent of the training set. The error $(Y - \hat f_D(x))^2$ is random through both $Y$ and $D$. Its expectation splits in two steps, each a result of the probability course.

Step 1: separate the noise. Add and subtract $m(x)$:

$$\begin{aligned} E\big[(Y - \hat f_D(x))^2\big] ={}& E\big[(Y - m(x))^2\big] + E_D\big[(m(x) - \hat f_D(x))^2\big] \\ &+ 2\,E\big[(Y - m(x))(m(x) - \hat f_D(x))\big] . \end{aligned}$$

The cross term is zero: $Y$ is independent of $D$, so the expectation of the product is the product of the expectations, and $E[Y - m(x)] = 0$ by the definition of $m$. The first term is $\sigma^2(x)$. This is the best-predictor decomposition of the conditional expectation note, applied to the random function $\hat f_D$ instead of a fixed $g$.

Step 2: split the estimation error. For fixed $x$, $\hat f_D(x)$ is an estimator of the number $m(x)$, random through $D$. Its mean squared error is its squared bias plus its variance, as derived in the estimators note (add and subtract $\bar f(x)$; the cross term has mean zero):

$$E_D\big[(\hat f_D(x) - m(x))^2\big] = \big(\bar f(x) - m(x)\big)^2 + \operatorname{Var}_D\big(\hat f_D(x)\big) .$$

Putting the steps together and averaging over the test point $x$:

$$\underbrace{E_{X, Y, D}\big[(Y - \hat f_D(X))^2\big]}_{\text{expected test error}} = \underbrace{E_X\big[\sigma^2(X)\big]}_{\text{noise}} + \underbrace{E_X\big[(\bar f(X) - m(X))^2\big]}_{\text{bias}^2} + \underbrace{E_X\big[\operatorname{Var}_D \hat f_D(X)\big]}_{\text{variance}} .$$

Three remarks keep the decomposition honest. It holds for squared error; for 0–1 loss there are decompositions with terms of the same meaning, but they do not add up so simply (Domingos, 2000). The bias and the variance belong to the method and the sample size: more data typically lower the variance of every method. The bias can fall too: with more points the $k$ nearest neighbours of a fixed $k$ come closer to $x$, and a fixed ridge penalty $\lambda$ weighs less against a growing $\sum_i x_i^2$. Only the bias of the best member of a model class, its approximation error, does not depend on the sample size. And in practice there is one training set, so the three terms are never observed; they are estimated by simulation, as in this note, or their sum is estimated by validation, as in the validation note.

Estimating the terms by simulation#

When the data distribution is known, the terms are easy to estimate: draw many training sets, fit the method on each, and at each test point take the mean and the variance of the predictions across training sets. For trees of growing depth on the problem of the previous note:

Figure

Bias and variance of regression trees

The bundle panels show what the two words mean. The shallow trees agree with each other and are wrong in the same way, coarsening the peaks of the sine; the deep trees are right on average and each one is wrong in its own way. The curve panel shows the trade-off: from depth 1 to depth 8 the squared bias falls from $0.072$ to $0.001$, the variance rises from $0.043$ to $0.094$ ($0.095$ for the fully grown tree), and their sum with the noise has its minimum at depth 3. The variance of the fully grown tree, $0.095$, is close to the noise variance $\sigma^2 = 0.09$ for a simple reason: a leaf with one training point predicts that point’s target, and the target carries the full noise. This needs dense training points: with 40 points on an interval, the point alone in the leaf of $x$ lies close to $x$, and $m$ barely differs between the two. A fully grown tree in this example therefore has expected test error close to $2\sigma^2$ plus its bias, here $0.187$.

What the decomposition says about each model#

Every model of the course has a knob that moves it along the same curve.

Nearest neighbours. With the training inputs held fixed and only the targets random, the $k$-nearest-neighbour prediction at $x$ is the average of $k$ independent noisy targets, so its variance is exactly $\sigma^2 / k$ (Hastie, Tibshirani and Friedman, 2009, Eq. 7.10): $0.25$, $0.05$ and $0.01$ for $\sigma = 0.5$ and $k = 1, 5, 25$. Its bias is the difference between $m(x)$ and the average of $m$ over the $k$ neighbours, which grows as the neighbours get farther away. Small $k$ is the flexible end, large $k$ the rigid one: on Friedman’s function in 10 dimensions (below), where the training inputs are redrawn with each training set and the neighbours move with them, the variance exceeds $\sigma^2/k$: $k = 1$ has variance $12.24$ against $\sigma^2 = 1$ and $k = 50$ has $0.41$, while the squared bias grows from $4.61$ to $13.00$, and $k = 5$ has the smallest total.

Ridge regression. A toy example with exact numbers: one feature, a fixed design with $\sum_i x_i^2 = 10$, no intercept, true slope $\beta = 1$, noise $\sigma^2 = 1$, prediction at $x_0 = 1$. The ridge slope is $\hat\beta_\lambda = \sum_i x_i y_i / (\sum_i x_i^2 + \lambda)$, so $E \hat\beta_\lambda = \beta \cdot 10/(10 + \lambda)$ and $\operatorname{Var} \hat\beta_\lambda = \sigma^2 \cdot 10/(10 + \lambda)^2$. Least squares ($\lambda = 0$) is unbiased with variance $0.1$. With $\lambda = 1$ the squared bias is $(1/11)^2 = 0.0083$ and the variance $10/121 = 0.0826$, a total of $0.0909$: a little bias buys a larger cut in variance, and the error falls below that of the unbiased estimator. With $\lambda = 5$ the bias takes over: $0.111 + 0.044 = 0.156$. Minimizing $\beta^2\lambda^2 + 10\sigma^2$ over $(10 + \lambda)^2$ gives the best penalty $\lambda = \sigma^2 / \beta^2$, the ratio of noise to signal, which is exactly the penalty of a Gaussian prior of variance $\beta^2$ in the MAP reading of ridge. On Friedman’s function, ridge on all 286 monomials of degree at most 3 shows the full curve: from $\lambda = 0.001$ to $\lambda = 1000$ the variance falls from $9.93$ to $0.39$ while the squared bias rises from $0.75$ to $6.81$, and the total is smallest near $\lambda = 10$.

Trees. Depth is the knob, as in the figure. On Friedman’s function, depth 2 has squared bias $9.79$ and variance $5.36$; a fully grown tree has squared bias $3.08$ and variance $9.34$, far above $\sigma^2 = 1$: in 10 dimensions 200 points are sparse, the training point alone in the leaf of $x$ lies far from $x$ and changes with each training set, so the prediction varies with $m$ at that point as well as with the noise. Deep trees are the textbook example of a low-bias, high-variance method.

Ensembles. Average $B$ models trained by the same randomized procedure. Their average has the same expectation as one of them, so the bias does not change; only the variance can fall. If the models were independent the variance would fall by a factor $B$; models trained on the same data are correlated, and the bagging note shows how much is left. The procedure here is a tree grown on a bootstrap sample, so the bias that averaging keeps is that of a bootstrap tree, not of a tree grown on the full training set: a bootstrap sample holds only about 63% of the distinct points, and its tree is slightly cruder. Fifty fully grown trees on bootstrap samples of the same training set have squared bias $3.61$, a little above the full-sample tree’s $3.08$ for that reason, and variance $1.32$ instead of $9.34$, and the lowest test error of the tree-based methods in the table below; on this smooth function, ridge on cubic features does better still. Boosting, in the gradient boosting note, attacks the other term: it combines high-bias models so that the sum has lower bias.

Method on Friedman’s function ($n = 200$, $\sigma^2 = 1$)	Bias²	Variance	Test error
kNN, $k = 1$ / $5$ / $50$	4.61 / 6.47 / 13.00	12.24 / 2.90 / 0.41	17.85 / 10.36 / 14.41
ridge on cubic features, $\lambda = 0.001$ / $10$ / $1000$	0.75 / 1.97 / 6.81	9.93 / 1.38 / 0.39	11.68 / 4.36 / 8.19
tree, depth 2 / 4 / full	9.79 / 4.74 / 3.08	5.36 / 5.98 / 9.34	16.15 / 11.72 / 13.42
50 bagged full trees	3.61	1.32	5.93

The probability course made the same point for estimators: in the estimators note a biased estimator of a uniform end point beat the unbiased one, and shrinkage beat the sample means. Regularization in every form, a ridge penalty, a minimum leaf size, early stopping, is a deliberate purchase of bias with variance.

The double-descent caveat#

The U-shaped curve is not a law. Belkin et al. (2019) observed that when the capacity of a model class grows past the point where it can fit the training data exactly (the interpolation threshold), the test error, after peaking near that threshold, can fall again, sometimes below the bottom of the U. They showed it for random-feature models fitted by minimum-norm least squares, for neural networks, and for random forests whose capacity is indexed first by the number of leaves and then by the number of fully grown trees. The decomposition itself is not contradicted: it is an identity. What fails is the assumption that variance must keep growing with the number of parameters, because among the many interpolating solutions the fitting procedure picks a particular one, of small norm or averaged over many trees, and that choice controls the variance. A small example in the Code section shows the peak: minimum-norm least squares on the first $p$ of 400 Gaussian features with $n = 40$ points has median excess error $0.26$ at $p = 10$, $23.7$ at the threshold $p = 40$, and $0.78$ at $p = 80$; in this example the second descent does not go back below the first minimum. Double descent matters for overparametrized models such as large networks. Most settings of the knobs of this section stay below the threshold, where the classical curve is the right picture. The flexible ends do not: 1-NN, the fully grown tree and ridge with $\lambda = 0.001$ on 286 features fit the 200 training points exactly, or for ridge almost exactly, and so sit at the threshold, while the 50 bagged full trees lie past it on the capacity axis of Belkin et al., which continues with the number of averaged fully grown trees.

What changed#

	One fitted model	A learning method
Random object	the test point	the test point and the training set
Error	$E[(Y - f(X))^2]$	noise + bias² + variance
Overfitting	train error ≪ test error	variance dominates
Underfitting	both errors high	bias dominates
Complexity knob		$k$, $\lambda$, depth, number of models averaged

Limitations#

• Squared error only. For classification losses the decomposition takes other forms, and bias and variance can interact (a high variance can even help a biased classifier).
• Not observable. The terms need the true $m$ and many training sets; on real data only their sum is estimated, by validation.
• Capacity is not one number. “Number of parameters” does not order models by variance when the fitting procedure selects among many exact fits, as double descent shows.

What you should remember#

• Expected test error over training sets = noise + bias² + variance; the cross terms vanish because the new target is independent of the training set and the average model is the mean of the fitted models.
• Bias: error of the average model. Variance: spread of the fitted models across training sets. Noise: the floor $E[\operatorname{Var}(Y \mid X)]$.
• Flexibility (small $k$, small $\lambda$, deep trees) lowers bias and raises variance; the best setting balances them and depends on $n$.
• kNN variance is $\sigma^2/k$ with the training inputs fixed; the best ridge penalty in one dimension is $\sigma^2/\beta^2$; a fully grown tree has variance about $\sigma^2$ when the training points are dense. With random, sparse inputs both variances are much larger.
• Averaging identically distributed models keeps the bias and lowers the variance; boosting lowers the bias.
• Past the interpolation threshold, test error can fall again (double descent); the decomposition still holds.

Code#

Exact ridge numbers, bias and variance by simulation for nearest neighbours, ridge, trees and bagging, and a double-descent peak:

import numpy as np
from sklearn.ensemble import BaggingRegressor
from sklearn.linear_model import Ridge
from sklearn.neighbors import KNeighborsRegressor
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler
from sklearn.tree import DecisionTreeRegressor

# ridge with one feature, fixed design (sum of x_i^2 = 10), beta = 1, sigma^2 = 1, at x0 = 1: exact bias^2 and variance
S, beta, s2 = 10.0, 1.0, 1.0
for lam in (0, 1, 5):
    b2, var = (beta * lam / (S + lam)) ** 2, s2 * S / (S + lam) ** 2
    print(lam, round(b2, 4), round(var, 4), round(b2 + var, 4))
# 0 0.0 0.1 0.1 | 1 0.0083 0.0826 0.0909 | 5 0.1111 0.0444 0.1556

# Friedman's function in 10 dimensions (5 of them irrelevant), n = 200, noise sigma = 1, 50 training sets
rng = np.random.default_rng(0)
f = lambda X: 10 * np.sin(np.pi * X[:, 0] * X[:, 1]) + 20 * (X[:, 2] - 0.5) ** 2 + 10 * X[:, 3] + 5 * X[:, 4]
n, d, sigma, R = 200, 10, 1.0, 50
Xtest = rng.random((500, d))

def decompose(make):
    """Squared bias and variance at the test points, over R independent training sets."""
    P = []
    for _ in range(R):
        X = rng.random((n, d)); y = f(X) + sigma * rng.normal(size=n)
        P.append(make().fit(X, y).predict(Xtest))
    P = np.array(P)
    return ((P.mean(0) - f(Xtest)) ** 2).mean(), P.var(0).mean()

models = [('kNN, k =', k, lambda k=k: KNeighborsRegressor(k)) for k in (1, 5, 20, 50)]
models += [('ridge on cubic features, lambda =', lam, lambda lam=lam: make_pipeline(PolynomialFeatures(3), StandardScaler(), Ridge(lam))) for lam in (0.001, 1, 10, 100, 1000)]
models += [('tree, depth =', depth, lambda depth=depth: DecisionTreeRegressor(max_depth=depth, random_state=0)) for depth in (2, 4, 6, None)]
models += [('bagged full trees, B =', 50, lambda: BaggingRegressor(DecisionTreeRegressor(), n_estimators=50, random_state=0))]
for name, knob, make in models:
    b2, v = decompose(make)
    print(f'{name} {knob}: bias2 {b2:.2f} variance {v:.2f} total {b2 + v + sigma ** 2:.2f}')
# kNN, k = 1: bias2 4.61 variance 12.24 total 17.85
# kNN, k = 5: bias2 6.47 variance 2.90 total 10.36
# kNN, k = 20: bias2 9.54 variance 0.90 total 11.44
# kNN, k = 50: bias2 13.00 variance 0.41 total 14.41
# ridge on cubic features, lambda = 0.001: bias2 0.75 variance 9.93 total 11.68
# ridge on cubic features, lambda = 1: bias2 1.08 variance 2.41 total 4.50
# ridge on cubic features, lambda = 10: bias2 1.97 variance 1.38 total 4.36
# ridge on cubic features, lambda = 100: bias2 4.30 variance 0.93 total 6.23
# ridge on cubic features, lambda = 1000: bias2 6.81 variance 0.39 total 8.19
# tree, depth = 2: bias2 9.79 variance 5.36 total 16.15
# tree, depth = 4: bias2 4.74 variance 5.98 total 11.72
# tree, depth = 6: bias2 3.56 variance 7.81 total 12.38
# tree, depth = None: bias2 3.08 variance 9.34 total 13.42
# bagged full trees, B = 50: bias2 3.61 variance 1.32 total 5.93

# double descent: minimum-norm least squares on the first p of 400 Gaussian features, n = 40, 200 training sets
n, P = 40, 400
w = rng.normal(size=P) / np.arange(1, P + 1); w /= np.linalg.norm(w)      # true weights, decaying, norm 1
for p in (10, 30, 40, 50, 80, 400):
    excess = []
    for _ in range(200):
        X = rng.normal(size=(n, P)); y = X @ w + 0.5 * rng.normal(size=n)
        b = np.linalg.pinv(X[:, :p]) @ y                 # least squares for p < n, minimum-norm interpolation for p > n
        excess.append(np.sum((b - w[:p]) ** 2) + np.sum(w[p:] ** 2))   # E(x'b - x'w)^2 for x ~ N(0, I)
    print(p, round(float(np.median(excess)), 2))
# 10 0.26 | 30 0.83 | 40 23.7 | 50 1.41 | 80 0.78 | 400 0.93

Next#

Among the tree methods the table has a clear winner: fifty trees grown on bootstrap samples and averaged keep the low bias of a deep tree and cut its variance by a factor of seven. Ridge on cubic features is better still on this smooth function, a reminder that the best family depends on the data. But the variance did not fall by a factor of fifty, because trees grown on resamples of the same data are correlated. The next note derives how much variance an average of correlated models keeps, and shows how random forests lower the correlation by giving each split only a random subset of the features.

References#

• S. Geman, E. Bienenstock, R. Doursat. Neural Networks and the Bias/Variance Dilemma. Neural Computation, 4(1), 1992.
• T. Hastie, R. Tibshirani, J. Friedman. The Elements of Statistical Learning, 2nd ed. Springer, 2009. Section 7.3 (the bias–variance decomposition; Eq. 7.10 for nearest neighbours, Eq. 7.11 for linear fits and ridge).
• P. Domingos. A Unified Bias-Variance Decomposition and its Applications. Proceedings of the 17th International Conference on Machine Learning (ICML), 2000.
• M. Belkin, D. Hsu, S. Ma, S. Mandal. Reconciling modern machine-learning practice and the classical bias–variance trade-off. PNAS, 116(32), 2019 (arXiv:1812.11118). Section 2 (random features and neural networks), Section 3 (decision trees and ensembles).
• J. H. Friedman. Multivariate Adaptive Regression Splines. Annals of Statistics, 19(1), 1991. The test function used here (Friedman’s first function, make_friedman1 in scikit-learn).
• E. Sokolov et al. Машинное обучение (HSE). Lecture notes, lecture 10: bagging, random forests and the bias–variance decomposition; practical session 10: the bias–variance decomposition.