Hyperparameter Search and Nested Cross-Validation

A search over 32 configurations returns one winner and a tempting number: its best validation score. That number helped choose the winner, so reporting it as though the configuration had been fixed in advance ignores selection. We need evaluation data that the search has not inspected.

A fitted model contains parameters learned from data, such as regression coefficients or neural-network weights. Many important choices are made before that fitting begins: the degree of a polynomial, the strength of a regularization penalty, the number of neighbors in a nearest-neighbor model, or the learning rate of an optimizer. These choices are hyperparameters. They define the model family, the fitting procedure, or both, and they are not updated by the inner optimization that fits the parameters.

A validation score is a noisy function of the hyperparameters, so choosing the best score is itself a learning procedure. If we use the same validation folds both to select hyperparameters and to report performance, the winner's score is usually too optimistic. The search has inspected many candidates and picked the one that happened to look best on those particular folds. Nested cross-validation separates the two jobs: an inner loop selects hyperparameters, and an outer loop estimates how well the entire selection-and-fit procedure generalizes.

In this chapter we compare grid search and random search under equal evaluation budgets, explain why selection bias appears, and build nested cross-validation from primitive loops. We then apply the ideas to a polynomial ridge regression model, assess it with outer folds, and compare equal-budget searches across repeated seeds. The chapter builds on the previous discussion of validation strategies and data leakage and on the classification metrics from the preceding chapter, where we saw that a score is only meaningful when we know exactly which data produced it.

Recording the search space, fold rule, budget, and random seeds makes the selection procedure explicit, so it can be inspected and rerun with the same data and code. Instead of trying a few settings by hand and remembering the results, we define a search space, a validation rule, and a budget. Every candidate is evaluated on the same folds, and the selected configuration can be traced back to a specific score. This is valuable when several people collaborate on a model, when a pipeline must be rerun later, or when a regulator or reviewer asks how the final configuration was chosen.

Grid search is simple and exhaustive over a finite grid. When the search space is small and mostly discrete, it is easy to audit and easy to parallelize. Random search is often more efficient when only a few hyperparameters matter, because it spends less effort varying irrelevant dimensions. Nested cross-validation adds an honest estimate of the complete procedure, including the search step. It also exposes instability. If different outer folds select very different hyperparameters, we learn that the selection rule is not stable, which is useful information that a single best score hides.

The computational cost of search is multiplicative. If we evaluate $B$ candidates with $K_{\mathrm{inner}}$ inner folds across $J$ outer folds, we fit roughly $J \times B \times K_{\mathrm{inner}}$ models before any final refit. Grid search grows exponentially with the number of tuned dimensions, and random search can waste draws in uninformative regions. Grid search finds the smallest observed validation loss on its finite grid, but that need not identify the lowest true risk over the wider search space. Random search may miss that configuration entirely. Neither method removes validation noise.

Nested cross-validation offers a way to assess selection when no untouched test set is available. If we already have a clean test set that was never used for development, a single outer holdout is often enough and is cheaper. Nested cross-validation is most useful when we need to use the same data for both selection and assessment, or when we want to compare complete model-selection procedures rather than one fixed model. Ordinary random-fold evaluation assumes independent observations representative of the target population. For time series, grouped data, or clustered observations, ordinary random folds can leak information and produce misleadingly optimistic estimates.

Search can overfit the inner validation folds, especially when the budget is large relative to the amount of validation data. A large search space also invites researcher degrees of freedom: if we keep expanding the space after seeing results, the reported score is no longer a clean estimate of a pre-specified procedure. Finally, preprocessing must be inside the cross-validation loop. Scaling, imputation, and data-driven feature selection learn from data and must be fitted within each training fold. A fixed polynomial feature map does not estimate population quantities, but its degree is selected inside the search, and the full pipeline keeps that selection and learned scaling within the correct splits.

It helps to view hyperparameter search as a second learning problem wrapped around the first. The inner problem fits model parameters for a fixed hyperparameter setting. The outer problem compares those settings using held-out data. The formulas below make both layers explicit, because the distinction between them is exactly what nested cross-validation is designed to protect.

Let the training data be

$$D=\{(x_i,y_i)\}_{i=1}^{n},$$

where $x_i$ is a feature vector and $y_i$ is a target. A learning algorithm $A$ takes a training set and a hyperparameter vector $\lambda$ and returns a fitted model $f_\lambda$. The hyperparameters live in a search space $\Lambda$, where:

• $\lambda$: a hyperparameter vector, an element of the search space
• $\Lambda$: the search space, the set of all allowed hyperparameter values
• $A$: the learning algorithm that maps a training set and $\lambda$ to a fitted model
• $f_\lambda$: the fitted model returned by $A$ for hyperparameters $\lambda$

In a polynomial ridge pipeline, for example,

$$\lambda=(d,\alpha),$$

where:

• $d$: the polynomial degree (a discrete hyperparameter)
• $\alpha$: the ridge penalty strength (a continuous hyperparameter)

The degree controls the feature map

$$\phi_d(x)=(x,x^2,\ldots,x^d),$$

where:

• $\phi_d(x)$: the polynomial feature map that sends a scalar input $x$ to a vector of its powers up to degree $d$
• $d$: the polynomial degree, the highest power included in the feature map
• $x$: a scalar input feature value

Note that $\phi_d(x)$ is a $d$-dimensional vector for scalar input $x$, with components $x, x^2, \ldots, x^d$, so the coefficient vector $\beta$ also has length $d$, ensuring the inner product is well defined:

