Overfitting & regularisation

A model that does well on the data it was trained on has proved nothing. This section is about the gap between that score and how it performs on anything else, what causes the gap to open, and how to see it happening in a learning curve rather than discovering it after launch. The fixes are a short list and they have an order, which is most of what makes this practical rather than theoretical.

The symptom#

The model performs far better on training data than on data it has not seen. It has memorised the examples, including their noise, rather than learning the pattern.

train accuracy 0.99    test accuracy 0.71     overfitting
train accuracy 0.72    test accuracy 0.70     underfitting
train accuracy 0.91    test accuracy 0.88     about right

That table is the summary. What is worth seeing once, properly, is the gap opening while you watch, because it does not open in the way most descriptions imply. Here is a gradient boosting model fitted to 700 rows of a problem where 15 per cent of the labels have been flipped, so there is a real pattern and real noise sitting on top of it. Every one of the 500 trees is scored on the training rows and on 3,000 rows the model never sees.

import numpy as np
from sklearn.datasets import make_classification
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.metrics import accuracy_score, log_loss

X, y = make_classification(n_samples=10_000, n_features=25, n_informative=6,
                           n_redundant=0, flip_y=0.15, random_state=0)
Xva, yva = X[7000:], y[7000:]          # 3000 rows, fixed, never trained on
Xtr, ytr = X[:700], y[:700]            # the training set we start with

m = GradientBoostingClassifier(n_estimators=500, max_depth=4, random_state=0).fit(Xtr, ytr)
tr = list(m.staged_predict_proba(Xtr))
va = list(m.staged_predict_proba(Xva))

print("trees  train acc  val acc    gap   val log loss")
for i in (9, 24, 49, 99, 199, 499):
    a = accuracy_score(ytr, tr[i].argmax(1))
    b = accuracy_score(yva, va[i].argmax(1))
    print(f"{i+1:5d}  {a:9.3f}  {b:7.3f}  {a-b:+.3f}   {log_loss(yva, va[i]):8.3f}")
best = int(np.argmin([log_loss(yva, p) for p in va]))
print("best validation log loss at tree", best + 1)

# trees  train acc  val acc    gap   val log loss
#    10      0.873    0.774  +0.099      0.521
#    25      0.913    0.787  +0.126      0.478
#    50      0.961    0.785  +0.176      0.476
#   100      0.999    0.783  +0.216      0.488
#   200      1.000    0.781  +0.219      0.528
#   500      1.000    0.788  +0.212      0.717
# best validation log loss at tree 37

Read the last column before the others. Validation accuracy is almost flat across the whole run, so a team tracking accuracy alone would see nothing wrong and ship the 500-tree model. Validation log loss, which scores how confident the model was rather than only whether it was right, bottoms out at tree 37 and then rises to 0.717 — half again as bad as where it started. That is what overfitting does first. The model does not begin getting more answers wrong; it begins being more certain about the ones it gets wrong, and that is the part that shows up as a badly calibrated probability in production long before it shows up as a drop in accuracy.

The other number to take from the table is 0.212. Train minus validation is the gap, and it is the quantity you monitor. A gap that grows while validation stands still is the signature; a gap near zero with both numbers poor is the opposite problem and needs the opposite treatment.

Bias and variance#

These two words get defined and then never used. They are worth making concrete, because the whole of the rest of this section is a decision about which of the two you have. Take a known function, sine over the unit interval, add noise of a known size, and fit polynomials of increasing degree. Because the true function is known, the error of a fitted model can be split into exactly three parts: how far the average fit sits from the truth (bias), how much the fit moves when the training sample changes (variance), and the noise, which nothing can remove.

rng = np.random.default_rng(0)
NOISE, GRID = 0.3, np.linspace(0, 1, 100)
truth = np.sin(2 * np.pi * GRID)

