Hyperparameters for statisticians

Statisticians who move into machine learning often find the word hyperparameter strange: it looks like a new concept, full of acronyms and “automatic optimisation” libraries. It is not new. The \(\lambda\) of ridge regression, the bandwidth of a kernel density estimator, the span of loess and the smoothing parameter of a spline are all hyperparameters, and statistics has been choosing them for more than 50 years — with cross-validation, GCV, information criteria and restricted maximum likelihood. What machine learning added was scale: models with ten or twenty hyperparameters, each fit taking minutes or hours. This note defines the concept in a statistician's language, shows why a hyperparameter cannot be estimated together with the parameters, compares the classical selection criteria in a simulation, presents the optimisation methods — grid, random search, Bayesian optimisation and successive halving — and ends with the bias that appears when you pick the best of many. Every number comes from a reproducible script.

§1Parameter or hyperparameter?

Take Hoerl and Kennard's (1970) ridge regression. For a fixed \(\lambda \ge 0\), the coefficients are

There are two layers here. \(\boldsymbol\beta\) is chosen by the fitting criterion: given \(\lambda\), minimising (1) gives a single \(\hat{\boldsymbol\beta}_\lambda\). But \(\lambda\) cannot be chosen by the same criterion: minimise (1) over \(\lambda\) as well and the answer is always \(\lambda = 0\), because the penalty can only add to the sum. That is the operational definition this note uses:

Formally, choosing a hyperparameter is a bilevel problem: at the inner level the parameters are fitted to the training data; at the outer level the hyperparameter is chosen to minimise an error the training data do not see,

\(R(\lambda)\) is the prediction risk, and it is unknown. All of §3 is about estimating it; all of §4 is about minimising a function that is only known with noise and is expensive to evaluate.

The word comes from Bayesian statistics. In a hierarchical model, such as those of Lindley and Smith (1972), the data depend on parameters and the parameters have a prior, whose own parameters — prior means and variances — came to be called hyperparameters. The link is more than a name: §3.3 shows that the ridge \(\lambda\) is the hyperparameter of a normal prior, and that this gives a way to estimate it.

The distinction is not about the nature of the number but about how it enters the criterion. The Box–Cox \(\lambda\), for instance, is estimated by maximum likelihood together with the coefficients — the profile likelihood has an interior maximum, so it is an ordinary parameter. The order \(p\) of an AR model, the number of components in principal-components regression and the number of knots of a regression spline, on the other hand, are hyperparameters: increasing them always improves the training fit.

Statisticians already have a long list of them; machine learning has another:

Two practical differences separate the two halves of the table. Classical statistical hyperparameters tend to be few, continuous and cheap: one or two numbers, and the fit for each value is a closed-form formula. Machine-learning hyperparameters tend to be many, of mixed types and expensive: continuous (learning rate), integer (depth), categorical (loss function), sometimes conditional (the dropout rate exists only if there is dropout), and every evaluation means training the model again. That is what makes hyperparameter optimisation a subject of its own.

§2Why it cannot be estimated jointly

The clearest way to see the problem is through the singular value decomposition \(\mathbf X = \mathbf U\mathbf D\mathbf V^\top\), with singular values \(d_1 \ge \dots \ge d_p\). The fitted values are

Each principal direction of the data is shrunk by the factor \(d_k^2/(d_k^2+\lambda)\): low-variance directions (small \(d_k\)) shrink the most. The effective degrees of freedom \(\operatorname{df}(\lambda)\) run from \(p\) (least squares, \(\lambda = 0\)) to 0 (\(\lambda \to \infty\)) and are the common currency for comparing models of continuous complexity. Training error is monotone in \(\lambda\), because \(\hat{\boldsymbol\beta}_0\) minimises the unconstrained sum of squares. The prediction risk has the familiar U shape:

where \(\boldsymbol\Sigma\) is the covariance of the predictors and \(\|\mathbf a\|^2_{\boldsymbol\Sigma} = \mathbf a^\top\boldsymbol\Sigma\mathbf a\). The figure shows both curves for a simulated dataset with \(n = 50\) observations and \(p = 40\) correlated predictors — a setting in which least squares is very unstable.