$$\phi_d(x_i)^\top\beta = \sum_{p=1}^{d} \beta_{p} x_i^{p}$$

The dot product is a scalar prediction, where:

• $\phi_d(x_i)$: the $d$-dimensional polynomial feature vector for observation $i$, with components $x_i, x_i^2, \ldots, x_i^d$
• $\beta$: the coefficient vector of length $d$ that weights each power
• $\beta_p$: the coefficient associated with the $p$-th power of the input
• $x_i$: the feature value for observation $i$ (a scalar in this univariate polynomial case)
• $d$: the polynomial degree, the highest power included in the sum

This feature map excludes a constant column. The implementation adds an unpenalized intercept and standardizes each power using training-fold means $\mu_p$ and positive scales $s_p$. Write the standardized power as $z_{ip}=(x_i^p-\mu_p)/s_p$ and the standardized feature vector as $z_i=(z_{i1},\ldots,z_{id})$, where:

• $z_{ip}$: the standardized version of the $p$-th power of input $x_i$
• $\mu_p$: the mean of $x^p$ computed on the training fold only
• $s_p$: the positive scale of $x^p$ (typically its standard deviation) computed on the training fold only
• $z_i$: the vector of standardized polynomial features for observation $i$

For a constant column, use scale one. Ridge trades residual error against coefficient size on this standardized scale:

$$\left(\hat\beta_0(d,\alpha),\hat\beta_\alpha\right)=\arg\min_{\beta_0,\beta}\left(\sum_{i=1}^{n}\left(y_i-\beta_0-z_i^\top\beta\right)^2+\alpha\lVert\beta\rVert_2^2\right),$$

where:

• $\hat\beta_0(d,\alpha)$: the unpenalized intercept estimate, a function of the degree $d$ and penalty $\alpha$
• $\hat\beta_\alpha$: the fitted (penalized) coefficient vector at penalty strength $\alpha$; it implicitly depends on the standardized feature basis determined by $d$
• $\beta$: the coefficient vector being optimized over, of dimension $d$ (the same length as the standardized feature vector $z_i$)
• $n$: the number of training observations
• $y_i$: the observed target for observation $i$
• $\beta_0$: the unpenalized intercept, not penalized by the ridge term
• $z_i$: the polynomial features standardized using only the training observations
• $\beta_0+z_i^\top\beta$: the model's predicted value for observation $i$
• $\alpha$: the ridge penalty strength (controls the size of the shrinkage)
• $\lVert\beta\rVert_2^2$: the squared L2 norm of the coefficients, $\sum_{p=1}^{d}\beta_{p}^2$

We will study linear regression and polynomial regression in detail later; for now, the important point is that $d$ changes the model family and $\alpha$ changes the fitting objective, but neither is learned by the inner least-squares step. Changing $\lambda$ changes the model family or the fitting rule. Note that the ridge objective above shows the unpenalized intercept $\beta_0$ alongside the penalized coefficients, so the prediction used in the loss includes both.

The quantity we would like to know is the generalization risk

$$R(\lambda\mid D)=\mathbb{E}_{(X,Y)\sim P}\left[L\left(Y,\,A(D;\lambda)(X)\right)\right],$$

where:

• $R(\lambda\mid D)$: the expected loss of the particular model fitted on $D$ with setting $\lambda$, where $\lambda$ is a hyperparameter vector in the search space $\Lambda$
• $D$: the training data, held fixed in this conditional risk
• $(X,Y)$: a new observation drawn from the data-generating distribution
• $P$: the data-generating distribution
• $L$: a loss such as squared error
• $A(D;\lambda)$: the fitted model returned by the learning algorithm on training data $D$ with setting $\lambda$
• $X$: a new input feature vector drawn from the data-generating distribution
• $A(D;\lambda)(X)$: the fitted model's prediction at the new input $X$

and the squared-error loss is

$$L(y,\hat y)=(y-\hat y)^2,$$

where:

• $L(y,\hat y)$: the squared-error loss for a true target $y$ and prediction $\hat y$
• $y$: the observed target
• $\hat y$: the model's predicted value

The new observation must be independent of the training data. In this conditional expression, fitting randomness is fixed; the expectation averages over new evaluation observations. A procedure-level risk additionally averages this conditional risk over fresh training samples of a specified size and over any fitting randomness. These are different targets: cross-validation evaluates fits trained on fewer rows than a fit on all $n$ observations. Below, $R(\lambda)$ abbreviates the procedure-level risk at the training size being evaluated, with that size held fixed across candidates. We estimate risk with held-out data, once per candidate.

For $K$-fold cross-validation, partition the training indices into folds $V_1,\ldots,V_K$. Let $D^{(-k)}$ be the data with fold $V_k$ removed, and let

$$f_\lambda^{(-k)}=A\left(D^{(-k)};\lambda\right)$$

be the model fitted without fold $k$. The cross-validation estimate averages the loss over all folds and, within each fold, over the held-out observations, so that every observation contributes once to the estimate:

$$\hat R_{\mathrm{CV}}(\lambda) = \frac{1}{K}\sum_{k=1}^{K} \frac{1}{|V_k|} \sum_{i\in V_k} L\left(y_i,f_\lambda^{(-k)}(x_i)\right)$$

where:

