Bias-Variance Tradeoff, in depth: deriving the decomposition, measuring it by simulation, and using learning curves to decide what to fix

The quantity being decomposed

Assume the data come from y = f(x) + e, where f is the unknown true function and e is noise with mean zero and variance sigma^2. You draw a training set D, fit a model, and get a predictor f_D. The training set is random, so f_D is random as well. Draw a different sample and you get a different fit.

Bias-variance analysis is about the expected squared error at a point x, averaged over both the noise in a fresh test label and the choice of training set: E_D,e[(y - f_D(x))^2]. This is not the error of the one model you trained. It is the average error of your whole procedure, meaning the model class, the training algorithm, the hyperparameters and the sample size. That is why the decomposition talks about recipes, not individual models.

Deriving the three terms

Write m(x) = E_D[f_D(x)] for the average prediction over many training sets. Then add and subtract both f(x) and m(x):

y - f_D(x) = (y - f(x)) + (f(x) - m(x)) + (m(x) - f_D(x))
             \_ noise _/   \___ bias ___/   \_ deviation _/

Square this and take expectations. There are three cross terms, and all of them vanish. The noise y - f(x) has mean zero and is independent of the training set. The deviation m(x) - f_D(x) has mean zero over training sets by the definition of m. And f(x) - m(x) is a constant at a fixed x. What remains is:

E[(y - f_D(x))^2] = (f(x) - m(x))^2      # bias^2: systematic error of the recipe
                  + E_D[(f_D(x) - m(x))^2]  # variance: wobble across training sets
                  + sigma^2                 # irreducible noise

Average over the test distribution of x to get the expected test MSE. All three terms are non-negative. Noise sets a floor that no model can go below. If your validation MSE is already close to sigma^2, more modelling work is wasted.

This clean additive split is specific to squared loss. For 0-1 classification loss there are several competing definitions of bias and variance, such as those of Kong and Dietterich, Breiman, and Domingos. In those versions the terms can interact, and extra variance can even reduce error at points where the average prediction is wrong. Cross-entropy has its own decomposition built on a different notion of the average prediction. Use the squared-loss version for intuition. Do not carry its exact arithmetic over to classifiers.

Measuring the terms by simulation

You cannot measure bias and variance on a real problem, because you do not know f and cannot draw many independent training sets. In a simulation you can. The script below fixes f(x) = sin(2 pi x), noise standard deviation 0.3 (so sigma^2 = 0.09), and 30 training points per set. It then fits polynomials of several degrees to 500 independent training sets.

import numpy as np

rng = np.random.default_rng(0)
f = lambda x: np.sin(2 * np.pi * x)
SIGMA, N_TRAIN, N_SETS = 0.3, 30, 500
x_test = np.linspace(0.05, 0.95, 200)

def fit_predict(deg, x, y, x_new, lam=0.0):
    X = np.vander(x, deg + 1, increasing=True)
    Xn = np.vander(x_new, deg + 1, increasing=True)
    if lam == 0:
        return Xn @ np.linalg.lstsq(X, y, rcond=None)[0]
    A = X.T @ X + lam * np.eye(deg + 1)
    w = np.linalg.solve(A, X.T @ y)
    return Xn @ w

def decompose(deg, lam=0.0, n=N_TRAIN):
    preds = np.empty((N_SETS, x_test.size))
    for s in range(N_SETS):
        x = rng.uniform(0, 1, n)
        y = f(x) + rng.normal(0, SIGMA, n)
        preds[s] = fit_predict(deg, x, y, x_test, lam)
    mean_pred = preds.mean(axis=0)
    bias2 = np.mean((mean_pred - f(x_test)) ** 2)
    var = np.mean(preds.var(axis=0))
    y_test = f(x_test) + rng.normal(0, SIGMA, (N_SETS, x_test.size))
    mse = np.mean((preds - y_test) ** 2)
    return bias2, var, mse

print("deg  bias2   var     noise  sum     measured")
for deg in (1, 3, 5, 9, 12):
    b, v, m = decompose(deg)
    print(f"{deg:>3}  {b:.4f}  {v:.4f}  {SIGMA**2:.3f}  {b+v+SIGMA**2:.4f}  {m:.4f}")
print("ridge deg 12")
for lam in (1e-6, 1e-4, 1e-2, 1e-1, 1.0):
    b, v, m = decompose(12, lam)
    print(f"{lam:g}  {b:.4f}  {v:.4f}  {m:.4f}")
print("deg 9, more data")
for n in (30, 100, 300):
    b, v, m = decompose(9, n=n)
    print(f"{n}  {b:.4f}  {v:.4f}  {m:.4f}")

The script calls decompose for degrees 1, 3, 5, 9 and 12, then for ridge-penalised degree 12, then for degree 9 at three sample sizes, all from one seeded generator, so call order matters. The measured column is a direct Monte Carlo estimate of test MSE. It is a check on the arithmetic, and it should agree with bias squared plus variance plus noise.

Worked example: reading the numbers

A straight line (degree 1) cannot follow a sine wave, so bias dominates and variance is small. Degree 3 removes almost all of the bias, and total error is close to the 0.09 floor. Degree 5 has even less bias, but its variance has doubled and the total is slightly worse. At degree 9 the variance is about 0.45, five times the noise. At degree 12, with 13 coefficients and only 30 points, some samples produce wild fits, and the averaged variance passes a thousand. In every row, sum and measurement agree to within half a percent, which is the decomposition at work.