Three things in the figure hold for any hyperparameter. First, training error (grey) is useless for choosing \(\lambda\): it recommends \(\lambda = 0\), the worst point of the true curve. Second, choosing well pays: least squares has more than twice the minimum error, and even “predict zero” beats it. Third, the estimated curves (cross-validation and GCV) follow the shape of the true one, but their minimum does not fall in the same place — the bottom of the U is flat, and estimation noise decides where the minimum lands. §3.4 measures how large that noise is.

§3What to optimise: estimating the prediction risk

There are three families of criteria for estimating \(R(\lambda)\) or a substitute for it. All are old friends of statistics.

3.1 Cross-validation, LOOCV and GCV

\(K\)-fold cross-validation (Stone, 1974) splits the data into \(K\) blocks, fits the model on \(K-1\) of them, measures the error on the block left out and averages. It is the most general criterion: it works for any model and any loss. The cost is \(K\) fits per value of \(\lambda\).

For linear smoothers — any method with \(\hat{\mathbf y} = \mathbf H_\lambda\mathbf y\), such as ridge, splines, loess and kernels — the case \(K = n\) (leave-one-out, LOOCV) has a closed form. Nothing needs refitting:

where \(h_{ii}\) is the leverage of observation \(i\). (It is the PRESS statistic divided by \(n\).) Craven and Wahba (1979) and Golub, Heath and Wahba (1979) replaced each \(h_{ii}\) by the average \(\operatorname{tr}(\mathbf H_\lambda)/n\), which gives generalised cross-validation:

GCV is invariant to rotations of the data (LOOCV is not) and reads simply: the residual sum of squares inflated by a factor that charges for the degrees of freedom spent.

3.2 Information criteria with effective degrees of freedom

The usual AIC counts parameters. For a smoother, the parameter count is replaced by \(\operatorname{df}(\lambda)\):

AIC is an asymptotic approximation that assumes \(\operatorname{df}\) is small relative to \(n\). When it is not — and with \(p = 40\) and \(n = 50\) it is not — AIC under-penalises complexity and chooses \(\lambda\) too small. The corrected version of Hurvich, Simonoff and Tsai (1998), \(\operatorname{AIC_c}\), was proposed precisely for choosing smoothing parameters and penalises much more as \(\operatorname{df}\) approaches \(n\). Mallows' \(C_p\), \(\operatorname{RSS}_\lambda + 2\hat\sigma^2\operatorname{df}(\lambda)\), belongs to the same family but needs an outside estimate of \(\sigma^2\).

3.3 Marginal likelihood: the hyperparameter as a variance component

This is the most “statistical” route and the least known outside statistics. Put a normal prior on the coefficients:

The posterior mean of \(\boldsymbol\beta\) is exactly \(\hat{\boldsymbol\beta}_\lambda\) from (1). But now \(\lambda = \sigma^2/\sigma^2_\beta\) is a variance ratio in a mixed model, and it can be estimated from the likelihood of the data with \(\boldsymbol\beta\) integrated out:

with \(\sigma^2\) profiled out at \(\hat\sigma^2_\lambda = \mathbf y^\top(\mathbf I + \mathbf X\mathbf X^\top/\lambda)^{-1}\mathbf y/n\). This is empirical Bayes (or “type II maximum likelihood”). When there are unpenalised fixed effects, such as the intercept and the linear parts of a spline, the right version is REML, which integrates those out too. It is what mgcv does by default since Wood (2011), and the same criterion fits the hyperparameters of a Gaussian process (§4.3). It needs no validation data, is smooth in \(\lambda\) and extends to many hyperparameters (one per model term) without the cost blowing up.

The price is one more assumption: (9) is optimal if the prior (8) is right. The simulation below tests what happens when it is and when it is not.

3.4 Comparing the criteria in a simulation

A thousand datasets with the design of the figure (n = 50, p = 40, AR(1) predictors with correlation 0.6) in two scenarios. In the first, the coefficients are fixed and decay geometrically with alternating signs — a vector “with structure”, nothing like a normal with common variance. In the second, the coefficients are drawn afresh for each dataset from \(N(0,\tau^2)\), exactly the prior (8), with \(\tau\) chosen to give the same signal variance. Because the data are simulated, the true risk \(R(\lambda)\) is known, and we can measure how much each criterion loses relative to the optimal \(\lambda\).