• $\hat R_{\mathrm{CV}}(\lambda)$: the $K$-fold cross-validated estimate of the risk at hyperparameters $\lambda$
• $K$: the number of folds
• $k$: the index of the fold currently held out
• $V_k$: the index set of observations in fold $k$
• $|V_k|$: the number of observations in fold $k$
• $f_\lambda^{(-k)}$: the model fitted on all folds except $k$, with hyperparameters $\lambda$
• $L$: the per-observation loss
• $y_i$: the observed target for observation $i$
• $x_i$: the feature vector for observation $i$

The inner sum averages loss within held-out fold $k$, and the outer sum gives each fold equal weight. If fold sizes differ, this need not equal the pooled per-observation mean loss. The superscript $(-k)$ is important: the model used to predict observations in $V_k$ never saw those observations during fitting. For a fixed $\lambda$, this is a standard validation estimate. The subtlety arises because we will compare many such estimates and keep the smallest one.

Parameters are the numbers that the fitting algorithm optimizes. In a linear model, they are the regression coefficients; in a neural network, they are the weights and biases. Hyperparameters are the numbers or structural choices that control how that optimization is set up. They include the polynomial degree, the ridge penalty, the number of folds in cross-validation, the learning rate of gradient descent, the number of trees in a forest, and the kernel width in a kernel method. The distinction is practical: parameters are estimated by minimizing a training objective, while hyperparameters are chosen by comparing validation performance or by using prior knowledge.

A useful way to keep the distinction straight is to ask what changes when a value changes:

• A parameter changes during fitting on a fixed training set. The fitting rule itself stays the same.
• A hyperparameter changes the fitting rule, the model family, or the search protocol. It is fixed while the parameters are optimized.
• A search space specifies which hyperparameter values are allowed. It is part of the procedure and must be fixed before the final evaluation.

The search space $\Lambda$ can contain discrete values, continuous intervals, categorical choices, or conditional choices. For a polynomial ridge model we might use the degree and the penalty:

$$d\in\{1,2,\ldots,8\}, \qquad \alpha\in[10^{-3},1],$$

where:

• $d$: the polynomial degree (a discrete hyperparameter)
• $\alpha$: the ridge penalty strength (a continuous hyperparameter)

For a positive penalty spanning several orders of magnitude, a logarithmic search allocates resolution to relative changes in $\alpha$. Moving from $0.001$ to $0.01$ changes the penalty by a factor of ten, just as moving from $0.1$ to $1$ does. The interval $[0,0.01]$ occupies only 1% of the linear range $[0,1]$. A coarse linear grid can therefore miss substantial relative differences among small positive penalties. Search spaces should be chosen before looking at the final test result, and they should be wide enough to include reasonable models but not so wide that the budget is spread too thin.

Once the search space is defined, the next question is how to spend the evaluation budget. Grid search evaluates a Cartesian product of candidate values. If the search space has $H$ hyperparameters and the $j$-th hyperparameter has $m_j$ candidate values, the grid is

$$\Lambda_{\mathrm{grid}} = \left\{ \left(\lambda_1^{(a_1)},\ldots,\lambda_H^{(a_H)}\right) : a_1\in\{1,\ldots,m_1\},\ldots,a_H\in\{1,\ldots,m_H\} \right\},$$

where:

• $\lambda_j^{(a_j)}$: the $a_j$-th candidate value for the $j$-th hyperparameter
• $m_j$: the number of candidate values for the $j$-th hyperparameter
• $H$: the number of hyperparameters being tuned
• $a_j$: the index selecting the candidate value along the $j$-th hyperparameter dimension

Each element of $\Lambda_{\mathrm{grid}}$ is a single combination, chosen by fixing one candidate index $a_j$ per dimension. The number of candidates is the size of the Cartesian product, obtained by multiplying the number of candidate values along each hyperparameter dimension:

$$B=\prod_{j=1}^{H} m_j,$$

where:

• $B$: the total number of grid candidates
• $H$: the number of hyperparameters being tuned
• $m_j$: the number of candidate values for the $j$-th hyperparameter
• $\prod_{j=1}^{H}$: the product over all $H$ hyperparameters
• $j$: the index identifying the hyperparameter whose candidate count is multiplied

If we tune eight polynomial degrees and four ridge penalties, the budget is $8\times 4=32$ candidate configurations. Each configuration is evaluated on the same inner folds. Grid candidate enumeration is deterministic and exhaustive over the specified grid. Randomness in model fitting or fold construction still needs to be controlled. Its weakness is the product: five hyperparameters with ten values each require $10^5=100{,}000$ candidates, which is often impractical. The grid also evaluates every combination, including combinations whose differences may be irrelevant to performance.

Random search draws candidates independently from a distribution $p(\lambda)$:

$$\lambda^{(b)}\sim p(\lambda), \qquad b=1,\ldots,B.$$

where:

• $\lambda^{(b)}$: the $b$-th candidate hyperparameter vector drawn by the search
• $p(\lambda)$: the search distribution over the hyperparameter space
• $B$: the total number of candidates (the budget)
• $b$: the index of the candidate draw

We then select

$$b^*=\arg\min_{b\in\{1,\ldots,B\}}\hat R_{\mathrm{CV}}\left(\lambda^{(b)}\right), \qquad \hat\lambda=\lambda^{(b^*)}.$$

where:

• $b^*$: the selected candidate index, using the first index to break an exact tie
• $\hat\lambda$: the selected candidate vector at that index
• $\hat R_{\mathrm{CV}}\left(\lambda^{(b)}\right)$: the $K$-fold cross-validation estimate of the risk for the $b$-th candidate

