The Mathematics Behind Model Validation

The mathematics under the scores: why training error is optimistic by a knowable amount, what cross-validation actually estimates, why its reported standard error is wrong, and how tuning inflates the winner.

Every number a model reports about itself is an estimate. The accuracy printed at the end of a script is not a property of the model; it is a random variable, computed from one particular sample, that happens to land near a quantity nobody can observe. Treating it as a fact rather than an estimate is the most common error in applied machine learning, and it is the reason models that scored 94 percent in a notebook return 81 percent in production.

Model validation is therefore a branch of estimation theory wearing engineering clothes. The questions it asks are the questions any statistician asks of any estimator. Is it biased, and in which direction? What is its variance? How does that variance shrink with more data? And what happens to all of the above when you use the same estimate a hundred times to choose between a hundred models? This article answers those four questions, and the answers explain why cross-validation exists, why the standard error people report for it is wrong, and why a tuned model’s validation score is a number you should not believe.

What we are trying to estimate

Assume the data comes from some fixed but unknown distribution P(x,y)P(x, y). The quantity we actually care about is the expected loss of our model on a fresh draw from that distribution, the generalisation error, also called the risk:

R(f^)=𝔼(x,y)P[L(y,f^(x))] R(\hat{f}) = \mathbb{E}_{(x, y) \sim P}\left[L\left(y, \hat{f}(x)\right)\right]

This is the number that decides whether the model is useful. It is also permanently unobservable, since computing it would require the distribution we do not have. Everything validation does is construct estimators of RR from the finite sample we possess.

The obvious candidate is the average loss on the training data, the empirical risk:

R^train(f^)=1ni=1nL(yi,f^(xi))\hat{R}{\text{train}}(\hat{f}) = \frac{1}{n}\sum{i=1}^{n} L\left(y_i, \hat{f}(x_i)\right)

and it is a catastrophically bad estimator, for a reason worth making precise rather than waving at.

Why training error lies, quantified

The model was chosen by minimising the loss on these specific points, so it has partly fitted their noise rather than the signal. Efron’s covariance formula makes the size of the resulting optimism exact:

𝔼[R]𝔼[R^train]=2ni=1nCov(y^i,;yi)\mathbb{E}\left[R\right] – \mathbb{E}\left[\hat{R}{\text{train}}\right] = \frac{2}{n}\sum{i=1}^{n} \operatorname{Cov}\left(\hat{y}_i, ; y_i\right)

Read the right-hand side carefully, because it is one of the most illuminating expressions in statistics. The optimism of the training score is exactly the degree to which each fitted value tracks its own observed target. A model that ignores yiy_i when predicting y^i\hat{y}_i has zero covariance and zero optimism. A model that chases yiy_i, meaning one flexible enough to bend toward each individual point, has high covariance, and its training score is downward-biased by precisely that amount.

For a linear model with d parameters under squared loss, the sum evaluates in closed form:

optimism=2dσ2n \text{optimism} = \frac{2 d \sigma^{2}}{n}

Every parameter costs 2σ2/n2\sigma^2/n of self-deception. This single expression is the origin of Mallows’ CpC_p, of AIC, and of the whole tradition of penalising model size, and it says something the modern practitioner should keep: training error understates true error by an amount that grows with model complexity and shrinks with sample size. It is not merely optimistic, it is optimistic in a way that rewards exactly the models you should distrust most.

Hence the split. Evaluate on data the fitting procedure never touched, and the covariance term is zero by construction, because y^i\hat{y}_i cannot correlate with a yiy_i it has never seen.

The bias-variance decomposition, derived

Before examining estimators of RR, it is worth seeing what RR is made of. Assume y=f(x)+ϵy = f(x) + \epsilon with 𝔼[ϵ]=0\mathbb{E}[\epsilon] = 0 and Var(ϵ)=σ2\operatorname{Var}(\epsilon) = \sigma^2, and let f^\hat{f} be trained on a random dataset D. Fix a test point x and take expectations over both the noise and the training set.

Start by splitting off the noise:

𝔼[(yf^(x))2]=𝔼[(f(x)+ϵf^(x))2] \mathbb{E}\left[(y – \hat{f}(x))^2\right] = \mathbb{E}\left[(f(x) + \epsilon – \hat{f}(x))^2\right]