1000 datasets per scenario, n = 50, p = 40, AR(1) predictors with ρ = 0.6, σ = 1. Left, fixed coefficients that decay geometrically; right, coefficients drawn from \(N(0, \tau^2)\) with \(\tau = 0.222\) — exactly the prior the marginal likelihood assumes. chosen log₁₀ λ: median (quartiles) of the chosen λ. Ratio: median of (prediction error at the chosen λ − σ²) / (error at the optimal λ − σ²); 1 is perfect. ≤ 1.1: share of datasets where the excess error was within 10% of the minimum.

Five lessons, all backed by the literature:

1. No criterion dominates. With fixed coefficients, GCV, LOOCV and 10-fold tie at the top; the marginal likelihood overestimates \(\lambda\) (median \(10^{1.7}\) against the optimum \(10^{0.85}\)) and loses more. When the prior is true the order reverses: the marginal likelihood is within 10% of the optimum in 65% of the datasets. This is Wahba's (1985) result comparing GCV and generalised maximum likelihood for splines: the likelihood is better when the Bayesian model is right; GCV is more robust when it is not.
2. Uncorrected AIC fails with so many degrees of freedom relative to \(n\): it picks \(\lambda\) more than an order of magnitude below the optimum and has double (or more) the excess error. AICc corrects in the right direction but here over-corrects.
3. The chosen \(\lambda\) is a random variable — a very variable one. With fixed coefficients the standard deviation of the 10-fold \(\log_{10}\hat\lambda\) is about 1.1 decades, against 0.4 decades for the optimal \(\lambda\). Worse, the correlation between the chosen and the optimal \(\lambda\) across datasets is slightly negative (about −0.1). Härdle, Hall and Marron (1988) documented the same for bandwidths: automatic selectors are right on average but do not track the sample-to-sample fluctuations of the optimum.
4. Even so, the cost of missing is small. The median ratio is about 1.1 for the best criteria: the bottom of the U is flat, so missing \(\lambda\) by half a decade costs little. That is why hyperparameter choice rarely needs to be perfect — it needs to avoid the extremes.
5. Closed forms are gold. GCV and LOOCV cost a single fit (one SVD) for the whole \(\lambda\) grid and tie with 10-fold, which costs ten. Whenever the model is a linear smoother, use them.

§4How to optimise: the problem as experimental design

Once a criterion is chosen, what is left is to minimise a function \(f(\boldsymbol\lambda)\) — validation error as a function of a vector of hyperparameters. For a statistician the best way to see this problem is as an experiment: each evaluation of \(f\) is an expensive run (training the model), the response is noisy (the validation set is finite), there are no derivatives, and the goal is to find the optimal region in few runs. It is, almost word for word, the problem of Box and Wilson's (1951) response surface methodology and of the computer experiments of Sacks, Welch, Mitchell and Wynn (1989). The methods below are experimental designs, from the simplest to the most adaptive.

Before any method, one rule matters more than the choice of method: parametrise on the right scale. Hyperparameters such as \(\lambda\), the learning rate and the bandwidth act multiplicatively — going from 0.001 to 0.01 changes as much as going from 0.1 to 1 — so they should be searched in \(\log\lambda\). A uniform search over \([0.0001, 1]\) spends 90% of its evaluations between 0.1 and 1.

4.1 Grid search

Pick \(L\) values for each of the \(d\) hyperparameters and evaluate all \(L^d\) combinations. It is a full factorial. It works very well with one or two hyperparameters and is easy to parallelise and to visualise. But the cost grows exponentially in \(d\): with 3 levels, 6 hyperparameters already need 729 fits, and a budget of 64 allows only 2 levels per axis.

4.2 Random search

Draw \(N\) independent configurations from the search region (on the log scale). Bergstra and Bengio (2012) showed this usually beats a grid with the same budget, through a low effective dimension argument: typically only a few hyperparameters matter much. A \(3\times3\times3\) grid tests 27 combinations but only 3 distinct values on each axis; if one of the three axes does not matter, it wastes two thirds of the budget repeating the same 9 combinations of the other two. Random search tests 27 distinct values on every axis.

There is also a simple guarantee. If the optimum has a “good enough” neighbourhood occupying a fraction \(q\) of the search space, the probability that at least one of \(N\) random configurations lands in it is