For the positive ridge penalty in this chapter, we choose a log-uniform distribution because it gives equal probability to equal-width intervals on the logarithmic scale, including each full order of magnitude within its support. Random search does not enumerate combinations. Independent draws often produce many different values along continuous dimensions; discrete values can repeat, and distinctness is not guaranteed. If only one or two hyperparameters strongly affect performance, random search can cover their important ranges much faster than a grid that also varies irrelevant dimensions. It is stochastic, however, so the result depends on the seed, and two runs with the same budget can select different configurations.

An equal-budget comparison must hold the data, inner folds, loss, and number of candidate evaluations fixed. Grid search and random search should not be compared when one method receives more fits or different validation splits. For random search, we should report the seed or repeat the search across several seeds, because the candidate draw is part of the procedure. The fair question is not "which method found a lower score once?" but "how does each complete search procedure behave under the same budget and protocol?"

Let $\Lambda_B$ denote the finite set of evaluated candidates within the larger allowed space $\Lambda$. After evaluating those candidates, we choose the one with the smallest validation loss:

$$\hat\lambda = \arg\min_{\lambda\in\Lambda_B} \hat R_{\mathrm{CV}}(\lambda).$$

The reported score of the winner is therefore

$$\hat R_{\mathrm{CV}}(\hat\lambda) = \min_{\lambda\in\Lambda_B} \hat R_{\mathrm{CV}}(\lambda).$$

This is a minimum of noisy estimates. Even if all candidates have the same true risk, the one that happens to receive a favorable validation sample will look better than the others. To see the mechanism, write the validation error of candidate $\lambda$ as

$$\hat R_{\mathrm{CV}}(\lambda) = R(\lambda)+\varepsilon_\lambda,$$

where:

• $\hat R_{\mathrm{CV}}(\lambda)$: the cross-validated estimate of the risk at hyperparameters $\lambda$
• $R(\lambda)$: the true risk of the model fitted with hyperparameters $\lambda$
• $\varepsilon_\lambda$: the estimation error on these particular folds (the gap between the validation estimate and the true risk)

The selected candidate has

$$\hat R_{\mathrm{CV}}(\hat\lambda) = R(\hat\lambda)+\varepsilon_{\hat\lambda},$$

where:

• $\hat R_{\mathrm{CV}}(\hat\lambda)$: the validation estimate of the selected candidate's risk
• $R(\hat\lambda)$: the true risk of the selected candidate
• $\varepsilon_{\hat\lambda}$: the estimation error of the selected candidate on these folds, the gap between its validation estimate and its true risk (defined as $\hat R_{\mathrm{CV}}(\lambda)-R(\lambda)$ evaluated at $\hat\lambda$)

Selection favors candidates with favorable estimation errors. Under the unbiased-candidate assumptions stated below, the resulting optimism is nonnegative in expectation; this is not a claim about every realized score. Defining optimism as the expected amount by which the selected model's validation estimate falls below its true risk (equivalently, the expected overstatement of performance), we have

$$\begin{aligned} \mathrm{optimism} &= \mathbb{E}\left[R(\hat\lambda)\right] - \mathbb{E}\left[\hat R_{\mathrm{CV}}(\hat\lambda)\right] \\ &= -\mathbb{E}\left[\varepsilon_{\hat\lambda}\right] \end{aligned}$$

where:

• $\mathrm{optimism}$: the expected overstatement of performance caused by selecting the best candidate
• $\hat\lambda$: the selected hyperparameters (a random quantity, depending on the validation folds)
• $R(\hat\lambda)$: the true risk of the model chosen by the search
• $\hat R_{\mathrm{CV}}(\hat\lambda)$: the validation estimate of the selected model's risk
• $\mathbb{E}[\cdot]$: the expectation over the randomness of the validation folds

For a fixed finite candidate set, suppose each candidate estimate is unbiased for its fixed target risk: $\mathbb{E}[\hat R_{\mathrm{CV}}(\lambda)]=R(\lambda)$. The minimum is a concave function of its arguments. Applying Jensen's inequality to this concave function gives $\mathbb{E}[\min_\lambda \hat R_{\mathrm{CV}}(\lambda)]\le\min_\lambda \mathbb{E}[\hat R_{\mathrm{CV}}(\lambda)]$, and substituting the unbiasedness $\mathbb{E}[\hat R_{\mathrm{CV}}(\lambda)]=R(\lambda)$ yields $\mathbb{E}[\min_\lambda\hat R_{\mathrm{CV}}(\lambda)]\le\min_\lambda R(\lambda)$. Also, over the fixed finite candidate set $\Lambda_B$, $\min_{\lambda\in\Lambda_B} R(\lambda)\le R(\hat\lambda)$ for every realization (since $\hat\lambda\in\Lambda_B$), so $\min_{\lambda\in\Lambda_B} R(\lambda)\le\mathbb{E}[R(\hat\lambda)]$. Combining these gives $\mathbb{E}[\min_{\lambda\in\Lambda_B}\hat R_{\mathrm{CV}}(\lambda)]\le\mathbb{E}[R(\hat\lambda)]$, so the expected optimism is nonnegative. Independence across candidates is unnecessary for this inequality. Actual CV estimates need not be unbiased for the risk of a model refitted on all data; the training-size difference is another source of discrepancy. The inequality is not a guarantee that a realized outer score exceeds its inner score.

