04 - Bias-Variance Tradeoff: The Central Dilemma of ML

A model goes wrong in two opposite ways, and the remedy for one makes the other worse. If it is too rigid, a straight line for data that draw a curve, it is wrong in the same way on any sample you give it: systematic error, in jargon bias. If it is too flexible, a degree-9 polynomial on thirty points, it chases the noise of the sample it has seen and changes shape with every new sample: instability error, in jargon variance. This article gives you a way to tell which of the two is ruining your model and what to do about it: it derives the decomposition of the error, measures it with a simulation in which the truth is known, and uses it to read learning curves, choose the regularisation and understand why averaging several models helps. You need chapter 02 (least squares and Ridge) and an idea of mean and variance.

Everything was run with numpy 2.4.6 and scikit-learn 1.9.0, on synthetic data with a fixed seed and a declared true function, so without downloads. Every output shown is identical across two runs.

The problem in an example

The true function is f(x) = sin(2πx) on [0, 1]; each measurement is y = f(x) + ε with Gaussian noise of standard deviation 0.3, hence variance σ² = 0.09. The noise is the part of the error that no model can remove: even knowing f perfectly, a new measurement would be off on average by 0.09 squared. We take thirty points and fit polynomials of degree 1, 3 and 9. Degree 1 will never manage to follow a sinusoid (underfitting); degree 9 has ten coefficients for thirty points and bends them onto the noise (overfitting): on the sample it has seen it errs little, on a new one a lot.

The decomposition of the error

Imagine repeating the experiment: each time you draw a different training sample D and get a different model f̂_D. At a fixed point x₀ the expected squared error on a new measurement splits into three terms: E[(y − f̂_D(x₀))²] = Var[f̂_D(x₀)] + (E[f̂_D(x₀)] − f(x₀))² + σ². The first is the variance, how much the prediction changes from one sample to another; the second is the squared bias, the distance between the mean prediction and the truth; the third is the irreducible noise. The textbook by James, Witten, Hastie and Tibshirani (An Introduction to Statistical Learning, section 2.2.2) gives the same decomposition. The proof is done by adding and subtracting f(x₀) and E[f̂_D(x₀)] inside the square: the cross products vanish because the noise of the new measurement is independent of the sample, and because the mean of the deviations from the mean is zero.

The decomposition has a practical consequence. With a single dataset you cannot measure bias and variance separately, because you need the truth and several samples. You can, however, reproduce the situation with a simulation where the truth is declared, and see whether the decomposition holds.

Measuring bias and variance with a simulation

The listing generates 500 samples of 30 points, fits a polynomial on each, and compares the 500 resulting curves at 49 points of the grid. The squared bias is the mean distance between the mean curve and f, the variance is the mean dispersion of the curves around their mean. As an independent check, the last computation evaluates the same curves against fresh, noisy labels: if the decomposition is right, the result must coincide with the sum of the three terms.

import numpy as np
from numpy.polynomial import Polynomial

rng = np.random.default_rng(0)
f = lambda x: np.sin(2 * np.pi * x)                 # true function
sigma, n, M = 0.3, 30, 500                          # noise, examples per dataset, number of datasets
xg = np.linspace(0.02, 0.98, 49)                    # points where the error is measured

print(f"irreducible noise sigma^2 = {sigma**2:.4f}")
for degree in (1, 3, 9):
    P = np.empty((M, len(xg)))                      # predictions of M models, each on a different dataset
    for m in range(M):
        x = rng.uniform(0, 1, n)
        y = f(x) + sigma * rng.standard_normal(n)
        P[m] = Polynomial.fit(x, y, degree)(xg)
    bias2 = np.mean((P.mean(axis=0) - f(xg)) ** 2)
    var = np.mean(P.var(axis=0))
    y_new = f(xg) + sigma * rng.standard_normal(P.shape)      # fresh labels, never seen
    test = np.mean((y_new - P) ** 2)
    print(f"degree {degree}: bias^2={bias2:.4f} variance={var:.4f} sum with sigma^2={bias2 + var + sigma**2:.4f}  measured test error={test:.4f}")

# Output:
# irreducible noise sigma^2 = 0.0900
# degree 1: bias^2=0.1820 variance=0.0231 sum with sigma^2=0.2951  measured test error=0.2951
# degree 3: bias^2=0.0037 variance=0.0158 sum with sigma^2=0.1095  measured test error=0.1096
# degree 9: bias^2=0.0220 variance=8.7130 sum with sigma^2=8.8250  measured test error=8.8381

Degree 1 has a squared bias of 0.1820 and a variance of 0.0231: the error is almost entirely systematic, and no additional sample removes it. Degree 3 has bias 0.0037 and variance 0.0158, and the expected error drops to 0.1095, close to the irreducible noise of 0.09. Degree 9 has a small bias (0.0220) but a variance of 8.7130, and the error becomes 8.825. The check column confirms the decomposition: 0.2951 against 0.2951, 0.1095 against 0.1096, 8.825 against 8.838; the small difference in the last case is sampling noise, because with tails this heavy a few extreme curves dominate the mean. The very high variance of degree 9 is not a defect of the measuring method: with 10 coefficients on 30 noisy points the polynomial is poorly constrained between the points.