Expanding the square gives three terms:

=𝔼[(f(x)f^(x))2]+𝔼[ϵ2]+2,𝔼[ϵ(f(x)f^(x))]= \mathbb{E}\left[(f(x) – \hat{f}(x))^2\right] + \mathbb{E}[\epsilon^2] + 2,\mathbb{E}\left[\epsilon\left(f(x) – \hat{f}(x)\right)\right]

The test noise is independent of the training set and has mean zero, so the cross term vanishes and 𝔼[ϵ2]=σ2\mathbb{E}[\epsilon^2] = \sigma^2:

=𝔼[(f(x)f^(x))2]+σ2= \mathbb{E}\left[(f(x) – \hat{f}(x))^2\right] + \sigma^2

Now treat the first term the same way. Write f(x)=𝔼D[f^(x)]\bar{f}(x) = \mathbb{E}_D[\hat{f}(x)] for the model’s average prediction across all possible training sets, and add and subtract it:

𝔼[(f(x)f(x)+f(x)f^(x))2]\mathbb{E}\left[\left(f(x) – \bar{f}(x) + \bar{f}(x) – \hat{f}(x)\right)^2\right]

Expanding again, the cross term contains𝔼D[f(x)f^(x)]=0 \mathbb{E}_D[\bar{f}(x) – \hat{f}(x)] = 0 by the definition of f\bar{f}, so it disappears and what remains is:

𝔼[(yf^(x))2]=(f(x)f(x))2bias2+𝔼[(f^(x)f(x))2]variance+σ2irreducible \mathbb{E}\left[(y – \hat{f}(x))^2\right] = \underbrace{\left(f(x) – \bar{f}(x)\right)^2}{\text{bias}^2} + \underbrace{\mathbb{E}\left[\left(\hat{f}(x) – \bar{f}(x)\right)^2\right]}{\text{variance}} + \underbrace{\sigma^2}_{\text{irreducible}}

The decomposition is an identity, not an approximation, and the twice-repeated trick, adding and subtracting a mean so a cross term dies, is the entire proof. Its message for validation is the third term. No estimator, no architecture, and no quantity of data drives the error below σ2\sigma^2. When a model’s score plateaus, the honest question is whether you are near the noise floor or near the limits of the hypothesis class, and those are different diagnoses with different treatments.

The holdout estimator and how much it wobbles

Split the data, fit on one part, score on the other. For classification accuracy, each held-out prediction is either right or wrong, so the score is an average of Bernoulli trials with success probability p, the true accuracy:

p^=1mi=1m𝟙[f^(xi)=yi] \hat{p} = \frac{1}{m}\sum_{i=1}^{m} \mathbb{1}\left[\hat{f}(x_i) = y_i\right]

This estimator is unbiased, which is its whole appeal, and its variance is the standard Bernoulli result:

Var(p^)=p(1p)m,SE=p(1p)m\operatorname{Var}(\hat{p}) = \frac{p(1 – p)}{m}, \qquad \text{SE} = \sqrt{\frac{p(1-p)}{m}}

Put numbers in it, because the numbers are sobering. A test set of 200 rows with an observed accuracy of 0.90 gives:

SE=0.9×0.12000.021\text{SE} = \sqrt{\frac{0.9 \times 0.1}{200}} \approx 0.021

and an approximate 95 percent interval of p^±1.96×SE\hat{p} \pm 1.96 \times \text{SE}, which is roughly:

[0.858,;0.942][0.858, ; 0.942]

The model you reported as “90 percent accurate” is, on this evidence, somewhere between 86 and 94. Two models scoring 0.90 and 0.92 on such a test set are statistically indistinguishable, and any narrative explaining why the second architecture is better is a story told about noise.

The m\sqrt{m} in the denominator sets the price of precision. Halving the interval requires quadrupling the test set:

SE1m\text{SE} \propto \frac{1}{\sqrt{m}}

which is why small test sets cannot be rescued by careful analysis. They simply do not contain the information.

Here the holdout’s dilemma appears. Every row moved into the test set to shrink that interval is a row removed from training, which worsens the model being measured. A single split forces a choice between measuring a good model badly and a worse model well.

K-fold cross-validation as an estimator