A simple calculation makes the size of the effect concrete. Suppose two candidates have the same true risk, and their validation estimates are independent normal random variables with the same mean $R_0$ and positive standard deviation $\sigma$. The normality and independence assumptions are stated explicitly because the exact value below depends on them; normality is an illustrative model assumption rather than a general property of validation estimates, and shared folds can make candidate estimates correlated. Write the standardized estimates as $Z_1, Z_2 \sim \mathcal{N}(0,1)$, independently. The expected minimum is

$$\begin{aligned} \mathbb{E}[\min(Z_1, Z_2)] &= \int_{-\infty}^{\infty} z \cdot 2\,\varphi(z)\,\Phi(-z)\, dz \\ &= -\frac{1}{\sqrt{\pi}} \\ &\approx -0.564 \end{aligned}$$

where:

• $z$: the value at which we evaluate whether one draw is the smaller
• $\varphi(z)$: the standard normal density
• $\Phi(z)$: the standard normal cumulative distribution function, the probability that a standard normal draw is at most $z$; its exceedance probability is $1-\Phi(z)=\Phi(-z)$
• the factor $2\,\varphi(z)\,\Phi(-z)$: the density that $z$ is the smaller of the two draws (one draw at $z$, the other above it, in either order)

To evaluate the integral, use the identity $\mathbb{E}[\min(Z_1, Z_2)] = -\mathbb{E}\lvert Z_1 - Z_2\rvert / 2$. With $W = Z_1 - Z_2 \sim \mathcal{N}(0, 2)$, we have $\sigma_W = \sqrt{2}$ and $\mathbb{E}\lvert W\rvert = \sqrt{2/\pi}\,\sigma_W$ for $W \sim \mathcal{N}(0, \sigma_W^2)$. This follows from the half-normal expectation

$$\begin{aligned} \mathbb{E}|W| &= 2\int_0^\infty w \cdot \frac{1}{\sqrt{2\pi}\sigma_W} e^{-w^2/(2\sigma_W^2)} \, dw \\ &= \sqrt{\frac{2}{\pi}}\,\sigma_W && \text{(half-normal mean for } W\sim\mathcal{N}(0,\sigma_W^2)\text{)} \\ &= \sqrt{\frac{2}{\pi}}\,\sqrt{2} && \text{(substitute } \sigma_W=\sqrt{2}\text{)} \\ &= \frac{2}{\sqrt{\pi}}. \end{aligned}$$

where:

• $w$: the integration variable running over nonnegative values
• $\sigma_W$: the standard deviation of the zero-mean normal $W$ (equal to $\sqrt{2}$ here, since $W = Z_1 - Z_2$)
• $e^{-w^2/(2\sigma_W^2)}$: the Gaussian kernel evaluated at $w$, with the normalizing factor $1/(\sqrt{2\pi}\sigma_W)$ outside
• the factor $2$: accounts for the two symmetric halves of the normal density when integrating only over $w \ge 0$

Therefore

$$\begin{aligned} \mathbb{E}[\min(Z_1, Z_2)] &= -\frac{\mathbb{E}|Z_1 - Z_2|}{2} \\ &= -\frac{1}{2} \cdot \sqrt{\frac{2}{\pi}} \cdot \sqrt{2} \\ &= -\frac{1}{2} \cdot \frac{2}{\sqrt{\pi}} \\ &= -\frac{1}{\sqrt{\pi}}. \end{aligned}$$

Rescaling back to the candidates' units, the winner's validation score is about $0.564\sigma$ below the common mean even though neither candidate is better. With more candidates, the minimum is pushed even lower. More precisely, if $M$ candidates share the same true risk and their validation estimates have standard deviation $\sigma$, the expected minimum score is

$$R_0+\sigma\, \mathbb{E}[\min(Z_1,\ldots,Z_M)],$$

where:

• $R_0$: the common true risk shared by all $M$ candidates
• $M$: the number of candidates sharing that common true risk
• $\sigma$: the standard deviation of each candidate's validation estimate
• $Z_m$: the standardized validation error of candidate $m$, defined as $Z_m=(\hat R_m-R_0)/\sigma$, an independent standard normal draw

This expression is obtained by writing each estimate as $\hat R_m = R_0 + \sigma Z_m$ with $Z_m=(\hat R_m-R_0)/\sigma$ standard normal, so the minimum is $R_0 + \sigma \min_m Z_m$ and linearity of expectation gives the expression above. Since $\mathbb{E}[\min(Z_1,\ldots,Z_M)]<0$ for $M\ge 2$ under these assumptions, the expected minimum lies below the common risk $R_0$.

This expression decreases as $M$ grows because the standardized minimum $\mathbb{E}[\min(Z_1,\ldots,Z_M)]$ decreases with $M$. This is the same "winner's curse" that can arise when we select the maximum of noisy measurements.

The winner's curse can be quantified without fitting a model. For $M$ independent candidate estimates with the same true risk, the expected best observed score is $R_0$ plus $\sigma$ times the expected minimum of $M$ standard normal draws. The next cell approximates that expected minimum by numerical integration over $[-8,8]$, and the figure plots the negative value as expected selection optimism in standard-deviation units.

Nested cross-validation addresses this by adding an outer loop. Partition the data into outer folds $V_1^{(o)},\ldots,V_J^{(o)}$, and let $D^{(-j)}$ denote the data with outer fold $V_j^{(o)}$ removed. For each outer fold $j$, run the entire inner search on $D^{(-j)}$ and obtain