Learning curves

In practice you do not know the truth, but you can look at how the error on the training set and on new data changes as the number of examples grows. The scikit-learn guide to validation and learning curves (learning_curve, validation_curve) describes the same tool. Here it is computed with the simulation: for each sample size n 200 models are fitted and the median of the training error and of the error on a test grid, noise included, is reported. The median is used because for degree 9 the mean is dominated by a few extreme samples.

import numpy as np
from numpy.polynomial import Polynomial

rng = np.random.default_rng(1)
f = lambda x: np.sin(2 * np.pi * x)
sigma, M = 0.3, 200
xt = np.linspace(0, 1, 400)                          # dense test grid
print("median over 200 datasets; n | degree 1: train  test | degree 9: train  test")
for n in (12, 20, 40, 80, 160, 320):
    row = []
    for degree in (1, 9):
        tr, te = [], []
        for _ in range(M):
            x = rng.uniform(0, 1, n)
            y = f(x) + sigma * rng.standard_normal(n)
            p = Polynomial.fit(x, y, degree)
            tr.append(np.mean((y - p(x)) ** 2))
            te.append(np.mean((f(xt) - p(xt)) ** 2) + sigma ** 2)
        row += [np.median(tr), np.median(te)]
    print(f"{n:4d} | {row[0]:12.4f} {row[1]:6.4f} | {row[2]:12.4f} {row[3]:6.4f}")

# Output:
# median over 200 datasets; n | degree 1: train  test | degree 9: train  test
#   12 |       0.2344 0.3182 |       0.0103 998.4158
#   20 |       0.2388 0.3117 |       0.0397 1.2177
#   40 |       0.2614 0.3001 |       0.0659 0.1297
#   80 |       0.2719 0.2935 |       0.0778 0.1057
#  160 |       0.2797 0.2909 |       0.0845 0.0962
#  320 |       0.2834 0.2891 |       0.0871 0.0927

The reading is the practical rule of diagnosis. With degree 1 the two columns approach each other and stop together at a high level: training error 0.2834 and test error 0.2891 at n = 320, well above the irreducible noise of 0.09. When both curves converge high, the model is too simple and adding data does not help: you need a more flexible model or other measurements. With degree 9 the training error is very low (0.0103 at n = 12, where the polynomial almost interpolates the points) while the test error is enormous (median 998.4); the gap closes only as n grows and at 320 the two errors are 0.0871 and 0.0927, both close to the limit of 0.09. A large gap between training and test indicates variance, and the cure comes from the data (more examples), from a simpler model, or from regularisation.

Worked exercise: how much data does it take to halve the variance?

For a model with a fixed number of parameters the variance of the predictions scales as 1/n: doubling the examples halves it, there is no need to quadruple them. Let us check it on degree 3, with 1000 samples for each n.

import numpy as np
from numpy.polynomial import Polynomial

rng = np.random.default_rng(4)
f = lambda x: np.sin(2 * np.pi * x)
sigma, M = 0.3, 1000
xg = np.linspace(0.02, 0.98, 49)
prev = None
for n in (30, 60, 120, 240):
    P = np.empty((M, len(xg)))
    for m in range(M):
        x = rng.uniform(0, 1, n)
        P[m] = Polynomial.fit(x, f(x) + sigma * rng.standard_normal(n), 3)(xg)
    var = np.mean(P.var(axis=0))
    print(f"n={n:3d}  variance (degree 3) = {var:.5f}" + (f"  ratio to the previous = {prev / var:.2f}" if prev else ""))
    prev = var

# Output:
# n= 30  variance (degree 3) = 0.01686
# n= 60  variance (degree 3) = 0.00657  ratio to the previous = 2.57
# n=120  variance (degree 3) = 0.00310  ratio to the previous = 2.12
# n=240  variance (degree 3) = 0.00156  ratio to the previous = 1.99

The variance goes from 0.01686 to 0.00156 when n goes from 30 to 240, a factor of about 11 for an 8-fold increase in data; the ratio between two successive levels is 2.57, 2.12 and 1.99, so it approaches 2 as n grows, while with little data the variance falls faster than expected. The bias, instead, does not depend on n: halving the variance is not enough if the error is systematic.

Regularisation: choosing where to sit on the curve

Regularisation adds to the least-squares problem a penalty on the size of the weights: the Ridge of chapter 02 minimises ||y − Xw||² + α||w||². A small α leaves the model free (little bias, a lot of variance); a large α stiffens it (more bias, less variance). The listing uses a degree-9 polynomial with standardised measurements, because the penalty treats all the weights in the same way, and measures bias and variance for six values of α.

import numpy as np
from sklearn.linear_model import Ridge
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import PolynomialFeatures, StandardScaler