Note what happens to bias at degree 12. The model class can represent the target, yet the measured bias is large. A few extreme fits drag the average prediction away from the truth, so high variance can leak into the bias estimate when the number of simulated sets is finite. The degree 9 rows show a related effect. A later call with the same settings in the same run printed a variance of 1.5476 instead of 0.4471. High-capacity fits have heavy-tailed errors, so variance estimates for them are noisy too. Expect that from real models, as well.

The knobs, and which term each one moves

Regularisation trades a little bias for a lot of variance. Using the same script with degree 12 and a ridge penalty lam gave:

Even a tiny penalty takes the degree-12 model from catastrophic to reasonable. Error is lowest around 1e-4. Beyond that point bias rises faster than variance falls. (This penalty also shrinks the intercept, which a production ridge usually leaves unpenalised.)

More data reduces variance without adding bias. For degree 9 the same run gave a variance of 1.5476 at 30 points, 0.0078 at 100 points and 0.0025 at 300 points. Measured MSE fell from 1.6451 to 0.0975 and then 0.0923, which is close to the floor. A model that is too complex for 30 points is fine with 300. Capacity is only "too high" relative to the amount of data.

Ensembling moves the terms in different ways. Bagging averages models trained on bootstrap resamples, which cuts variance when the members' errors are not perfectly correlated, and leaves bias roughly unchanged. That is why random forests use deep, low-bias trees. Boosting fits members in sequence to the remaining errors, which mainly reduces bias. Early stopping in gradient-trained networks behaves much like regularisation. Better features can reduce bias without adding capacity, by moving the target closer to something the model class represents easily.

Diagnosing a real model with learning curves

In practice you have one dataset, not 500. The usable signal is how training and validation error change with training set size. scikit-learn's learning_curve computes this with cross-validation:

import numpy as np
from sklearn.model_selection import learning_curve, KFold

sizes, train_scores, val_scores = learning_curve(
    model, X, y,
    train_sizes=np.linspace(0.1, 1.0, 8),
    cv=KFold(n_splits=5, shuffle=True, random_state=0),
    scoring="neg_mean_squared_error",
)
train_mse = -train_scores.mean(axis=1)
val_mse, val_spread = -val_scores.mean(axis=1), val_scores.std(axis=1)
for n, tr, va, sd in zip(sizes, train_mse, val_mse, val_spread):
    print(f"n={n:6d}  train={tr:.4f}  val={va:.4f} +/- {sd:.4f}  gap={va - tr:.4f}")

Read the result like this. If both curves converge early to a level well above your noise estimate, with a small gap, the model is bias-limited. More data will not help. Add capacity or features, or reduce regularisation. If training error is low and there is a large gap that is still shrinking at full size, the model is variance-limited. Get more data, regularise, or ensemble. A large standard deviation across folds is a direct sign of variance as well. Splits must respect groups and time, or leakage will make the variance look smaller than it is. The train, validation and test split guide covers this, and the Spark cross-validation deep dive shows fold design at scale.

Here is an illustrative case, with made-up numbers. A gradient-boosted demand model shows training MSE of 0.8 and validation MSE of 2.9 at full size, the gap is still shrinking, and fold spread is wide. This model is variance-limited. The first things to try are a lower learning rate with more rounds and early stopping, shallower trees and row subsampling, in each case compared on the same folds. If the gap closes and validation error settles near the training error but well above the noise estimate, the problem has become bias. Then the next step is features, such as holiday calendars or lagged demand, not more regularisation.

To estimate the noise floor, look for repeated measurements or duplicate inputs with different labels, or measure inter-annotator disagreement. Without a floor estimate you cannot tell "high bias" apart from "the labels are noisy".

Where the classical picture breaks

The U-curve assumes that capacity keeps adding variance. Heavily over-parameterised models, such as wide neural networks or minimum-norm interpolating regressors, often do not follow that pattern. As parameters grow past the point where training error reaches zero, test error can peak and then fall again. This is double descent. The implicit bias of the optimiser, such as gradient descent's tendency to find low-norm solutions, acts as regularisation that does not appear in the parameter count. The decomposition still holds as an identity. What fails is the assumption that parameter count is a good measure of capacity. For deep models, judge capacity by validation behaviour, not by the parameter count.

Failure modes and trade-offs

• Leakage looks like low variance. If near-duplicates cross the split, validation error is optimistic and the gap is hidden. Deduplicate before you diagnose.
• Tuning on validation data adds variance you cannot see. A large hyperparameter search overfits the validation set. Keep a test set you look at once.
• Distribution shift is neither bias nor variance. The decomposition assumes test data come from the training distribution. A model that degrades in production may be facing a new f, not overfitting the old one. Monitor for drift separately, as the evaluation frameworks guide describes.
• Picking a point on the curve is a business decision. A slightly higher-bias model can be more stable from retrain to retrain, which matters when predictions drive decisions people notice. Lower average error is not the only goal.

What to do next

1. Run the simulation above and change one thing at a time: noise level, sample size, degree. Predict each result before you run it.
2. Estimate the noise floor for your own task from duplicate labels or repeated measurements.
3. Plot learning curves for your current model with group-aware or time-aware splits, and decide whether it is bias-limited or variance-limited.
4. If it is variance-limited, try regularisation strength, more data and bagging, in that order, and compare against the same folds.
5. If it is bias-limited, add features or capacity, and stop collecting more of the same data.
6. Record the fold-to-fold spread alongside the mean in every model comparison.