$$\hat\lambda_j = \arg\min_{\lambda\in\Lambda_B} \hat R_{\mathrm{CV}}^{(j)}(\lambda),$$

where:

• $\hat\lambda_j$: the hyperparameters selected by the inner search when outer fold $j$ is held out
• $\hat R_{\mathrm{CV}}^{(j)}(\lambda)$: the $K$-fold inner cross-validation estimate of candidate $\lambda$ computed only on $D^{(-j)}$

For random search, $\Lambda_B$ here means the candidates drawn for this outer training run. Candidate generation must also exclude the outer test targets.

Then fit the selected model on the outer training set,

$$f_j=A\left(D^{(-j)};\hat\lambda_j\right),$$

and evaluate it on the outer test fold. The nested estimate averages the outer-fold losses, each computed with the model whose hyperparameters were selected by that fold's inner search:

$$\hat R_{\mathrm{nested}} = \frac{1}{J}\sum_{j=1}^{J} \frac{1}{|V_j^{(o)}|} \sum_{i\in V_j^{(o)}} L\left(y_i,f_j(x_i)\right)$$

where:

• $\hat R_{\mathrm{nested}}$: the nested cross-validation estimate of the procedure's risk
• $J$: the number of outer folds
• $D^{(-j)}$: the data with outer fold $j$ removed, used to run the inner search and select $\hat\lambda_j$
• $V_j^{(o)}$: the index set of observations in outer fold $j$
• $|V_j^{(o)}|$: the number of observations in outer fold $j$
• $f_j$: the model refitted on the outer training set using the selected hyperparameters $\hat\lambda_j$
• $L$: the per-observation loss
• $y_i$: the observed target for observation $i$
• $x_i$: the feature vector for observation $i$
• $i$: the index running over observations within an outer test fold
• $j$: the index of the outer fold

The outer test fold was not used to choose $\hat\lambda_j$, so it gives a cleaner view of how the selection procedure performs on unseen data. Notice that $\hat\lambda_j$ can differ from fold to fold. The nested estimator is not the cross-validation score of one fixed model; it is an estimate of the risk of the complete procedure that takes a training set, runs a search, and returns a fitted model. With independent identically distributed observations and a procedure fixed before inspecting outer results, outer held-out losses evaluate that procedure at the outer-training sample size, approximately $(1-1/J)n$. They do not directly estimate the conditional risk of the final model refitted on all $n$ rows. Finite data and overlapping outer training sets still create variability. To see why the distinction matters, note that $\hat R_{\mathrm{nested}}$ averages over outer folds each of which selected its own $\hat\lambda_j$; keeping each outer test fold out of its entire selection process is what avoids reusing its scores for selection, but it also means the nested score is not the score of any single deployable model.

After nested cross-validation, we usually want a final model for deployment. The standard recipe is to run the search on all available development data, obtain $\hat\lambda$, and refit

$$\hat f=A(D;\hat\lambda).$$

The nested estimate describes the procedure; the final model is the procedure applied to all development data. If an untouched test set exists, we evaluate $\hat f$ on it once. That test score is often the cleanest estimate of final performance, and in that situation nested cross-validation is not mandatory. Nested cross-validation is the appropriate tool when we need to estimate the procedure using the same data that were available for development, or when we want to compare several complete procedures under a common protocol.

Reporting should include the search method, the search space, the budget, the inner and outer fold counts, the metric, the random seeds, and the selected hyperparameters across outer folds. A single number hides too much. It is useful to report the mean outer score, the standard deviation across outer folds, and a standard error across repeated nested runs. If we repeat the entire nested procedure $M\ge2$ times with different outer splits, we can summarize the repetition-level scores with a mean and an estimated standard error. The mean of the per-repetition nested estimates is

$$\bar R=\frac{1}{M}\sum_{m=1}^{M}\hat R_m,$$

and the estimated standard error of that mean is

$$\widehat{\mathrm{SE}}(\bar R) = \frac{1}{\sqrt{M}} \sqrt{\frac{1}{M-1} \sum_{m=1}^{M} \left(\hat R_m-\bar R\right)^2}$$

where:

• $M$: the number of repeated nested runs (each with a fresh outer split)
• $m$: the index of the repetition
• $\hat R_m$: the nested estimate produced by the $m$-th repetition
• $\bar R$: the mean nested estimate across repetitions
• $\widehat{\mathrm{SE}}(\bar R)$: the estimated standard error of $\bar R$ across repetitions
• $\left(\hat R_m-\bar R\right)^2$: the squared deviation of repetition $m$'s estimate from the mean
• $\sqrt{\frac{1}{M-1}\sum_{m=1}^{M}\left(\hat R_m-\bar R\right)^2}$: the sample standard deviation of the repetition-level scores

The formula is the sample standard deviation of the repetition-level scores divided by $\sqrt{M}$: the inner square root is the sample standard deviation across repetitions, and dividing by $\sqrt{M}$ converts that spread into a standard error for their mean.

With independently drawn partition and search seeds, this is a Monte Carlo standard error for the mean over that randomization, conditional on this fixed dataset. It can decrease as we add repetitions without adding observations. The repetition-level standard deviation describes split/search variability; neither quantity measures how results vary over new datasets or supplies a population-risk confidence interval. If the selected degree or penalty changes drastically from fold to fold, the procedure is unstable, and a simpler model or more data may be preferable.

Grid and random search treat validation loss as a black-box objective and do not use derivatives. A discrete degree cannot be updated by an ordinary continuous gradient. Some continuous hyperparameters, including ridge penalties, admit differentiable optimization under suitable conditions, but that is a different method from the searches implemented here. A validation surface need not be convex, and finite-data noise can affect which candidate wins.