Sixty random evaluations find, with 95% probability, a configuration in the best 5% — whatever the dimension. Designs that fill space more evenly, such as the Latin hypercube (McKay, Beckman and Conover, 1979) or Sobol sequences, keep this advantage and avoid clumps.

4.3 Bayesian optimisation

Grid and random search are non-adaptive designs: they fix every point before seeing any result. Bayesian optimisation is sequential: after each evaluation it fits a statistical model for \(f\) — the surrogate — and uses it to choose the next point. The most common surrogate is a Gaussian process, which a geostatistician will recognise as kriging:

where \(\mu_t\) is the kriging prediction after \(t\) evaluations and \(s_t\) its standard error — small near points already evaluated, large far from them. The next point maximises an acquisition criterion. The most used is the expected improvement of Jones, Schonlau and Welch (1998). With \(f^\ast\) the smallest value observed so far and \(z = (f^\ast - \mu_t)/s_t\),

The two terms are the balance between exploitation and exploration: the first is large where the prediction is good (low \(\mu_t\)); the second where uncertainty is large (high \(s_t\)). Expected improvement is exactly zero at points already evaluated and positive everywhere else, so the method never stops exploring.

Two details interest the statistician. First, the Gaussian process has hyperparameters of its own — the kernel's correlation lengths, one per dimension — and they are estimated by the marginal likelihood of §3.3. They are hyperparameters chosen by a Bayesian criterion in order to choose other hyperparameters. Second, there are alternatives to the Gaussian process: the TPE of Bergstra, Bardenet, Bengio and Kégl (2011), the default in Optuna, models the density of the hyperparameters among good configurations and among bad ones and picks the next point by the ratio of the two densities; it handles categorical and conditional variables better. Snoek, Larochelle and Adams (2012) popularised the Gaussian process for tuning neural networks.

4.4 Cheap evaluations: successive halving and Hyperband

For many models you can get an approximate evaluation for a fraction of the cost: train on part of the data, for fewer iterations or fewer epochs. Successive halving (Jamieson and Talwalkar, 2016) turns this into a tournament. With factor \(\eta = 3\):

1. Draw 81 configurations and evaluate each on a small budget (here, 36 of the 324 training points).
2. Keep the best third (27) and evaluate again with three times the budget (108 points).
3. Keep the best 9 and evaluate on the full data. Pick the winner.

Each round spends the same number of training points (\(81\times36 = 27\times108 = 9\times324\)), so the whole tournament costs the same as 27 full fits — but it examines 81 configurations instead of 27. The risk is eliminating early a configuration that only shines with a large budget (a small learning rate, for instance, starts slowly). Hyperband (Li, Jamieson, DeSalvo, Rostamizadeh and Talwalkar, 2018) guards against this by running several tournaments with different starting points. For a statistician the closest analogy is sequential designs with stopping and screening experiments: many treatments in a cheap first phase, only the promising ones in the expensive phase.

4.5 Other families

When the hyperparameters are continuous and the model is differentiable, the gradient of the validation error with respect to them can be computed by implicit differentiation of the optimality condition of the inner level of (2), and used for gradient descent. Evolutionary algorithms and population-based training keep a population of configurations that recombine and mutate. They are useful in niches (very large networks, many continuous hyperparameters) but beyond the scope of this note.

4.6 The methods side by side

To compare the methods on a problem with more than one hyperparameter, the script uses kernel ridge regression: \(n = 324\) training observations with 5 inputs, a non-linear true function \(f(\mathbf x) = \sin(1.5x_1) + 0.5x_2x_3 + 0.6x_4\) (the fifth input is noise) and error variance 0.25. The kernel mixes an RBF and a linear kernel,

and there are three hyperparameters: \(\lambda\) (penalty, \(10^{-5}\) to \(10\)), \(\gamma\) (inverse kernel width, \(10^{-3}\) to \(10\)) and \(a\) (linear-kernel weight, 0 to 1). Each method has a budget of 27 full fits, picks the configuration with the lowest error on 200 validation points, and that configuration is scored on a test set of 5,000 points. Everything is repeated on 100 datasets.