Cross-validation refuses the choice. Partition the data into k folds, train k times leaving one fold out each time, and score each held-out fold:

CVk=1kj=1kR^j\text{CV}k = \frac{1}{k}\sum{j=1}^{k} \hat{R}_j

Equivalently, and more usefully for analysis, every point is predicted exactly once by a model that did not see it:

CVk=1ni=1nL(yi,;f^κ(i)(xi)) \text{CV}k = \frac{1}{n}\sum{i=1}^{n} L\left(y_i, ; \hat{f}^{-\kappa(i)}(x_i)\right)

where κ(i)\kappa(i) is the fold containing point i. Every row serves as test data once and as training data k-1 times, so the estimate uses the full sample without any point ever scoring the model that fitted it.

What does this estimator estimate? Not quite what people assume. Each fold’s model is trained on n(k1)/kn(k-1)/k points, so CVk\text{CV}_k is an almost-unbiased estimate of the error of a model trained on that reduced amount, not on all n. Since more data usually means less error, the estimate is pessimistic for the model you eventually ship, which is fitted on everything:

𝔼[CVk];;R(f^n)typically\mathbb{E}\left[\text{CV}_k\right] ; \ge ; R\left(\hat{f}_n\right) \quad \text{typically}

The size of the bias is governed by the steepness of the learning curve at that sample size. With k = 5 each model sees 80 percent of the data, and with k = 10 it sees 90 percent, so ten folds are less biased than five. On a flat learning curve the difference is negligible; on a steep one, meaning small n, it is not. This is the bias half of the fold-count trade-off, and it points toward large k.

Why the reported standard error is wrong

Here is where practice and mathematics part company. Everyone computes scores.std() across the folds and divides by k\sqrt{k}, treating the fold scores as an independent sample. They are not independent, and the reason is structural: any two training sets in a k-fold scheme overlap in (k2)/(k1)(k-2)/(k-1) of their rows. The models are trained on nearly the same data, so their errors are correlated.

For an average of k quantities each with variance σ2\sigma^2 and average pairwise correlation ρ\rho, the variance is:

Var(CVk)=σ2k+k1k,ρ,σ2 \operatorname{Var}\left(\text{CV}_k\right) = \frac{\sigma^2}{k} + \frac{k-1}{k},\rho,\sigma^2

This is the same formula that governs random forests, and it is doing the same work here. The naive standard error keeps only the first term and silently sets ρ=0\rho = 0, which understates the true variance whenever the folds are positively correlated, and they always are. Reported confidence intervals around cross-validated scores are therefore too narrow, and the significance you think you have detected between two models is partly an artefact of ignoring the covariance.

The situation is worse than a correctable underestimate. Bengio and Grandvalet proved that there is no unbiased estimator of the variance of k-fold cross-validation, because the covariance structure depends on quantities the procedure cannot recover from the data it has. This is not a gap awaiting a clever fix; it is an impossibility result. The practical consequence is modest but firm: report the spread of the fold scores as a description of stability, which it honestly is, and do not dress it up as a confidence interval, which it is not.

Leave-one-out, and why more folds stop helping

Push k to its maximum, k = n, and every model trains on n-1 points. The bias essentially disappears, since each fold’s model is trained on almost the full dataset. Tempting, and yet leave-one-out cross-validation is notorious for high variance. The correlation formula explains why exactly.

Two LOOCV training sets differ by two rows out of n. The fitted models are nearly identical, so ρ1\rho \to 1, and substituting into the variance expression:

Var(CVn);;σ2as ρ1\operatorname{Var}\left(\text{CV}_n\right) ; \longrightarrow ; \sigma^2 \quad \text{as } \rho \to 1

The averaging has bought nothing. Averaging reduces variance only to the extent the things averaged are different, and LOOCV averages n nearly-identical models. The n-fold structure looks like a huge sample and behaves like one observation, which is the same lesson the ensembles article drew from the same formula: correlation, not count, sets the floor.

So fold count trades bias against variance from both directions. Small k means each model sees less data, raising bias, but the folds are less correlated, lowering variance. Large k means less bias and more correlation. Five and ten are conventional because they sit near the bottom of that trade for typical learning curves, not because of any theorem.

Stratification, and repetition