Random search has a simple coverage property. The following formula makes the probability of covering the good region explicit for an independent-draw search. Suppose the good region has probability mass $q$ under the chosen search distribution. Each independent draw then has probability $q$ of producing an acceptable candidate; the value of $q$ is simply the probability mass the search distribution assigns to the good region, and a uniform distribution over the search region makes it equal to the good region's volume (or length) divided by the total volume (or length). The probability that $B$ independent draws contain at least one good candidate is

$$P(\text{at least one hit})=1-(1-q)^B,$$

where:

• $P(\text{at least one hit})$: the probability that at least one of the $B$ candidates falls in the good region
• $q$: the per-draw probability of an acceptable candidate under the search distribution
• $B$: the number of independent candidate draws
• $(1-q)^B$: the probability that all $B$ independent draws miss the good region, since each draw misses with probability $1-q$ and draws are independent

For a target hit probability $p$ with $0<p<1$ and per-draw hit probability $q$ with $0<q<1$, the smallest integer budget giving hit probability at least $p$ is

$$B_{\min}=\left\lceil\frac{\log(1-p)}{\log(1-q)}\right\rceil,$$

where:

• $B_{\min}$: the smallest integer number of independent draws needed to reach hit probability at least $p$, obtained by rounding the real-valued lower bound up to the next integer
• $p$: the target lower bound on the probability of at least one hit (with $0<p<1$)
• $q$: the probability mass of the good region under the search distribution (with $0<q<1$)
• $\lceil\cdot\rceil$: the ceiling function, which rounds its argument up to the smallest integer greater than or equal to it

To derive the budget bound, start from the requirement that the hit probability reach $p$ and rearrange the steps one at a time:

$$\begin{aligned} 1-(1-q)^B &\ge p && \text{(require hit probability at least } p \text{)} \\ (1-q)^B &\le 1-p && \text{(rearrange)} \\ B\log(1-q) &\le \log(1-p) && \text{(take natural log of both sides; } \log \text{ is increasing)} \\ B &\ge \dfrac{\log(1-p)}{\log(1-q)} && \text{(divide by } \log(1-q)<0 \text{, reversing the inequality)} \end{aligned}$$

Because $1-q\in(0,1)$, we have $\log(1-q)<0$, so the division used in the last step reverses the inequality direction. Because $1-p\in(0,1)$ as well, $\log(1-p)<0$, and the ratio of two negative quantities is positive, so the lower bound on $B$ is a positive number. The ceiling rounds the real-valued lower bound upward. The limiting cases confirm the bound: if $q=0$, no finite positive budget can hit; if $q=1$, one draw already suffices since $1-(1-q)^B=1\ge p$ for any $B\ge 1$. This formula explains why random search can be effective when the good region is not too small. It also shows the limits: if $q$ is tiny, the required budget becomes enormous.

The coverage formula can be read as a budget-planning curve. The next cells evaluate $1-(1-q)^B$ over a range of budgets for three per-draw hit probabilities, then plot the resulting hit probabilities.

Each outer fold performs $B$ candidates times $K_{\mathrm{inner}}$ inner fits, plus one refit on the outer training set, giving

$$\begin{aligned} N_{\mathrm{fits}} &= J\left(BK_{\mathrm{inner}}+1\right) && \text{(total fits including one refit per outer fold)} \\ C &\approx N_{\mathrm{fits}}\,T_{\mathrm{fit}} && \text{(approximate total cost)} \\ &= J\left(BK_{\mathrm{inner}}+1\right)\,T_{\mathrm{fit}} && \text{(substitute } N_{\mathrm{fits}}\text{)} \end{aligned}$$

where:

• $N_{\mathrm{fits}}$: the number of fits for nested assessment, including one outer refit per fold, given by $J(BK_{\mathrm{inner}}+1)$
• $C$: the approximate fitting cost of the nested procedure, $C\approx N_{\mathrm{fits}}T_{\mathrm{fit}}$
• $J$: the number of outer folds
• $B$: the number of candidate configurations searched per outer fold
• $K_{\mathrm{inner}}$: the number of inner folds used to score each candidate
• $T_{\mathrm{fit}}$: the average time to fit one model
• $\left(BK_{\mathrm{inner}}+1\right)$: the per-outer-fold fit count, namely $B$ candidates times $K_{\mathrm{inner}}$ inner fits plus one refit

The $JBK_{\mathrm{inner}}$ search fits and $J$ outer refits make up this count. A final search and refit on all development data add $BK_{\mathrm{inner}}+1$ fits. With $J=5$, $B=32$ and $K_{\mathrm{inner}}=4$, nested assessment uses $645$ fits per repetition and a final search/refit adds $129$. Different candidates can have different fit times, so equal fit counts need not mean equal wall-clock time. Whether this cost is practical depends on the candidate budget, fold counts, per-fit cost, available compute, and the value of assessing the complete selection procedure. Parallelizing over candidates or folds can reduce wall-clock time, but it does not reduce the total number of fits.

We first import the libraries used for the arithmetic and the later search implementation. The hidden setup seeds the shared plotting environment once.

The first block builds one row per outer fold. Each row contains the three inner scores, the selected degree, the best inner score, and the outer score. The selection step is an argmin: in each row, we compare the three inner scores and keep the smallest. For example, if the inner scores are $0.48$, $0.42$, and $0.41$, the selected degree is the third one and the best inner score is $0.41$. That is the choice the inner loop would make.