print("degree    bias^2  variance   noise     sum   test MSE")
for deg in (0, 1, 3, 5, 7, 9):
    preds = []
    for _ in range(300):                      # 300 independent training sets of 30 points
        x = rng.random(30)
        yy = np.sin(2 * np.pi * x) + rng.normal(0, NOISE, 30)
        preds.append(np.polyval(np.polyfit(x, yy, deg), GRID))
    P = np.array(preds)
    bias2 = float(((P.mean(0) - truth) ** 2).mean())
    var = float(P.var(0).mean())
    test = truth + rng.normal(0, NOISE, (300, 100))   # fresh noisy labels, as a check
    mse = float(((P - test) ** 2).mean())
    print(f"{deg:6d} {bias2:9.3f} {var:9.3f} {NOISE**2:7.3f} {bias2+var+NOISE**2:7.3f} {mse:10.3f}")

# degree    bias^2  variance   noise     sum   test MSE
#      0     0.495     0.023   0.090   0.608      0.607
#      1     0.204     0.021   0.090   0.315      0.314
#      3     0.005     0.017   0.090   0.112      0.113
#      5     0.000     0.074   0.090   0.165      0.164
#      7     0.004     0.573   0.090   0.667      0.669
#      9     0.022     6.499   0.090   6.611      6.607

Bias falls to nothing by degree 3 and stays there — a cubic can already trace a sine closely enough. Variance is 0.017 at degree 3 and 6.499 at degree 9, a factor of nearly four hundred, from thirty training points that differ only in which thirty they are. The best total error is at degree 3, and it is 0.112 against a noise floor of 0.090, so 80 per cent of the remaining error is irreducible. Nothing in this section can beat 0.090, and a project that has not worked out its own equivalent of that number will keep spending on the part that cannot move.

• High bias: the model is too simple to capture the pattern. Both scores are poor. Underfitting.
• High variance: the model is complex enough to fit noise. Training score far exceeds test. Overfitting.
• More capacity reduces bias and increases variance. Finding the balance for your data volume is the bias-variance tradeoff.
• Above degree 5 here, bias is flat and total error is rising: every unit of extra capacity is buying variance and nothing else. That is the exact condition under which reducing capacity helps and nothing else will.

Learning curves#

Plot training and validation performance against training set size, or against epochs. The shape tells you which problem you have and what will help.

Curves converged, both poor          more data will not help
                                      you need a better model

Large gap, validation flat           more data will help
                                      or reduce complexity

Validation rising then falling       stop earlier

That is the rule. Here are both cases measured, so you can see what each actually looks like before you have to read one under pressure. Case A is a random forest on a thirty-column problem; case B is logistic regression on two interleaved crescents, a shape no straight line can separate.

from sklearn.datasets import make_moons
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import learning_curve

Xa, ya = make_classification(n_samples=20000, n_features=30, n_informative=20,
                             n_redundant=0, flip_y=0.05, class_sep=0.8, random_state=0)
Xb, yb = make_moons(n_samples=20000, noise=0.3, random_state=0)

for name, model, XX, yy in [
    ("A  random forest, 30 columns", RandomForestClassifier(n_estimators=100, random_state=0), Xa, ya),
    ("B  logistic regression, two moons", LogisticRegression(), Xb, yb)]:
    n, t, v = learning_curve(model, XX, yy, train_sizes=[200, 500, 1500, 5000, 15000],
                             cv=4, scoring="accuracy", shuffle=True, random_state=0)
    print(name)
    for i, s in enumerate(n):
        print(f"   n={s:6d}   train {t[i].mean():.3f}   val {v[i].mean():.3f}   gap {t[i].mean()-v[i].mean():+.3f}")

# A  random forest, 30 columns
#    n=   200   train 1.000   val 0.697   gap +0.303
#    n=   500   train 1.000   val 0.796   gap +0.204
#    n=  1500   train 1.000   val 0.863   gap +0.137
#    n=  5000   train 1.000   val 0.910   gap +0.090
#    n= 15000   train 1.000   val 0.934   gap +0.066
# B  logistic regression, two moons
#    n=   200   train 0.867   val 0.845   gap +0.022
#    n=   500   train 0.857   val 0.848   gap +0.009
#    n=  1500   train 0.851   val 0.850   gap +0.001
#    n=  5000   train 0.851   val 0.849   gap +0.002
#    n= 15000   train 0.850   val 0.849   gap +0.001