The figure shows the first replication and two things the table below confirms. The surface has a diagonal valley: \(\lambda\) and \(\gamma\) interact, and tuning one at a time (fix \(\gamma\), find the best \(\lambda\), then fix that \(\lambda\) and find the best \(\gamma\)) can stop far from the optimum. And the third hyperparameter barely matters: in a main-effects fit over 4,860 random configurations, \(\lambda\) explains 10% of the variation in log validation error, \(\gamma\) 13% and \(a\) 0.3% (the rest is the \(\lambda\)–\(\gamma\) interaction). That is why the 27-point grid collapses to 9 spots on the plane that matters.

To see Bergstra and Bengio's argument properly at work, the second block of the experiment uses an RBF kernel with one width per input, \(\exp\{-\sum_j\gamma_j(x_j - x'_j)^2\}\): 6 hyperparameters (\(\lambda\) and five \(\gamma_j\)), a budget of 64 fits and 40 replications. Now the grid can only have 2 levels per axis.

Kernel ridge regression, n = 324 training, 200 validation, 5,000 test; the lowest possible error (noise) is 0.249. Three-hyperparameter block (λ, γ, a): 100 replications, 27 full fits. Successive halving uses as many training points as 27 full fits. Six-hyperparameter block (λ and one width per input): 40 replications, 64 fits. “Beats random search”: share of replications whose test error was below random search's in the same replication.

The results tell a story worth reading carefully:

• With 3 hyperparameters and 27 evaluations, the methods nearly tie. Median test errors range from 0.350 to 0.359, a difference of about 2%, against a gap of more than 40% to the irreducible error (0.249). Bayesian optimisation is best but beats random search in only 62% of the replications. The grid did as well as random search here — because the optimum fell near one of its points, something you cannot know in advance. Successive halving examined three times as many configurations and landed in the middle.
• With 6 hyperparameters and 64 evaluations, the order is clear. The grid (2 levels per axis) is worst and loses to random search in 72% of the replications. Bayesian optimisation beats random search in all 40 replications, and its median error (0.292) is well below random search's (0.341). In a 6-dimensional space, learning from earlier evaluations makes a difference.
• Bayesian optimisation rediscovered the structure of the problem. The chosen widths (median \(\log_{10}\gamma_j\)) were −0.9 for \(x_1\) (the sine input, which wants a narrow kernel), about −1.9 for \(x_2\) and \(x_3\), −2.6 for \(x_4\) (which enters linearly, so it wants a wide kernel) and −3.0 for \(x_5\) — the lower limit of the search, which in practice switches off the pure-noise input. That is “automatic relevance determination”, obtained as a by-product of the optimisation.

§5The winner's curse

There is a hidden bias in this whole process, and a statistician will spot it from afar. The validation error of each configuration is a noisy estimate. Taking the minimum of many noisy estimates gives an estimate biased downwards: the winner tends to be a configuration that is good and lucky on the validation set. If \(\hat f_1,\dots,\hat f_k\) are unbiased estimates of \(f_1,\dots,f_k\),

and the gap grows with \(k\) and with the variance of each \(\hat f_j\). It is the same phenomenon as post-selection inference and regression to the mean. In the experiment:

Means over 100 replications; the k configurations are the first k of one random search of 243. With k = 1 there is no selection and the optimism is zero within Monte Carlo error (the “±” is one standard error). It grows with k: the more configurations compete, the more the winner owes its win to noise in the validation set.

The optimism is small here — about 6% of the error with 243 configurations — because the validation set has 200 points and the configurations are highly correlated with each other. With smaller validation sets, thousands of configurations or optimisers that search the surface more aggressively, it can be much larger. Cawley and Talbot (2010) showed cases where the selection bias is of the same order as the difference between the algorithms being compared, and Varma and Simon (2006) showed that reporting the cross-validation error of the selected model can make a classifier with no signal at all look useful.

The remedy is the usual one: the error you report must come from data that played no part in the choice. In practice, either a separate test set, touched once at the end, or nested cross-validation: an outer \(K\)-fold loop to estimate the error and, inside each outer fold, a complete inner loop that chooses the hyperparameters using only that fold's training data. The cost is multiplicative (\(K_{\text{outer}}\times K_{\text{inner}}\times\) number of configurations), but it is the only way to estimate the error of the whole procedure, tuning included.