The second step records, for each outer fold, the outer MSE of the selected configuration refitted on the outer training set and scored on the outer test fold; those values are supplied in the outer_mse column. The key distinction is summarized here:

1. The inner score answers: "Among these candidates, which one looks best on the inner validation folds?"
2. The outer score answers: "After choosing that candidate, how well does the fitted model perform on data that were not used for the choice?"

The printed table lets us compare those two answers fold by fold. For each outer fold $j$, let $\hat R_{\mathrm{CV},j}(\lambda)$ be the inner cross-validation score of candidate $\lambda$ computed on the training data for outer fold $j$, and let $\hat R_{\mathrm{outer},j}$ be the model's score on outer test fold $j$. The inner mean is the average of the best inner scores, where the minimum is over the evaluated candidate set $\Lambda_B$ and the inner cross-validation is computed only on the training portion of outer fold $j$:

$$\bar R_{\mathrm{inner}} = \frac{1}{J}\sum_{j=1}^{J} \min_{\lambda\in\Lambda_B} \hat R_{\mathrm{CV},j}(\lambda)$$

where:

• $\bar R_{\mathrm{inner}}$: the mean of the best inner cross-validation scores across outer folds
• $J$: the number of outer folds
• $\Lambda_B$: the finite set of evaluated candidate hyperparameters
• $\hat R_{\mathrm{CV},j}(\lambda)$: the inner cross-validation score for candidate $\lambda$ on the training portion of outer fold $j$

The outer mean is the average of the outer test scores:

$$\bar R_{\mathrm{outer}} = \frac{1}{J}\sum_{j=1}^{J} \hat R_{\mathrm{outer},j}.$$

where:

• $\bar R_{\mathrm{outer}}$: the mean of the outer test-fold scores across outer folds
• $J$: the number of outer folds
• $\hat R_{\mathrm{outer},j}$: the model's score on outer test fold $j$

The difference between the outer and inner means estimates the inner-outer gap for this table:

$$\widehat{\mathrm{gap}} = \bar R_{\mathrm{outer}}-\bar R_{\mathrm{inner}},$$

where the outer mean uses each fold's held-out score while the inner mean uses that fold's best-selected candidate, so their difference is a descriptive inner-outer gap:

• $\widehat{\mathrm{gap}}$: the estimated gap between the outer and inner means
• $\bar R_{\mathrm{inner}}$: the mean of the best inner-fold scores across outer folds
• $\bar R_{\mathrm{outer}}$: the mean of the outer-fold scores across outer folds
• $J$: the number of outer folds
• $\hat R_{\mathrm{CV},j}(\lambda)$: the inner cross-validation score for candidate $\lambda$ on outer fold $j$, computed only on the training portion of outer fold $j$
• $\hat R_{\mathrm{outer},j}$: the model's score on outer test fold $j$

The difference is positive when the outer mean exceeds the selected inner mean. Selection, different fitting sample sizes, and finite-sample variation all contribute; this gap does not isolate selection bias. The output reports this gap. It also reports the outer standard deviation and the standard error

$$\widehat{\mathrm{SE}}=\frac{s_{\mathrm{outer}}}{\sqrt{J}}$$

where:

• $\widehat{\mathrm{SE}}$: the independence-based standard-error benchmark for this illustrative table
• $s_{\mathrm{outer}}$: the sample standard deviation of the outer-fold scores
• $J$: the number of outer folds

The fold standard deviation describes spread among these five values. Dividing it by $\sqrt{J}$ would be a standard-error estimate for independent identically distributed fold scores, but overlapping training sets invalidate that justification for ordinary nested CV. Treat the printed outer SE as an independence benchmark, not a calibrated uncertainty interval.

With the bookkeeping clear, we can implement the complete procedure on a real synthetic dataset. We generate 180 observations from a cubic polynomial with Gaussian noise, then fit a polynomial ridge regression model. The pipeline is

$$\text{polynomial features} \;\longrightarrow\; \text{standardization} \;\longrightarrow\; \text{ridge regression}.$$

Standardization is inside the pipeline, so it is fit only on the training portion of each split. This is not a cosmetic detail: if we standardized the full dataset before splitting, the validation folds would influence the training preprocessing through the mean and scale. Keeping the scaler inside the pipeline makes the entire preprocessing step part of the model that is refit within every fold.

The hyperparameters are the polynomial degree $d\in\{1,\ldots,8\}$ and the ridge penalty $\alpha\in\{10^{-3},10^{-2},10^{-1},1\}$. The grid therefore has $8\times 4=32$ candidates. We use 5 outer folds, 4 inner folds, and 10 repeated outer splits. The functions below implement the search from primitive loops so that every selection step is visible: candidate generation, inner cross-validation, selection, outer refitting, and outer scoring.

We use a SeedSequence hierarchy to give each repetition its own child sequence, then separate outer splits, inner splits, primitive candidate draws, and library candidate draws within that repetition. NumPy designs these child streams to be independent with very high probability; they are reproducible pseudorandom streams, not a proof of physical independence. We rebuild the same hierarchy when comparing methods, so grid and random search share outer and inner folds. Candidate samplers have separate streams. We pass derived 32-bit integer seeds to the scikit-learn splitters and random search. The finite seed protocol is checked for collisions; distinct hierarchy paths alone would not guarantee distinct 32-bit integers.