Case A: the gap starts at 0.303 and is still closing at 15,000 rows, and validation has climbed from 0.697 to 0.934 without flattening. Every additional row is still buying accuracy, so the answer is more data, and the curve even tells you roughly how much — the last tripling bought 2.4 points, so the next one will buy less than that and you can decide whether it is worth the labelling cost before spending it.

Case B: the gap is 0.022 at two hundred rows and 0.001 at fifteen thousand. The model is not memorising anything; it has extracted everything it is capable of extracting and both curves have converged at 0.850. Seventy-five times more data moved validation accuracy by four thousandths. More data is not the fix and neither is regularisation, which only ever removes capacity from a model that does not have enough. A random forest on the same rows reaches 0.904, so the fix here is a model that can bend.

The levers#

1. 1.More data. Always the best fix when available, and usually the least available.
2. 2.Early stopping. Free, and stops training at the point validation was best.
3. 3.Regularisation. L2 penalises large weights; L1 drives some to zero and selects features.
4. 4.Dropout, in neural networks. Randomly ignore units during training so the model cannot rely on any one.
5. 5.Reduce capacity. Fewer parameters, shallower trees, fewer features.
6. 6.Data augmentation, where the domain allows it.

The order is not arbitrary: it runs from the lever that raises the ceiling to the levers that only lower the gap. Applied to the overfitting model from the top of this section, four of them do something quite different from what the list implies.

import time

def score(model, Xt=Xtr, yt=ytr):
    t0 = time.perf_counter()
    model.fit(Xt, yt)
    return (accuracy_score(yt, model.predict(Xt)),
            accuracy_score(yva, model.predict(Xva)), time.perf_counter() - t0)

print("lever                          train    val     gap   fit s")
for label, model in [
    ("nothing (500 trees)", GradientBoostingClassifier(n_estimators=500, max_depth=4, random_state=0)),
    ("early stopping (37 trees)", GradientBoostingClassifier(n_estimators=37, max_depth=4, random_state=0)),
    ("less capacity (depth 2)", GradientBoostingClassifier(n_estimators=500, max_depth=2, random_state=0)),
    ("slower, subsampled", GradientBoostingClassifier(n_estimators=500, max_depth=4,
                                                      subsample=0.5, learning_rate=0.02, random_state=0)),
]:
    a, b, s = score(model)
    print(f"{label:30s} {a:.3f}  {b:.3f}  {a-b:+.3f}  {s:5.1f}")
a, b, s = score(GradientBoostingClassifier(n_estimators=500, max_depth=4, random_state=0), X[:7000], y[:7000])
print(f"{'10x the data (7000 rows)':30s} {a:.3f}  {b:.3f}  {a-b:+.3f}  {s:5.1f}")

# lever                          train    val     gap   fit s
# nothing (500 trees)            1.000  0.788  +0.212    1.8
# early stopping (37 trees)      0.946  0.786  +0.159    0.1
# less capacity (depth 2)        0.999  0.765  +0.234    0.9
# slower, subsampled             1.000  0.790  +0.210    0.9
# 10x the data (7000 rows)       0.974  0.842  +0.132   22.2

Only one lever moved validation accuracy, and it is the one at the top of the list: ten times the data bought 5.4 points, from 0.788 to 0.842, for thirteen times the fit time. Early stopping gave back 0.002 of accuracy and took eighteen times less time to train — and, as the first table showed, it halved the log loss, which is the reason to do it. Slowing the learning rate and subsampling gained 0.002, inside the noise. Reducing depth to 2 lost 2.3 points: the model was not too complex for the problem, it was too complex for 700 rows, and taking capacity away made it worse at both.

That last row is the one worth remembering, because it is the standard mistake. Seeing a gap of 0.212 and reaching for less capacity treats the gap as the thing to minimise. The gap is a symptom. What you are minimising is validation error, and a change that closes the gap by making the training score worse has usually done nothing for you. Check both numbers after every lever, and keep the one that moved the second.

