The Dimensionality Challenge in Statistical Overfitting — Epoche C1
The claim, and the ratio that drives it Fit an ordinary linear regression to 100 observations using 90 predictors, with the linear model exactly correct and the predictors genuinely informative, and the fitted equation will predict fresh data about ten times worse than the unavoidable noise floor. Nothing about that model is complex: it is a straight-line rule, the simplest object in the statistician's catalogue, and it is not misspecified. What ruins it is the ratio of predictors to observations. That is the phenomenon this essay is about, and it is the reason the standard account of overfitting — that overfitting is what happens when a model is too flexible for its data, and that a simpler model would generalise better — is incomplete rather than wrong. Some vocabulary first, since the argument is quantitative. Overfitting is the gap between how well a fitted rule does on the data used to fit it and how well it does on data drawn afresh from the same source. Write $N$ for the number of observations and $p$ for the number of predictors, or features, measured on each. The quantity that will do all the work below is their ratio, $$\gamma = \frac{p}{N},$$ the number of parameters estimated per observation available. High-dimensional means $\gamma$ is not small — it is a statement about the ratio, not about $p$ alone. Ten predictors and ten observations is a high-dimensional problem; ten predictors and a million observations is not. The thesis is that $\gamma$, and not the functional form of the model, is what governs the variance of the fitted rule, and that this holds even for the least flexible model available. The bias–variance decomposition, stated correctly Before the high-dimensional behaviour can be derived, the classical decomposition it modifies has to be written down exactly, because the usual quick statement of it contains an error worth flagging. Suppose the response is generated as $Y = f^{*}(x) + \varepsilon$, where $f^{*}$ is the true regression function — the conditional mean of $Y$ given the predictors, the best any rule could do — and $\varepsilon$ is noise with mean zero and variance $\sigma^2$, independent of the predictors. Let $\hat{f}_{\mathcal{D}}$ be the rule fitted to a random training sample $\mathcal{D}$ of size $N$. Then at a fixed point $x$, averaging over both the randomness of the training sample and the noise in a new observation at $x$, $$E\big[(Y - \hat{f}_{\mathcal{D}}(x))^2\big] = \sigma^2 + \big(E[\hat{f}_{\mathcal{D}}(x)] - f^{*}(x)\big)^2 + \operatorname{Var}\big(\hat{f}_{\mathcal{D}}(x)\big).$$ The three terms are the irreducible error, the squared bias, and the variance. The derivation is two applications of the same trick: write $Y - \hat{f} = \varepsilon + (f^{*} - E[\hat{f}]) + (E[\hat{f}] - \hat{f})$, square, and take expectations. Every cross term vanishes — the first because $\varepsilon$ has mean zero and is independent of the training sample, the third because $E[\hat f] - \hat f$ has mean zero by construction — leaving the three squares. The correction is this. The bias must be measured against $f^{*}(x)$, the true regression function, and not against $y$, a realised observation. Writing the bias term as $(E[\hat f(x)] - y)^2$, as compressed statements of the decomposition often do, both mixes a random quantity into a term meant to be deterministic and double-counts the noise, since $y$ already contains the $\varepsilon$ that the $\sigma^2$ term accounts for. The identity as written above is the one that is true. The classical reading of this identity is that bias falls and variance rises as a model is made more flexible, so there is an interior optimum. That reading is sound. What it does not say — and what the rest of this essay supplies — is where the variance term actually comes from. It is tempting to read "variance" as "flexibility of the functional form". For linear least squares the variance term can be computed in closed form, and flexibility of form is not what appears in the answer. Why the variance of a linear fit is governed by $p/N$ This section computes the variance term for the least flexible model there is, and finds $\gamma$ sitting in it. Take the linear model $Y = X\beta + \varepsilon$, where $X$ is the $N \times p$ matrix of predictor values, $\beta$ the vector of $p$ unknown coefficients, and the noise has variance $\sigma^2$ in each coordinate and is uncorrelated across observations. The ordinary least squares estimate and its covariance are $$\hat{\beta} = (X^{\top}X)^{-1}X^{\top}Y, \qquad \operatorname{Var}(\hat{\beta}) = \sigma^2 (X^{\top}X)^{-1}.$$ The covariance formula follows by substituting $Y = X\beta + \varepsilon$ into the estimator: the systematic part reproduces $\beta$ exactly, and the noise contributes $(X^{\top}X)^{-1}X^{\top}\varepsilon$, whose covariance is $(X^{\top}X)^{-1}X^{\top}\cdot \sigma^2 I \cdot X(X^{\top}X)^{-1} = \sigma^2 (X^{\top}X)^{-1}$. Now sum the variance of the fitted values over the training points. The vector of fitted values is $\hat{Y} = HY$ with $H = X(X^{\top}X)^{-1}X^{\top}$, the matrix that projects onto the $p$-dimensional column space of $X$. A projection matrix has every eigenvalue equal to $0$ or $1$, and the number of ones is the dimension of the space projected onto, so $\operatorname{tr}(H) = p$. Since $\operatorname{Var}(\hat{Y}) = \sigma^2 H$, summing the diagonal gives $$\sum_{i=1}^{N}\operatorname{Var}(\hat{y}_i) = \sigma^2\operatorname{tr}(H) = \sigma^2 p,$$ and dividing by $N$ gives the average variance per fitted value: $\sigma^2 p / N = \sigma^2\gamma$. The factor $p$ is there because a projection onto a $p$-dimensional space keeps exactly $p$ of the $N$ independent noise directions; the division by $N$ is there because the total is spread over $N$ points. This is the standard result behind Mallows' criterion and behind the treatment of optimism in Hastie, Tibshirani and Friedman (2009): the expected error of a correctly specified linear fit, measured at the