§6A practical guide

1. Start from the defaults. The defaults of mature libraries were chosen to work reasonably on many problems. Tune first what the literature or documentation names as most important (in boosting: learning rate and number of trees; in random forests: almost nothing).
2. Use the log scale for anything positive that acts multiplicatively, and choose wide limits. If the optimum lands on the edge of the region, widen the region.
3. If there is a closed-form criterion, use it. GCV, LOOCV via (5) and REML cost one fit. That covers ridge, splines, GAMs and mixed models.
4. Choose the optimiser by budget and dimension. One or two cheap hyperparameters: a fine grid. Several hyperparameters, moderate evaluations: random search or a Latin hypercube with a few dozen points. Expensive evaluations and up to about 20 continuous dimensions: Bayesian optimisation. The ability to evaluate on part of the data or part of the iterations: successive halving or Hyperband.
5. Prefer the simplest model within the noise. The one-standard-error rule (Breiman, Friedman, Olshen and Stone, 1984) picks the most regularised model whose cross-validation error is within one standard error of the minimum. It acknowledges that the minimum is noisy (§3.4) and trades a small loss in fit for stability. It is lambda.1se in glmnet.
6. Report the error on data that did not see the choice (§5) and document the search region, budget, criterion and seeds. Tuning is part of the model; a result without this information is not reproducible.

In code, with the most common tools:

# R: ridge/lasso with CV and the one-standard-error rule; GAM with REML
library(glmnet)
cv <- cv.glmnet(X, y, alpha = 0, nfolds = 10)     # alpha = 0: ridge
coef(cv, s = "lambda.1se")

library(mgcv)
m <- gam(y ~ s(x1) + s(x2) + s(x3), method = "REML")   # one lambda per term, by REML

# Python: random search, successive halving and Bayesian optimisation (TPE)
from scipy.stats import loguniform
from sklearn.experimental import enable_halving_search_cv  # noqa: F401
from sklearn.model_selection import RandomizedSearchCV, HalvingRandomSearchCV, cross_val_score
from sklearn.kernel_ridge import KernelRidge

space = {"alpha": loguniform(1e-5, 10), "gamma": loguniform(1e-3, 10)}   # log scale
search = RandomizedSearchCV(KernelRidge(kernel="rbf"), space, n_iter=60, cv=5,
                            scoring="neg_mean_squared_error", random_state=0).fit(X, y)
tourney = HalvingRandomSearchCV(KernelRidge(kernel="rbf"), space, factor=3,
                                resource="n_samples", cv=5, random_state=0).fit(X, y)

import optuna
def objective(trial):
    m = KernelRidge(kernel="rbf", alpha=trial.suggest_float("alpha", 1e-5, 10, log=True),
                    gamma=trial.suggest_float("gamma", 1e-3, 10, log=True))
    return -cross_val_score(m, X, y, cv=5, scoring="neg_mean_squared_error").mean()
study = optuna.create_study(sampler=optuna.samplers.TPESampler(seed=0))
study.optimize(objective, n_trials=60)

The script behind this note implements everything from scratch in numpy — ridge via the SVD, the six criteria of §3, kernel ridge regression, the Gaussian process with expected improvement and successive halving — so every formula can be checked against the code.

§7Summary

1. A hyperparameter defines which estimator is used and cannot be chosen by the fitting criterion, which would push it to the most complex extreme. Choosing it is a bilevel problem: fit on training data, evaluate outside it.
2. Statistics already has many: the \(\lambda\) of ridge and lasso, bandwidths, spans, smoothing parameters, \(k\). The word comes from Bayesian hierarchical models — and the ridge \(\lambda\) literally is the hyperparameter of a normal prior.
3. To estimate prediction risk there is cross-validation (general, expensive), LOOCV and GCV (closed form for linear smoothers), information criteria with effective degrees of freedom (AIC fails when \(\operatorname{df}\) is large relative to \(n\)) and marginal likelihood / REML (optimal when the prior is right).
4. In the simulation no criterion dominated: GCV and CV were best with fixed coefficients; the marginal likelihood with coefficients drawn from the prior. The chosen \(\lambda\) varied much more than the optimum and barely correlated with it — but because the bottom of the curve is flat, the typical loss was about 10% of the excess error.
5. Hyperparameter optimisation is experimental design: grid (factorial), random search (better when few hyperparameters matter), Bayesian optimisation (kriging with expected improvement) and successive halving (sequential screening). With 3 hyperparameters the methods tied; with 6, Bayesian optimisation beat random search in 100% of the replications and the grid came last.
6. The winner's validation error is optimistic (winner's curse). Report the error on an untouched test set or with nested cross-validation.