Two questions decide whether a model generalises: what stops it memorising the training data, and whether the validation score you are reading is honest. This section covers both. First the regularisers, an L2 or L1 penalty, early stopping, dropout and augmentation, and then the split, which is where most inflated scores come from: the same customer on both sides, a scaler fitted before the split, or time-ordered rows that were shuffled. A validation number produced by a careless split is worse than no number, because you act on it.

Penalise large weights#

An L2 penalty adds the sum of squared weights to the loss, so the model only uses large weights when they clearly pay for themselves. L1 adds the sum of absolute weights, which pushes many weights to exactly zero and so also selects features.

from sklearn.datasets import make_classification
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score

# 1,400 rows, 200 columns, 15 of which carry any signal at all.
X, y = make_classification(n_samples=1400, n_features=200, n_informative=15,
                           n_redundant=0, flip_y=0.1, random_state=0)

for C in [0.01, 0.1, 1, 10, 100]:          # C is the inverse of penalty strength
    model = LogisticRegression(C=C, max_iter=2000)
    scores = cross_val_score(model, X, y, cv=5, scoring="roc_auc")
    print(f"C={C:<6} AUC {scores.mean():.3f} +/- {scores.std():.3f}")

# C=0.01   AUC 0.764 +/- 0.019
# C=0.1    AUC 0.740 +/- 0.018
# C=1      AUC 0.734 +/- 0.019
# C=10     AUC 0.733 +/- 0.019
# C=100    AUC 0.733 +/- 0.019

Extending it downwards is worth doing properly, because it reveals the thing the two penalties do differently. Below is the same data split into 400 training rows and a thousand held out, swept from a penalty so strong the model can barely move to one so weak it is effectively absent, for both L2 and L1, with the number of columns the fitted model actually uses.

import numpy as np
from sklearn.metrics import roc_auc_score
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

Xtr, ytr, Xva, yva = X[:400], y[:400], X[400:], y[400:]

def fit(**kw):
    p = make_pipeline(StandardScaler(), LogisticRegression(max_iter=5000, **kw)).fit(Xtr, ytr)
    return (roc_auc_score(ytr, p.predict_proba(Xtr)[:, 1]),
            roc_auc_score(yva, p.predict_proba(Xva)[:, 1]),
            int((np.abs(p[-1].coef_) > 1e-6).sum()))

print("penalty      C   train AUC  val AUC     gap   columns used")
for C in (1e-5, 1e-3, 1e-2, 1e-1, 1, 10):
    t, v, nz = fit(C=C)
    print(f"L2   {C:8} {t:9.3f} {v:8.3f} {t-v:+8.3f} {nz:9d}")
for C in (1e-2, 1e-1, 1, 10):
    t, v, nz = fit(C=C, penalty="l1", solver="liblinear")
    print(f"L1   {C:8} {t:9.3f} {v:8.3f} {t-v:+8.3f} {nz:9d}")

# penalty      C   train AUC  val AUC     gap   columns used
# L2      1e-05     0.900    0.648   +0.252       200
# L2      0.001     0.911    0.647   +0.264       200
# L2       0.01     0.941    0.644   +0.297       200
# L2        0.1     0.974    0.637   +0.337       200
# L2          1     0.998    0.628   +0.370       200
# L2         10     1.000    0.623   +0.377       200
# L1       0.01     0.500    0.500   +0.000         0
# L1        0.1     0.876    0.728   +0.148        54
# L1          1     0.991    0.646   +0.345       167
# L1         10     1.000    0.632   +0.368       179

L2 flattens out at 0.648 and cannot be pushed past it. Shrinking every weight towards zero still leaves 185 noise columns contributing a small amount each, and a small amount each, across 185 columns, is most of the model. L1 at the right strength deletes them outright and reaches 0.728. The top row of the L1 block is the other half of the lesson: at C=0.01 the penalty deletes everything, the model has no columns left and scores exactly 0.500, which is a coin flip. Both failure modes are a badly chosen strength rather than a badly chosen penalty, which is why the strength is swept rather than picked.

