Skip to main content

Nested cross-validation

Examples

Choosing a model and estimating its error are two different jobs. If we do both with the same data, the answer is optimistic by an amount we can predict.

On this page we measure that amount.

Setup

import numpy as np
import matplotlib.pyplot as plt

RNG = np.random.default_rng(5)
SIGMA = 0.15
DEGREES = np.arange(1, 11)


def f(x):
    return np.sin(2 * x) + 2 * (x - 0.5) ** 3 - 0.5 * x


def sample(n, rng):
    x = rng.uniform(0, 1, n)
    return x, f(x) + rng.normal(0, SIGMA, n)


def design(x, d):
    return np.vander(x, d + 1, increasing=True)


def fit(x, y, d):
    coef, *_ = np.linalg.lstsq(design(x, d), y, rcond=None)
    return coef


def mse(coef, x, y):
    return float(np.mean((design(x, len(coef) - 1) @ coef - y) ** 2))


def cross_val_mse(x, y, d, K, rng):
    folds = np.array_split(rng.permutation(len(x)), K)
    tot = cnt = 0
    for i in range(K):
        va = folds[i]
        tr = np.concatenate([folds[j] for j in range(K) if j != i])
        tot += mse(fit(x[tr], y[tr], d), x[va], y[va]) * len(va)
        cnt += len(va)
    return tot / cnt

Two procedures

Flat cross-validation cross-validates every degree and reports the smallest score.

Nested cross-validation splits off an outer test fold first. Inside the remaining data it cross-validates to pick a degree. It then refits that degree on the inner data and scores it on the outer fold, which took no part in the choice. This is repeated for every outer fold.

def flat_cv(x, y, K=10, rng=None):
    """Select and report with the same data. Optimistic."""
    rng = rng or np.random.default_rng()
    scores = [cross_val_mse(x, y, d, K, rng) for d in DEGREES]
    j = int(np.argmin(scores))
    return DEGREES[j], scores[j]


def nested_cv(x, y, K_outer=5, K_inner=10, rng=None):
    """Select inside, score outside."""
    rng = rng or np.random.default_rng()
    outer = np.array_split(rng.permutation(len(x)), K_outer)
    scores, chosen = [], []
    for i in range(K_outer):
        test = outer[i]
        rest = np.concatenate([outer[j] for j in range(K_outer) if j != i])
        inner = [cross_val_mse(x[rest], y[rest], d, K_inner, rng)
                 for d in DEGREES]
        d = DEGREES[int(np.argmin(inner))]
        scores.append(mse(fit(x[rest], y[rest], d), x[test], y[test]))
        chosen.append(d)
    return chosen, float(np.mean(scores))
1
This fold is put aside. Nothing after this line may look at it.
2
The whole selection runs on rest only.
3
We refit the chosen degree on all of rest, then score once on test.

How large is the gap

We repeat both procedures on fresh data sets and compare them to the truth. The truth here is the error of the selected model, measured on a large independent test set.

N = 80
REPS = 120

xt = RNG.uniform(0, 1, 20_000)
yt = f(xt) + RNG.normal(0, SIGMA, 20_000)

flat_scores, nested_scores, truth_scores = [], [], []

for _ in range(REPS):
    x, y = sample(N, RNG)
    d_flat, s_flat = flat_cv(x, y, K=10, rng=RNG)
    _, s_nested = nested_cv(x, y, K_outer=5, K_inner=10, rng=RNG)
    truth = mse(fit(x, y, d_flat), xt, yt)
    flat_scores.append(s_flat)
    nested_scores.append(s_nested)
    truth_scores.append(truth)

flat_scores = np.array(flat_scores)
nested_scores = np.array(nested_scores)
truth_scores = np.array(truth_scores)

print(f"{'':22s}{'mean':>9s}{'bias vs truth':>16s}")
print(f"{'flat CV (reported)':22s}{flat_scores.mean():9.4f}{flat_scores.mean() - truth_scores.mean():+16.4f}")
print(f"{'nested CV':22s}{nested_scores.mean():9.4f}{nested_scores.mean() - truth_scores.mean():+16.4f}")
print(f"{'true test error':22s}{truth_scores.mean():9.4f}{0.0:+16.4f}")
1
The benchmark. The degree that flat CV picked, refitted on all the training data, scored on 20 000 fresh points. 
                           mean   bias vs truth
flat CV (reported)       0.0232         -0.0014
nested CV                0.0260         +0.0014
true test error          0.0246         +0.0000
fig, ax = plt.subplots(figsize=(6.4, 3.0))
bins = np.linspace(min(flat_scores.min(), truth_scores.min()) * 0.9,
                   max(nested_scores.max(), truth_scores.max()) * 1.02, 34)
ax.hist(flat_scores, bins=bins, alpha=0.75, label="flat CV (reported)")
ax.hist(nested_scores, bins=bins, alpha=0.6, label="nested CV")
ax.axvline(truth_scores.mean(), color="#a8452f", ls="--", lw=1.4,
           label="mean true test error")
ax.set(xlabel="reported MSE", ylabel="count")
ax.set_title("What each procedure reports")
ax.legend()
plt.show()
Figure 27.1: What each procedure reports, over 120 repetitions. Flat CV sits to the left of the truth. Nested CV sits slightly to the right.

Nested cross-validation does not land exactly on the truth either. It sits a little above it. Each outer fold trains on only four fifths of the data, so the models it scores are slightly worse than the one we will finally fit on everything. Nested CV is mildly pessimistic and flat CV is optimistic. Of the two, being pessimistic is the safer mistake.

The gap is small here, because the degree of a polynomial is a search over ten options. It grows with the size of the search. A grid over three hyper-parameters with ten values each is a thousand comparisons, and the smallest of a thousand noisy numbers lies well below the truth.

The same experiment with a growing search. The blue bars are what flat cross-validation reports, the orange bars are the true test error of the model it picked. The gap grows with the number of models compared.

The same thing with scikit-learn

GridSearchCV is the inner loop. Putting it inside cross_val_score gives the outer loop.

from sklearn.linear_model import LinearRegression
from sklearn.model_selection import GridSearchCV, KFold, cross_val_score
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import PolynomialFeatures

x, y = sample(N, RNG)
X = x.reshape(-1, 1)

pipe = Pipeline([("poly", PolynomialFeatures()), ("lin", LinearRegression())])
grid = {"poly__degree": DEGREES}

inner = GridSearchCV(pipe, grid, cv=KFold(10, shuffle=True, random_state=0),
                     scoring="neg_mean_squared_error")

flat = -inner.fit(X, y).best_score_
outer = -cross_val_score(inner, X, y, cv=KFold(5, shuffle=True, random_state=1),
                         scoring="neg_mean_squared_error").mean()

print(f"flat   (best_score_):  {flat:.4f}")
print(f"nested (cross_val):    {outer:.4f}")
1
best_score_ is the score of the winner. It is not a test error.
2
Here the whole grid search is refitted inside every outer fold. 
flat   (best_score_):  0.0225
nested (cross_val):    0.0208

The rule

Any data we used to make a choice cannot also be used to estimate how good that choice is.

Choices include the model family, the hyper-parameters, the features, which points to drop as outliers, and which of three attempts to write up.

Standardization is also a choice. Fit the scaler inside the fold, on the training part only. A scaler fitted on all the data has seen the mean of the validation fold.

For the project, cross-validate to choose, and report the score from data that was held out before you started choosing. If you have tuned against your test set even once, it is a validation set now, and you need a new one.