§8Further reading

• A. E. Hoerl and R. W. Kennard (1970). Ridge regression: biased estimation for nonorthogonal problems. Technometrics, 12(1), 55–67.
• D. V. Lindley and A. F. M. Smith (1972). Bayes estimates for the linear model. Journal of the Royal Statistical Society B, 34(1), 1–41.
• M. Stone (1974). Cross-validatory choice and assessment of statistical predictions. Journal of the Royal Statistical Society B, 36(2), 111–147.
• P. Craven and G. Wahba (1979). Smoothing noisy data with spline functions. Numerische Mathematik, 31, 377–403.
• G. H. Golub, M. Heath and G. Wahba (1979). Generalized cross-validation as a method for choosing a good ridge parameter. Technometrics, 21(2), 215–223.
• G. Wahba (1985). A comparison of GCV and GML for choosing the smoothing parameter in the generalized spline smoothing problem. Annals of Statistics, 13(4), 1378–1402.
• W. Härdle, P. Hall and J. S. Marron (1988). How far are automatically chosen regression smoothing parameters from their optimum? Journal of the American Statistical Association, 83(401).
• C. M. Hurvich, J. S. Simonoff and C.-L. Tsai (1998). Smoothing parameter selection in nonparametric regression using an improved Akaike information criterion. Journal of the Royal Statistical Society B, 60(2), 271–293.
• S. N. Wood (2011). Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. Journal of the Royal Statistical Society B, 73(1), 3–36.
• G. E. P. Box and K. B. Wilson (1951). On the experimental attainment of optimum conditions. Journal of the Royal Statistical Society B, 13(1), 1–45.
• M. D. McKay, R. J. Beckman and W. J. Conover (1979). A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics, 21(2), 239–245.
• J. Sacks, W. J. Welch, T. J. Mitchell and H. P. Wynn (1989). Design and analysis of computer experiments. Statistical Science, 4(4), 409–423.
• D. R. Jones, M. Schonlau and W. J. Welch (1998). Efficient global optimization of expensive black-box functions. Journal of Global Optimization, 13, 455–492.
• J. Bergstra and Y. Bengio (2012). Random search for hyper-parameter optimization. Journal of Machine Learning Research, 13, 281–305.
• J. Bergstra, R. Bardenet, Y. Bengio and B. Kégl (2011). Algorithms for hyper-parameter optimization. NeurIPS 24.
• J. Snoek, H. Larochelle and R. P. Adams (2012). Practical Bayesian optimization of machine learning algorithms. NeurIPS 25.
• K. Jamieson and A. Talwalkar (2016). Non-stochastic best arm identification and hyperparameter optimization. AISTATS.
• L. Li, K. Jamieson, G. DeSalvo, A. Rostamizadeh and A. Talwalkar (2018). Hyperband: a novel bandit-based approach to hyperparameter optimization. Journal of Machine Learning Research, 18(185), 1–52.
• S. Varma and R. Simon (2006). Bias in error estimation when using cross-validation for model selection. BMC Bioinformatics, 7, 91.
• G. C. Cawley and N. L. C. Talbot (2010). On over-fitting in model selection and subsequent selection bias in performance evaluation. Journal of Machine Learning Research, 11, 2079–2107.
• C. E. Rasmussen and C. K. I. Williams (2006). Gaussian Processes for Machine Learning, ch. 5. MIT Press. Marginal likelihood for hyperparameters.
• T. Hastie, R. Tibshirani and J. Friedman (2009). The Elements of Statistical Learning, 2nd ed., ch. 7. Springer.
• M. Feurer and F. Hutter (2019). Hyperparameter optimization. In F. Hutter, L. Kotthoff and J. Vanschoren (eds), Automated Machine Learning, ch. 1. Springer. A complete, free review.