The rule that follows is short. When you have more columns than rows and expect most of them to be junk, L1 is doing a different job from L2 and usually the job you want. When the columns are all plausibly relevant and correlated with each other, L2 is better, because L1 picks one of a correlated pair arbitrarily and drops the other, which makes the fitted model unstable across refits even when its accuracy is fine.

Stop before it memorises#

For models trained in steps, such as neural networks and gradient boosting, validation loss usually falls, flattens, then rises as the model starts memorising. Early stopping keeps the weights from the best validation point. It is the cheapest regulariser there is.

from sklearn.ensemble import GradientBoostingClassifier
from sklearn.metrics import log_loss

model = GradientBoostingClassifier(n_estimators=1, max_depth=3,
                                   warm_start=True, random_state=0)
best, patience, bad, kept = float("inf"), 5, 0, 0
for epoch in range(1, 201):
    model.set_params(n_estimators=epoch)
    model.fit(X[:400], y[:400])                 # warm_start adds one more tree
    val = log_loss(y[400:], model.predict_proba(X[400:]))
    if val < best:
        best, bad, kept = val, 0, epoch         # save(model) belongs here
    else:
        bad += 1
        if bad >= patience:
            break
print(f"stopped at epoch {epoch}, keeping epoch {kept}, val log loss {best:.4f}")

# stopped at epoch 22, keeping epoch 17, val log loss 0.6071

Dropout and augmentation#

Dropout randomly zeroes activations during training, so the network cannot rely on any single path. Augmentation creates new training examples by transforming existing ones: flips and crops for images, noise for audio. Augmentation is often the strongest regulariser for images, because it adds real information about what should not change the label.

Cross-validation, and its traps#

K-fold cross-validation trains on k minus 1 parts and validates on the remaining one, rotating, so every example is used for validation once. It gives a mean and a spread. The spread matters: if two settings differ by less than it, you cannot tell them apart.

That sentence is easy to nod at and hard to feel, so here it is as a result. Two forest depths are compared on the same 600 rows, eight times over. Nothing changes between the eight runs except the random seed that decides where the fold boundaries fall.

from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import KFold

Xs, ys = make_classification(n_samples=600, n_features=20, n_informative=8,
                             flip_y=0.15, random_state=0)
print("seed   depth 6   depth 10   winner")
for seed in range(8):
    cv = KFold(5, shuffle=True, random_state=seed)
    a = cross_val_score(RandomForestClassifier(max_depth=6, n_estimators=200, random_state=0),
                        Xs, ys, cv=cv, scoring="roc_auc").mean()
    b = cross_val_score(RandomForestClassifier(max_depth=10, n_estimators=200, random_state=0),
                        Xs, ys, cv=cv, scoring="roc_auc").mean()
    print(f"{seed:4d}  {a:8.4f}  {b:9.4f}   depth {6 if a > b else 10}")

# seed   depth 6   depth 10   winner
#    0    0.9218     0.9194   depth 6
#    1    0.9238     0.9220   depth 6
#    2    0.9242     0.9250   depth 10
#    3    0.9251     0.9219   depth 6
#    4    0.9271     0.9279   depth 10
#    5    0.9251     0.9235   depth 6
#    6    0.9267     0.9260   depth 6
#    7    0.9242     0.9253   depth 10

Depth 6 wins five times and depth 10 wins three, and the margins are in the third decimal place. A tuning run that reported the seed-0 result would say depth 6 by 0.0024 and a run that happened to use seed 4 would say depth 10 by 0.0008, and both would be presented as a finding. The honest summary of this experiment is that the two settings are indistinguishable on 600 rows. If you need to tell them apart, repeat the cross-validation with several seeds and compare the spread of the means to the gap between them, or get more data — those are the only two moves, and picking the winner of one split is not one of them.

• Classes are imbalanced: use stratified folds, so each fold keeps the class ratio.
• Several rows belong to one patient, user or document: use group folds, so one group never appears on both sides. Otherwise the model is tested on people it has already seen.
• Rows are ordered in time: train on the past, validate on the future, never shuffle.
• Any preprocessing that learns from data, such as scaling or target encoding, must be fitted inside each fold. A sklearn Pipeline does this for you.