Two variance reductions are nearly free. Stratified folds preserve the class proportions in every split, removing the component of fold-to-fold variance that comes purely from one fold randomly receiving more positives than another. That variability is noise about the partition, not information about the model, and eliminating it makes the estimate more stable without touching its bias. On imbalanced data the effect is large, and on a rare-class problem an unstratified fold may contain no positives at all, at which point the fold’s recall is undefined rather than merely noisy.

Repeated k-fold does the other obvious thing: run the whole scheme several times with different random partitions and average. This reduces the variance contributed by the arbitrary choice of partition, though it cannot touch the variance from the finite sample itself. Ten repeats of five-fold is not equivalent to a fifty-times-larger dataset, and no resampling scheme conjures information the sample does not contain.

The multiple comparisons trap

Now the failure that ruins more projects than all the others combined. Suppose you use a validation estimate to choose among M candidate models. Each candidate’s score is its true performance plus noise. You select the maximum, and then report that maximum as the chosen model’s performance.

The winner won partly on merit and partly on luck, and the reported score contains both. For M estimates with standard error s, the expected maximum of the noise alone grows like:

𝔼[max1jMZj];;2lnM \mathbb{E}\left[\max_{1 \le j \le M} Z_j\right] ; \approx ; \sqrt{2 \ln M}

so the selection bias in the winner’s reported score is approximately:

bias;;s2lnM\text{bias} ; \approx ; s\sqrt{2 \ln M}

Numbers again. A grid search over 100 configurations, with each score carrying a standard error of 0.02, produces:

0.02×2ln1000.02×3.030.0610.02 \times \sqrt{2 \ln 100} \approx 0.02 \times 3.03 \approx 0.061

Six accuracy points of pure illusion, manufactured by the act of choosing. Your best configuration is not six points better than your average one; it is roughly average and lucky. And note the villainy of the 2lnM\sqrt{2 \ln M} term: it grows with the number of things you tried, so the more thorough your search, the more inflated the winner. Diligence is punished.

The validation set has become training data. Not for the model’s parameters, but for the choice of model, and Efron’s covariance is back in a new costume, since the selection now covaries with the validation noise.

The fix is nested cross-validation. An inner loop selects hyperparameters; an outer loop, whose data the selection never touched, estimates performance:

Nested CV=1Kj=1KR^j(model selected using only the other K1 folds)\text{Nested CV} = \frac{1}{K}\sum_{j=1}^{K} \hat{R}_j\left(\text{model selected using only the other } K-1 \text{ folds}\right)

The cost is K×M×K \times M \times inner folds fits, which is why it is skipped, and the alternative is honest and cheap: a final test set opened exactly once, after all decisions are made. The discipline that ends every capstone script, keep the test set locked until the last line, is not fastidiousness. It is the only defence against s2lnMs\sqrt{2 \ln M}.

What the mathematics tells you to do

Never quote training error, because 2dσ2/n2d\sigma^2/n tells you exactly how much it is lying and in which direction. Size your test set with the standard error formula before you run anything, because 200 rows buy an interval of roughly four points and no analysis afterwards will narrow it. Prefer five or ten folds, since LOOCV’s ρ1\rho \to 1 means its extra folds average nearly identical models. Stratify always, because it removes variance that carries no information. Report the fold spread as stability rather than as a confidence interval, because no unbiased variance estimator for k-fold exists. And treat every tuned score as inflated by roughly s2lnMs\sqrt{2 \ln M} until an untouched test set says otherwise.

The unifying idea is a change in posture. A score is not a measurement of the model, it is a draw from a sampling distribution, and the questions worth asking are about that distribution rather than the draw. Bias asks where it is centred, variance asks how wide it is, correlation between folds asks why averaging helps less than you hoped, and selection asks what you did to it by choosing. A practitioner who reads the number as a fact will chase noise for weeks. One who reads it as an estimate with a standard error will know, before the meeting, which of two models is genuinely better and which pair is a coin flip dressed as a result.

See you soon, Andrei.

View Comments (2)

Leave a Reply

Subscribe to My Newsletter

Subscribe to my email newsletter to get the latest posts delivered right to your email. Pure inspiration, zero spam.

Discover more from Discuss Data Science, Machine Learning and Analytics

Subscribe now to keep reading and get access to the full archive.

Continue reading