rng = np.random.default_rng(2)
f = lambda x: np.sin(2 * np.pi * x)
sigma, n, M = 0.3, 30, 300
xg = np.linspace(0.02, 0.98, 49)[:, None]
print("alpha     bias^2  variance  sum+sigma^2")
for alpha in (1e-6, 1e-3, 1e-2, 1e-1, 1.0, 10.0):
    P = np.empty((M, len(xg)))
    for m in range(M):
        x = rng.uniform(0, 1, (n, 1))
        y = f(x[:, 0]) + sigma * rng.standard_normal(n)
        mod = make_pipeline(PolynomialFeatures(9, include_bias=False), StandardScaler(), Ridge(alpha=alpha)).fit(x, y)
        P[m] = mod.predict(xg)
    b2 = np.mean((P.mean(axis=0) - f(xg[:, 0])) ** 2)
    v = np.mean(P.var(axis=0))
    print(f"{alpha:7.0e}  {b2:6.4f}  {v:8.4f}  {b2 + v + sigma ** 2:8.4f}")

# Output:
# alpha     bias^2  variance  sum+sigma^2
#   1e-06  0.0002    0.1588    0.2490
#   1e-03  0.0007    0.1763    0.2671
#   1e-02  0.0031    0.0281    0.1211
#   1e-01  0.0318    0.0320    0.1538
#   1e+00  0.0694    0.0261    0.1855
#   1e+01  0.1306    0.0205    0.2411

The variance drops a lot, from 0.1763 to 0.0281, going from alpha=1e-3 to 1e-2, while the squared bias grows from 0.0007 to 0.0031 and then to 0.1306 for alpha=10. The total error has a minimum at α = 0.01 (0.1211). Two remarks. Even alpha=1e-6 is already a regularisation: compared with the unpenalised degree 9 of the previous simulation (variance 8.713), the variance here is 0.1588, because the penalty and the standardisation stabilise the system. Moreover the minimum obtained (0.1211) is close to, but not better than, the unpenalised degree 3 (0.1095): regularising a model that is too large does not replace the choice of the right complexity. How to choose α without knowing the truth is done with the cross-validation of chapter 05.

Reducing the variance by averaging

The second way against variance is to average several models. Bagging (bootstrap aggregating) trains many models, each on a sample obtained by drawing with replacement from the training set, and averages their predictions. A tree with no depth limit is the extreme case: almost zero bias, high variance. The listing compares a single tree and the bagging of 50 trees on the same problem, with 100 samples of 60 points.

import numpy as np
from sklearn.ensemble import BaggingRegressor
from sklearn.tree import DecisionTreeRegressor

rng = np.random.default_rng(3)
f = lambda x: np.sin(2 * np.pi * x)
sigma, n, M = 0.3, 60, 100
xg = np.linspace(0.02, 0.98, 49)[:, None]
models = {
    "one unbounded tree": lambda: DecisionTreeRegressor(random_state=0),
    "bagging of 50 trees": lambda: BaggingRegressor(DecisionTreeRegressor(), n_estimators=50, random_state=0),
}
for name, factory in models.items():
    P = np.empty((M, len(xg)))
    for m in range(M):
        x = rng.uniform(0, 1, (n, 1))
        y = f(x[:, 0]) + sigma * rng.standard_normal(n)
        P[m] = factory().fit(x, y).predict(xg)
    b2 = np.mean((P.mean(axis=0) - f(xg[:, 0])) ** 2)
    print(f"{name:24s} bias^2={b2:.4f} variance={np.mean(P.var(axis=0)):.4f}")

# Output:
# one unbounded tree       bias^2=0.0012 variance=0.0936
# bagging of 50 trees      bias^2=0.0007 variance=0.0430

The variance is halved (from 0.0936 to 0.0430) and the squared bias does not get worse (from 0.0012 to 0.0007). It is not divided by 50, because the 50 trees are trained on reshuffled samples of the same training set and so their errors are correlated; there is a reduction, but it has a limit. This is the idea that random forests push further by decorrelating the trees, the subject of chapter 06.

Beyond the classical picture

The "U" shape of the test error as a function of complexity is not the last word. Belkin and colleagues, in "Reconciling modern machine learning practice and the bias-variance trade-off" (arXiv:1812.11118), discuss the behaviour of models with more parameters than examples, in which the test error can fall again after the interpolation peak (double descent). This article does not reproduce the effect and does not report its results: the experiments above stay in the classical regime, with models that are small relative to the data.

What you can do now

You can write the decomposition of the squared error and check it with a simulation, tell a rigid model from an unstable one with learning curves, choose between more data, more regularisation and a less complex model according to what you see, and you know why averaging trees reduces the variance and not the bias. What remains is to understand how to estimate the test error when data are scarce: 05 - Cross-validation and model selection.

The same series also includes: 01 - Introduction to Machine Learning: Taxonomy and Workflow, 03 - Logistic Regression and Binary Classification and 12 - ML Deployment in Production: Serving and Monitoring.