import numpy as np
from sklearn.model_selection import GroupKFold, TimeSeriesSplit
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler

# patient_ids is the id column of your own data, one entry per row.
patient_ids = np.repeat(np.arange(280), 5)      # 280 patients, 5 visits each

pipe = make_pipeline(StandardScaler(), LogisticRegression(max_iter=2000))
cross_val_score(pipe, X, y, cv=GroupKFold(n_splits=5), groups=patient_ids)
cross_val_score(pipe, X, y, cv=TimeSeriesSplit(n_splits=5))

The second item on that list is the one that costs the most and is spotted the least, so it is worth putting a number on. Two hundred patients, eight visits each, sixteen hundred rows. Every visit is drawn from its patient's underlying profile plus a little noise, and the label belongs to the patient rather than to the visit — which is how patient data, user data and document data are actually shaped.

rng = np.random.default_rng(0)
n_patients, per_patient, d = 200, 8, 12
profile = rng.normal(size=(n_patients, d))              # one latent profile per patient
label = (profile[:, 0] + rng.normal(0, 1.0, n_patients) > 0).astype(int)
groups = np.repeat(np.arange(n_patients), per_patient)
Xg = profile[groups] + rng.normal(0, 0.4, (n_patients * per_patient, d))
yg = label[groups]                                      # every visit shares its patient's label

m = RandomForestClassifier(n_estimators=200, random_state=0)
a = cross_val_score(m, Xg, yg, cv=KFold(5, shuffle=True, random_state=0), scoring="roc_auc")
b = cross_val_score(m, Xg, yg, cv=GroupKFold(5), groups=groups, scoring="roc_auc")
print(f"shuffled k-fold   AUC {a.mean():.3f} +/- {a.std():.3f}")
print(f"group k-fold      AUC {b.mean():.3f} +/- {b.std():.3f}")

# shuffled k-fold   AUC 0.947 +/- 0.016
# group k-fold      AUC 0.735 +/- 0.044

Twenty-one points of AUC, and the wrong number is the one that looks better and has the tighter spread. With shuffled folds the model does not need to generalise to a new patient; for almost every validation row there are seven near-copies of it in training, so it is being asked to recognise people it has already met. Group folds ask the question production will ask, which is whether the model works on a patient it has never seen, and the answer is 0.735.

Notice the spread as well as the mean. Group k-fold reports plus or minus 0.044 against 0.016, nearly three times wider, because there are only 200 independent units in this dataset rather than 1,600, and the honest error bar reflects that. An inflated score usually arrives with an inflated confidence, which is why a suspiciously tight cross-validation spread on entity-shaped data is worth checking before the mean is.

Data leakage#

Every trap above is one species of the same animal. Data leakage is information reaching the model that would not be available at prediction time, so the offline score measures something the production system will never get to do. It is the most common reason a model that validated at 0.95 serves at 0.70, and it is worth seeing rather than being warned about, because the warning never feels urgent enough.

import numpy as np
from sklearn.feature_selection import SelectKBest, f_classif
from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import cross_val_score
from sklearn.pipeline import make_pipeline

rng = np.random.default_rng(0)
X = rng.normal(size=(400, 800))        # 800 columns of pure noise
y = rng.integers(0, 2, 400)            # a coin flip: nothing is predictable

# WRONG: pick the columns that correlate with y using all the data,
# then cross-validate on the survivors. The selection has already seen
# every validation fold, so the folds are no longer held out.
corr = np.array([abs(np.corrcoef(X[:, j], y)[0, 1]) for j in range(800)])
best = np.argsort(-corr)[:20]
print(round(cross_val_score(LogisticRegression(), X[:, best], y, cv=5).mean(), 3))

# RIGHT: the selection is a step in the pipeline, so it is refitted on
# the training part of each fold and never sees the validation part.
pipe = make_pipeline(SelectKBest(f_classif, k=20), LogisticRegression())
print(round(cross_val_score(pipe, X, y, cv=5).mean(), 3))

# 0.685   <- on data with no signal whatsoever
# 0.503   <- the truth