Hyperparameter Optimization: Foundations, Algorithms, Best Practices and Open Challenges

Supervised ML addresses the problem of inferring a model from labeled training data that is then used to predict data from the same underlying distribution with minimal error. Let $\mathcal{D}$ be a labeled data set with $n$ observations, where each observation $(\mathbf{x}^{(i)},y^{(i)})$ consists of a $p$-dimensional feature vector²²2More generally, with a slight abuse of notation, the feature vector $\mathbf{x}^{(i)}$ could be taken as a tensor (for example when the feature is an image). $\mathbf{x}^{(i)}\in\mathcal{X}$ and its label $y^{(i)}\in\mathcal{Y}$. Hence, we define the data set³³3More precisely, $\mathcal{D}$ is an indexed tuple, but we will continue to use common terminology and call it a data set. $\mathcal{D}=\left(\left(\mathbf{x}^{(1)},y^{(1)}\right),\ldots,\left(\mathbf{x}^{(n)},y^{(n)}\right)\right)$. We assume that $\mathcal{D}$ has been sampled i.i.d. from an underlying, unknown distribution, so $\mathcal{D}\sim(\mathds{P}_{xy})^{n}$. Each dimension of a $p$-dimensional $\mathbf{x}$ will usually be of numerical, integer, or categorical type. While some ML algorithms can handle all of these data types natively (e.g., tree-based methods), others can only work on numeric features and require encoding techniques for categorical types. The most common supervised ML tasks are regression and classification, where $\mathcal{Y}=\mathds{R}$ for regression and $\mathcal{Y}$ is finite and categorical for classification with $|\mathcal{Y}|=g$ classes. Although we mainly discuss HPO in the context of regression and classification, all covered concepts easily generalize to other supervised ML tasks, such as Poisson regression, survival analysis, cost-sensitive classification, multi-output tasks, and many more. An ML model is a function $f:\mathcal{X}\rightarrow\mathds{R}^{g}$ that assigns a prediction in $\mathds{R}^{g}$ to a feature vector from $\mathcal{X}$. For regression, $g$ is $1$, while in classification the output represents the $g$ decision scores or posterior probabilities of the $g$ candidate classes. Binary classification is usually simplified to $g=1$, with a single decision score in $\mathds{R}$ or only the posterior probability for the positive class. The function space – usually parameterized – to which a model belongs is called the hypothesis space and denoted as $\mathcal{H}$.

The goal of supervised ML is to fit a model given $n$ observations sampled from $\mathds{P}_{xy}$, so that it generalizes well to new observations from the same data generating process. Formally, an ML learner or inducer $\mathcal{I}$ configured by HPs $\bm{\lambda}\in\Lambda$ maps a data set $\mathcal{D}$ to a model $\hat{f}$ or equivalently to its associated parameter vector $\bm{\hat{\theta}}$, i.e.,

where $\mathds{D}:=\bigcup_{n\in\mathds{N}}(\mathcal{X}\times\mathcal{Y})^{n}$ is the set of all data sets. While model parameters $\bm{\hat{\theta}}$ are an output of the learner $\mathcal{I}$, HPs $\bm{\lambda}$ are an input. We also write $\mathcal{I}_{\bm{\lambda}}$ for $\mathcal{I}$ or $\hat{f}_{\mathcal{D},\bm{\lambda}}$ for $\hat{f}$ if we want to stress that the inducer was configured with $\bm{\lambda}$ or that the model was learned on $\mathcal{D}$ by an inducer configured by $\bm{\lambda}$. A loss function $L:\mathcal{Y}\times\mathds{R}^{g}\rightarrow\mathds{R}^{+}_{0}$ measures the discrepancy between the prediction and the true label. Many ML learners use the concept of empirical risk minimization (ERM) in their training routine to produce their fitted model $\hat{f}$, i.e., they optimize $\mathcal{R}_{\text{emp}}(f)$ or $\mathcal{R}_{\text{emp}}(\bm{\theta})$ over all candidate models $f\in\mathcal{H}$

on the training data $\mathcal{D}$ (c.f. Figure 1). This empirical risk is only a stochastic proxy for what we are actually interested in, namely the theoretical risk or true generalization error $\mathcal{R}(f):=\mathds{E}_{(\mathbf{x},y)\sim\mathds{P}_{xy}}[L\left(y,f(\mathbf{x})\right)]$. For many complex hypothesis spaces, $\mathcal{R}_{\text{emp}}(f)$ can become considerably smaller than its true risk $\mathcal{R}(f)$. This phenomenon is known as overfitting, which in ML is usually addressed by either constraining the hypothesis space or regularized risk minimization, i.e., adding a complexity penalty $J(\bm{\theta})$ to (2).

After training an ML model $\hat{f}$, a natural follow-up step is to evaluate its future performance given new, unseen data. We seek to use an unbiased, high-quality statistical estimator, which numerically quantifies the performance of our model when it is used to predict the target variable for new observations drawn from the same data-generating process $\mathds{P}_{xy}$.

Most learners are highly configurable by HPs, and their generalization performance usually depends on this configuration in a non-trivial and subtle way. HPO algorithms automatically identify a well-performing HPC $\bm{\lambda}\in\tilde{\Lambda}$ for an ML algorithm $\mathcal{I}_{\bm{\lambda}}$. The search space $\tilde{\Lambda}\subset\Lambda$ contains all considered HPs for optimization and their respective ranges:

where $\tilde{\Lambda}_{i}$ is a bounded subset of the domain of the $i$-th HP $\Lambda_{i}$, and can be either continuous, discrete, or categorical. This already mixed search space can also contain dependent HPs, leading to a hierarchical search space: An HP $\lambda_{i}$ is said to be conditional on $\lambda_{j}$ if $\lambda_{i}$ is only active when $\lambda_{j}$ is an element of a given subset of $\Lambda_{j}$ and inactive otherwise, i.e., not affecting the resulting learner [135]. Common examples are kernel HPs of a kernelized machine such as the SVM, when we tune over the kernel type and its respective hyperparameters as well. Such conditional HPs usually introduce tree-like dependencies in the search space, and may in general lead to dependencies that may be represented by directed acyclic graphs.

The general HPO problem as visualized in Figure 2 is defined as:

where $\bm{\lambda}^{*}$ denotes the theoretical optimum, and $c(\bm{\lambda})$ is a shorthand for the estimated generalization error when $\mathcal{I},\mathcal{J},\rho$ are fixed. We therefore estimate and optimize the generalization error $\widehat{\mathrm{GE}}(\mathcal{I},\mathcal{J},\rho,\bm{\lambda})$ of a learner $\mathcal{I}_{\bm{\lambda}}$, w.r.t. an HPC $\bm{\lambda}$, based on a resampling split $\mathcal{J}=\left((J_{\mathrm{train},1},J_{\mathrm{test},1}),\dots,(J_{\mathrm{train},B},J_{\mathrm{test},B})\right)$.⁵⁵5Note again that optimizing the resampling error will result in biased estimates, which are problematic when reporting the generalization error; use nested CV for this, see Section 4.4. Note that $c(\bm{\lambda})$ is a black-box, as it usually has no closed-form mathematical representation, and hence no analytic gradient information is available. Furthermore, the evaluation of $c(\bm{\lambda})$ can take a significant amount of time. Therefore, the minimization of $c(\bm{\lambda})$ forms an expensive black-box optimization problem.

Taken together, these properties define an optimization problem of considerable difficulty. Furthermore, they rule out many popular optimization methods that require gradients or entirely numerical search spaces or that must perform a large number of evaluations to converge to a well-performing solution, like many meta-heuristics. Furthermore, as $c(\bm{\lambda})=\widehat{\mathrm{GE}}(\mathcal{I},\mathcal{J},\rho,\bm{\lambda})$, which is defined via resampling and evaluates $\bm{\lambda}$ on randomly chosen validation sets, $c$ should be considered a stochastic objective – although many HPO algorithms may ignore this fact or simply handle it by assuming that we average out the randomness through enough resampling replications.

We can thus define the HP tuner $\tau:(\mathcal{D},\mathcal{I},\tilde{\Lambda},\rho)\mapsto\hat{\bm{\lambda}}$ that proposes its estimate $\hat{\bm{\lambda}}$ of the true optimal configuration $\bm{\lambda}^{*}$ given a dataset $\mathcal{D}$, an inducer $\mathcal{I}$ with corresponding search space $\tilde{\Lambda}$ to optimize, and a target measure $\rho$. The specific resampling splits $\mathcal{J}$ used can either be passed into $\tau$ as well or are internally handled to facilitate adaptive splitting or multi-fidelity optimization (e.g., as done in [82]).

All HPO algorithms presented here work by the same principle: they iteratively propose HPCs $\bm{\lambda}^{+}$ and then evaluate the performance on these configurations. We store these HPs and their respective evaluations successively in the so-called archive $\mathcal{A}=((\bm{\lambda}^{(1)},c(\bm{\lambda}^{(1)})),(\bm{\lambda}^{(2)},c(\bm{\lambda}^{(2)})),\dots)$, with $\mathcal{A}^{[t+1]}=\mathcal{A}^{[t]}\cup(\bm{\lambda}^{+},c(\bm{\lambda}^{+}))$ if a single point is proposed by the tuner.

Many algorithms can be characterized by how they handle two different trade-offs: a) The exploration vs. exploitation trade-off refers to how much budget an optimizer spends on either attempting to directly exploit the currently available knowledge base by evaluating very close to the currently best candidates (e.g., local search) or exploring the search space to gather new knowledge (e.g., random search). b) The inference vs. search trade-off refers to how much time and overhead is spent to induce a model from the currently available archive data in order to exploit past evaluations as much as possible. Other relevant aspects that HPO algorithms differ in are: Parallelizability, i.e., how many configurations a tuner can (reasonably) propose at the same time; global vs. local behavior of the optimizer, i.e., if updates are always quite close to already evaluated configurations; noise handling, i.e., if the optimizer takes into account that the estimated generalization error is noisy; multifidelity, i.e., if the tuner uses cheaper evaluations, for example on smaller subsets of the data, to infer performance on the full data; search space complexity, i.e., if and how hierarchical search spaces as introduced in Section 5 can be handled.

As discussed in Section 4.1, HPO produces an approximately optimal HPC $\hat{\bm{\lambda}}=\tau(\mathcal{D},\mathcal{I},\tilde{\Lambda},\rho)$ by optimizing it w.r.t. the resampled performance $c(\bm{\lambda})=\widehat{\mathrm{GE}}(\mathcal{I},\mathcal{J},\rho,\bm{\lambda})$. This is still risk minimization w.r.t. (hyper)parameters, where we search for optimal parameters $\hat{\bm{\lambda}}$ so that the risk of our predictor $\hat{f}_{\hat{\bm{\lambda}},\bm{\hat{\theta}}}$ becomes minimal when measured on validation data via

where $\hat{f}_{\hat{\bm{\lambda}},\bm{\hat{\theta}}}=\mathcal{I}(\mathcal{D}_{\text{train}},\hat{\bm{\lambda}})$ and $\bm{F}_{\hat{\bm{\lambda}},\bm{\hat{\theta}}}$ is the prediction matrix of $\hat{f}$ on validation data, for a pointwise loss function. The above is formulated for a single holdout split $(J_{\mathrm{train}},J_{\mathrm{test}})$ in order to demonstrate the tight connection between (first level) risk minimization and HPO; Eq. $\eqref{eq:ges}$ provides the generalization for arbitrary resampling with multiple folds. This is somewhat obfuscated and complicated by the fact that we cannot evaluate Eq. (8) in one go, but must rather fit one or multiple models $\hat{f}$ during its computation (hence also its black-box nature). It is useful to conceptualize this as a bilevel inference mechanism; while the parameters $\bm{\hat{\theta}}$ of $f$ for a given HPC are estimated in the first level, in the second level we infer the HPs $\hat{\bm{\lambda}}$. However, both levels are conceptually very similar in the sense that we are optimizing a risk function for model parameters which should be optimal for the the data distribution at hand. In case of the second level, this risk function is not $\mathcal{R}_{\text{emp}}(\bm{\theta})$, but the harder-to-evaluate generalization error $\widehat{\mathrm{GE}}$. An intuitive, alternative term for HPO is second level inference [55], visualized in Figure 8.

There are mainly two reasons why such a bilevel optimization is preferable to a direct, joint risk minimization of parameters and HPs [55]:

• •

  Typically, learners are constructed in such a way that optimized first-level parameters can be more efficiently computed for fixed HPs, e.g., often the first-level problem is convex, while the joint problem is not.

• •

  Since the generalization error is eventually optimized for the bilevel approach, the resulting model should be less prone to overfitting.

Thus, we can define a learner with integrated tuning as a mapping $\mathcal{T}_{\mathcal{I},\tilde{\Lambda},\rho,\mathcal{J}}:\mathds{D}\rightarrow\mathcal{H}$, $\mathcal{D}\rightarrow\mathcal{I}_{\tau(\mathcal{D},\mathcal{I},\tilde{\Lambda},\rho)}(\mathcal{D})$, which maps a data set $\mathcal{D}$ to the model $\hat{f}_{\hat{\bm{\lambda}}}$ that has the HPC set to $\hat{\bm{\lambda}}$ as optimized by $\tau$ on $\mathcal{D}$ and is then itself trained on the whole of $\mathcal{D}$; all for a given inducer $\mathcal{I}$, performance measure $\rho$, and search space $\tilde{\Lambda}$. Algorithmically, this learner has a 2-step training procedure (see Figure 8), where tuning is performed before the final model fit. This ?self-tuning? learner $\mathcal{T}$ ?shadows? the tuned HPs of its search space $\tilde{\Lambda}$ from the original learner and integrates their configuration into the training procedure ⁸⁸8The self-tuning learner actually also adds new HPs that are the control parameters of the HPO procedure.. If such a learner is cross-validated, we naturally arrive at the concept of nested CV, which is discussed in the following Section 4.4.

As discussed in Section 3.2, the evaluation of any learner should always be performed via resampling on independent test sets to ensure non-biased estimation of its generalization error. This is necessary because evaluating $\hat{f}$ on the data set $\mathcal{D}$ that was used for its construction would lead to an optimistic bias. In the general HPO problem as in Eq. (10), we already minimize this generalization error by resampling:

If we simply report the estimated $c(\hat{\bm{\lambda}})$ value of the returned best HPC, this also creates an optimistically biased estimator of the generalization error, as we have violated the fundamental ?untouched test set? principle by optimizing on the test set(s) instead.

To better understand the necessity of an additional resampling step, we consider the following example in Figure 9, introduced by [21]. Assume a balanced binary classification task and an inducer $\mathcal{I}_{\bm{\lambda}}$ that ignores the data. Hence, $\bm{\lambda}$ has no effect, but rather ?predicts? the class labels in a balanced but random manner. Such a learner always has a true misclassification error of $\mathrm{GE}(\mathcal{I},\bm{\lambda},n_{\mathrm{train}},\rho)=0.5$ (using $\rho_{CE}$ as a metric), and any normal CV-based estimator will provide an approximately correct value as long as our data set is not too small. We now ?tune? this learner, for example, by RS – which is meaningless, as $\bm{\lambda}$ has no effect. The more tuning iterations are performed, the more likely it becomes that some model from our archive will produce partially correct labels simply by random chance, and the (only randomly) ?best? of these is selected by our tuner at the end. The more we tune, the smaller our data set, or the more variance our GE estimator exhibits, the more expressed this optimistic bias will be.

To avoid this bias, we introduce an additional outer resampling loop around this inner HPO-resampling procedure – or as discussed in Section 4.3, we simply regard this as cleanly cross-validating the self-tuned learner $\mathcal{T}_{\mathcal{I},\tilde{\Lambda},\rho,\mathcal{J}}$. This is called nested resampling, which is illustrated in Figure 10.

The procedure works as follows: In the outer loop, an outer model-building or training set is selected, and an outer test set is cleanly set aside. Each proposed HPC $\bm{\lambda}^{+}$ during tuning is evaluated via inner resampling on the outer training set. The best performing HPC $\hat{\bm{\lambda}}$ returned by tuning is then used to fit a final model for the current outer loop on the outer training set, and this model is then cleanly evaluated on the test set. This is repeated for all outer loops, and all outer test performances are aggregated at the end.

Some further comments on this general procedure: (i) Any resampling scheme is possible on the inside and outside, and these schemes can be flexibly combined based on statistical and computational considerations. Nested CV and nested holdout are most common. (ii) Nested holdout is often called the train-validation-test procedure, with the respective terminology for the generated three data sets resulting from the 3-way split. (iii) Many users often wonder which ?optimal? HPC $\hat{\bm{\lambda}}$ they are supposed to report or study if nested CV is performed, with multiple outer loops, and hence multiple outer HPCs $\hat{\bm{\lambda}}$. However, the learned HPs that result from optimizations within CV are considered temporary objects that merely exist in order to estimate $\widehat{\mathrm{GE}}(\mathcal{I},\mathcal{J},\rho,\bm{\lambda})$. The comparison to first-level risk minimization from Section 4.3 is instructive here: The formal goal of nested CV is simply to produce the performance distribution on outer test sets; the $\hat{\bm{\lambda}}$ can be considered as the fitted HPs of the self-tuned learner $\mathcal{T}$. If the parameters of a final model are of interest for further study, the tuner $\mathcal{T}$ should be fitted one final time on the complete data set. This would imply a final tuning run on the complete data set for second-level inference.

Nested resampling ensures unbiased outer evaluation of the HPO process, but, as CV for the first level, it is only a process that is used to estimate performance – it does not directly help in constructing a better model. The biased estimation of performance values is not a problem for the optimization itself, as long as all evaluated HPCs are still ranked correctly. But after a considerably large amount of evaluations, wrong HPCs might be selected due to stochasticity or overfitting to the splits of the inner resampling. This effect has been called either overtuning, meta-overfitting or oversearching [103, 112]. At least parts of this problem seem directly related to the problem of multiple hypothesis testing. However, it has not been analysed very well yet, and unlike regularization for (first level) ERM, not many counter measures are currently known for HPO.

Most classifiers do not directly output class labels, but rather probabilities or real-valued decision scores, although many metrics require predicted class labels. A score is converted to a predicted label by comparing it to a threshold $t$ so that $\hat{y}=[f(\mathbf{x})\geq t]$, where we use the Iverson bracket $[]$ with $\hat{y}=1$ if $f(\mathbf{x})\geq t$ and $\hat{y}=0$ in all other cases. For binary classification, the default thresholds are $t=0.5$ for probabilities and $t=0$ for scores. However, depending on the metric and the classifier, different thresholds can be optimal. Threshold tuning [122] is the practice of optimizing the classification threshold $t$ for a given model to improve performance. Strictly speaking, the threshold constitutes a further HP that must be chosen carefully. However, since the threshold can be varied freely after a model has been built, it does not need to be tuned jointly with the remaining HPs and can be optimized in a separate, subsequent, and cheaper step for each proposed HPC $\bm{\lambda}^{+}$.

After models have been fitted and predictions for all test sets have been obtained when $c(\bm{\lambda}^{+})$ is computed via resampling, the vector of joint test set scores $\tilde{\bm{F}}$ can be compared against the joint vector of test set labels $\mathbf{y}$ to optimize for an optimal threshold $\hat{t}=\operatorname*{arg\,min}_{t}\rho(\tilde{\mathbf{y}},[\tilde{\bm{F}}\geq t])$, where the Iverson bracket is evaluated component-wise and $t\in[0,1]$ for a binary probabilistic classifier and $t\in\mathds{R}$ for a score classifier. Since $t$ is scalar and evaluations are cheap, $\hat{t}$ can be found easily via a line search. This two-step approach ensures that every HPC is coupled with its optimal threshold. $c(\bm{\lambda}^{+})$ is then defined as the optimal performance value for $\bm{\lambda}^{+}$ in combination with $\hat{t}$.

The procedure can be generalized to multi-class classification. Class probabilities or scores $\hat{\pi}_{k}(\mathbf{x}),k=1,\ldots,g$ are divided by threshold weights $w_{k}$, $k=1,\ldots,g$, and the $k$ that yields the maximal $\frac{\hat{\pi}_{k}(\mathbf{x})}{w_{k}}$ is chosen as the predicted class label. The weights $w_{k},k=1,\ldots,g$ are optimized in the same way as $t$ in the binary case. Generally, threshold tuning can be performed jointly with any HPO algorithm. In practice, threshold tuning is implemented as a post-processing step of an ML pipeline (Sections 5.1 & 5.2).

In the following, we will extend the HPO problem for a specific learner towards configuring a full pipeline including preprocessing. A linear ML pipeline is a succession of preprocessing methods followed by a learner at the end, all arranged as nodes in a linear graph. Each node acts in a very similar manner as the learner, and has an associated training and prediction procedure. During training, each node learns its parameters (based on the training data) and sends a potentially transformed version of the training data to its successor. Afterwards, a pipeline can be used to obtain predictions: The nodes operate on new data according to their model parameters. Figure 11 shows a simple example.

Pipelines are important to properly embed the full model building procedure, including preprocessing, into cross-validation, so every aspect of the model is only inferred from the training data. This is necessary to avoid overfitting and biased performance evaluation [21, 63], as it is for basic ML. As each node represents a configurable piece of code, each node can have HPs, and the HPs of the pipeline are simply the joint set of all HPs of its contained nodes. Therefore, we can model the whole pipeline as a single HPO problem with the combined search space $\tilde{\Lambda}=\tilde{\Lambda}_{\text{op},1}\times\cdots\times\tilde{\Lambda}_{\text{op},k}\times\tilde{\Lambda}_{\mathcal{I}}$.

More flexible pipelining, and especially the selection of appropriate nodes in a data-driven manner via HPO, can be achieved by representing our pipeline as a directed acyclic graph. Usually, this implies a single source node that accepts the original data and a single sink node that returns the predictions. Each node represents a preprocessing operation, a learner, a postprocessing operation, or a directive operation that directs how the data is passed to the child node(s).

One instance of such a pipeline is illustrated in Figure 12.

Here, we consider the choice between multiple mutually exclusive preprocessing steps, as well as the choice between different ML algorithms. Such a choice is represented by a branching operator, which can be configured through a categorical parameter and can determine the flow of the data, resulting in multiple ?modeling paths? in our graph.

The HP space induced by a flexible pipeline is potentially more complex. Depending on the setting of the branching HP, different nodes and therefore different HPs are active, resulting in a very hierarchical search space. If we build a graph that includes a sufficiently large selection of preprocessing steps combined with sufficiently many ML models, the result can be flexible enough to work well on a large number of data sets – assuming it is correctly configured in a data-dependent manner. Combining such a graph with an efficient tuner is the key principle of AutoML [47, 100, 104].

A resampling procedure is usually chosen based on two fundamental properties of the available data: (i) the number of observations, i.e., to what degree do we face a small sample size situation, and (ii) whether the i.i.d. assumption for our data sampling process is violated.

For smaller data sets, e.g., $n<500$, repeated CV with a high number of repetitions should be used to reduce the variance while keeping the pessimistic bias small [21]. The larger the data set, the fewer splits are necessary. Consequently, for data sets of ?medium? size with $500\leq n\leq 50000$, usually 5- or 10-fold CV is recommended. Beyond that, simple holdout might be sufficient. Note that even for large $n$, sample size problems can occur, for example, when data is imbalanced. In that case, repeated resampling might still be required to obtain properly accurate performance estimates. Stratified sampling, which ensures that the relative class frequencies for each train/test split are consistent with the original data set, helps in such a case.

One fundamental assumption about our data is that observations are i.i.d., i.e., $\left(\mathbf{x}^{(i)},y^{(i)}\right)\overset{i.i.d}{\sim}\mathds{P}_{xy}.$ This assumption is often violated in practice. A typical example is repeated measurements, where observations occur in ?blocks? of multiple, correlated data, e.g., from different hospitals, cities or persons. In such a scenario, we are usually interested in the ability of the model to generalize to new blocks. We must then perform CV with respect to the blocks, e.g., ?leave one block out?. A related problem occurs if data are collected sequentially over a period of time. In such a setting, we are usually interested in how the model will generalize in the future, and the rolling or expanding window forecast must be used for evaluation [9] instead of regular CV. However, discussing these special cases is out of scope for this work.

In HPO, resampling strategies must be specified for the inner as well as the outer level of nested resampling. The outer level is simply regular ML evaluation, and all comments from above hold. We advise readers to study further material such as [72]. The inner level concerns the evaluation of $c(\bm{\lambda})$ through resampling during tuning. While the same general comments from above apply, in order to reduce runtime, repetitions can also be scaled down. We are not particularly interested in very accurate numerical performance estimates at the inner level, and we must only ensure that HPCs are properly ranked during tuning to achieve correct selection, as discussed in Section 4.4. Hence, it might be appropriate to use a 10-fold CV on the outside to ensure proper generalization error estimation of the tuned learner, but to use only 2 folds or simple holdout on the inside. In general, controlling the number of resampling repetitions on the inside should be considered an aspect of the tuner and should probably be automated away from the user (without taking away flexible control in cases of i.i.d. violations, or other deviations from standard scenarios). However, not many current tuners provide this, although racing is one of the attractive exceptions.

The choice of the performance measure should be guided by the costs that suboptimal predictions by the model and subsequent actions in the real-word context of applying the model would incur. Often, popular but simple measures like accuracy do not meet this requirement. Misclassification of different classes can imply different costs. For example, failing to detect an illness may have a higher associated cost than mistakenly admitting a person to the hospital. There exists a plethora of performance measures that attempt to emphasize different aspects of misclassification with respect to prediction probabilities and class imbalances, c.f. [72]and many listed in Table 1 in Appendix B. For other applications, it might be necessary to design a performance measure from scratch or based on underlying business key performance indicators (KPIs).

While a further discussion of metrics is again out of scope for this article, two pieces of advice are pertinent. First, as HPO is a black-box, no real constraints exist regarding the mathematical properties of $\rho$ (or the associated outer loss). For first-level risk minimization, on the other hand, we usually require differentiability and convexity of $L$. If this is not fulfilled, we must approximate the KPI with a more convenient version. Second, for many applications, it is quite unclear whether there is a single metric that captures all aspects of model quality in a balanced manner. In such cases, it can be preferable to optimize multiple measures simultaneously, resulting in a multi-criteria optimization problem [62].

For HPO, it is necessary to define the search space $\tilde{\Lambda}$ over which the optimization is to be performed. This choice can have a large impact on the performance of the tuned model. A simple search space is a (lower dimensional) Cartesian product of individual HP sets that are either numeric (continuous or integer-valued) or categorical. Encoding categorical values as integers is a common mistake that degrades the performance of optimizers that rely on information about distances between HPCs, such as BO. The search intervals of numeric HPs typically must be bounded within a region of plausibly well-performing values for the given method and data set.

Many numeric HPs are often either bounded in a closed interval (e.g., $[0,1]$) or bounded from below (e.g., $[0,\infty)$). The former can usually be tuned without modifications. HPs bounded by a left-closed interval should often be tuned on a logarithmic scale with a generous upper bound, as the influence of larger values often diminishes. For example, the decision whether $k$-NN uses $k=2$ vs. 3 neighbors will have a larger impact than whether it uses $k=102$ vs. $k=103$ neighbors. The logarithmic scale can either be defined in the tuning software or must be set up manually by adjusting the algorithm to use transformations: If the desired range of the HP is $[a,b]$, the tuner optimizes on $[\log a,\log b]$, and any proposed value is transformed through an exponential function before it is passed to the ML algorithm. The logarithm and exponentiation must refer to the same base here, but which base is chosen does not influence the tuning process.

The size of the search space will also considerably influence the quality of the resulting model and the necessary budget of HPO. If chosen too small, the search space may not contain a particularly well-performing HPC. Choosing too wide HP intervals or including inadequate HPs in the search space can have an adverse effect on tuning outcomes in the given budget. If $\tilde{\Lambda}$ is simply too large, it is more difficult for the optimizer to find the optimum or promising regions within the given budget. Furthermore, restricting the bounds of an HP may be beneficial to avoid values that are a priori known to cause problems due to unstable behavior or large resource consumption. If multiple HPs lead to poor performance throughout a large part of their range – for example, by resulting in a degenerate ML model or a software crash – the fraction of the search space that leads to fair performance then shrinks exponentially in the number of HPs with this problem. This effect can be viewed as a further manifestation of the so-called curse of dimensionality.

Due to this curse of dimensionality and the considerable runtime costs of HPO, we would like to tune as few HPs as possible. If no prior knowledge from earlier experiments or expert knowledge exists, it is common practice to leave other HPs at their software default values with the assumption that the developers of the algorithm chose values that work well under a wide range of conditions, which is not necessarily given and its often not documented how these defaults were specified. Recent approaches have studied how to empirically find optimal default values, tuning ranges and HPC prior distributions based on extensive meta-data [142, 141, 47, 109, 138, 251, 108, 53].

It is possible to optimize several learning algorithms in combination, as described in Section 5.2, but this introduces HP dependencies. The question then arises of which of the large number of ML algorithms (or preprocessing operations) should be considered. However, [193] showed that in many cases, only a small but diverse set of learners is necessary to choose one ML algorithm that performs sufficiently well.

The number of HPs considered for optimization has a large influence on the difficulty of an HPO problem. Particularly large search spaces typically arise from optimizing over large pipelines or multiple learners (Section 5). With very few HPs (up to about 2-3) and well-functioning discretization, GS may be useful due to its interpretable, deterministic, and reproducible nature; however, it is not recommended beyond this [11]. BO with GPs work well for up to around 10 HPs. However, more HPs typically require more function evaluations – which in turn is problematic runtime-wise, since GPs scale cubically with the number of dimensions. On the other hand, BO with RFs have been used successfully on search spaces with hundreds of HPs [135] and can usually handle mixed hierarchical search spaces better. Pure sampling-based methods, such as RS and Hyperband, work well even for very large HP spaces as long as the ?effective? dimension (i.e., the number of HPs that have a large impact on performance) is low, which is often observed in ML models [11]. Evolutionary algorithms (and those using similar metaheuristics, such as racing) can also work with truly large search spaces, and even with search spaces of arbitrarily complex structure if one is willing do use non-custom mutation and crossover operators. Evolutionary algorithms may also require fewer objective evaluations than RS. Therefore, they occupy a middle ground between (highly complex, sample-efficient) BO and (very simple but possibly wasteful) RS.

Another property of algorithms that is especially relevant to practitioners with access to large computational resources is parallelizability, which is discussed in Subsection 6.7.2. Furthermore, HPO algorithms differ in their simplicity, both in terms of algorithmic principle and in terms of usability of available implementations, which can often have implications for their usefulness in practice. While more complex optimization methods, such as those based on BO, are often more sample-efficient or have other desirable properties compared to simple methods, they also have more components that can fail. When performance evaluations are cheap and the search space is small, it may therefore be beneficial to fall back on simple and robust approaches such as RS, Hyperband, or any tuner with minimal inference overhead. The availability of any implementation at all (and the quality of that implementation) is also important; there may be optimization algorithms based on beautiful theoretical principles that have a poorly maintained implementation or are based on out-of-date software. The additional cost of having to port an algorithm to the software platform being used, or even implement it from scratch, could be spent on running more evaluations with an existing algorithm.

One might wish to select an HPO algorithm that performed best on previous benchmarks. However, no single benchmark exists which includes all relevant scenarios and whose results generalize to all possible applications. Specific benchmark results can therefore only indicate how well an algorithm works for a selected set of data sets, a predefined budget, specific parallelization, specific learners, and search spaces. Even worse, extensive comparison studies are missing in the current literature, although efforts have been made to establish unified benchmarks. [35] showed that (i) BO with GPs were a strong approach for small continuous spaces with few evaluations, (ii) BO with Random Forests performed well for mixed spaces with a few hundred evaluations, and (iii) for large spaces and many evaluations, ES were the best optimizers.

In addition to the choice of algorithm, the choice of implementation is also important in practice. We note that a thorough comparison of all available software packages and frameworks is not within the scope of this article. Nevertheless, we provide a brief list of popular implementations and frameworks as reasonable starting points for newcomers.

Simple and established tuning algorithms for HPO are usually shipped with general purpose ML frameworks such as Scikit-learn [248], Tensorflow⁹⁹9https://www.tensorflow.org/, PyTorch [247], MXNet [176], mlr3 [226, 163], tidymodels¹⁰¹⁰10https://www.tidymodels.org/ or h2o.ai¹¹¹¹11https://github.com/h2oai/h2o-3. Although modern state-of-the-art algorithms often build on, extend, or connect to such an ML framework, they are usually developed in independent software projects.

For Python, there exist a plethora of HPO toolkits, e.g., Spearmint [262], SMAC¹²¹²12https://github.com/automl/SMAC3 [213], BoTorch [158], Dragonfly [77], or Oríon¹³¹³13https://github.com/Epistimio/orion. Multiple HPO methods are supported by toolkits like Hyperopt [160], Optuna [153] or Weights & Biases¹⁴¹⁴14https://www.wandb.com/. A popular framework that combines modern HPO approaches with the Scikit-learn toolbox for ML in Python is auto-sklearn [194].

The availability of HPO implementations in R is comparably smaller. However, more established HPO tuners are either already shipped with the ML framework, or can be found, e.g., in the packages mlrMBO¹⁵¹⁵15https://github.com/mlr-org/mlrMBO, irace [234], DiceOptim [254], or rBayesianOptimization¹⁶¹⁶16https://github.com/yanyachen/rBayesianOptimization.

See Appendix C for more information

Choosing an appropriate budget for HPO or dynamically terminating HPO based on collected archive data is a difficult practical challenge that is not readily solved. The user typically has these options: (i) A certain amount of runtime is specified before HPO is started, solely based on practical considerations and intuition. This is what currently happens nearly exclusively in practice. A simple rule-of-thumb might be scaling the budget of HPC evaluations with search space dimensionality in a linear fashion like $50\times l$ or $100\times l$, which would also define the number of full budget units in a multi-fidelity setup. The downside is that while it is usually simpler to define when a result is needed, it is nearly impossible to determine whether this computational budget is enough to solve the HPO problem reasonably well. (ii) The user can set a lower bound regarding the required generalization error, i.e., what model quality is deemed good enough in order to put the resulting model to practical use. Notably, if this lower bound is set too optimistically, the HPO process might never terminate or simply take too long. (iii) If no (significant) progress is observed for a certain amount of time during tuning, HPO is considered to have converged or stagnated and stopped. This criterion bears the risk of terminating too early, especially for large and complex search spaces. Furthermore, this approach is conceptually problematic, as it mainly makes sense for iterative, local optimizers. No HPO algorithm from this article, except for maybe ES, belongs to this category, as they all contain a strong exploratory, global search component. This often implies that search performance can flatline for a longer time and then abruptly change. (iv) For BO, a small maximum value of the $\operatorname{EI}$ acquisition function (see Eq. 11) indicates that further progress is estimated to be unlikely [76]. This criterion could even be combined with other, no-BO type tuners. It implies that the surrogate model is trustworthy enough for such an estimate, including its local uncertainty estimate. This is somewhat unlikely for large, complex search spaces. In practice, it is probably prudent to set up a combination of all mentioned criteria (if possible) as an ?and? or ?or? combination, with somewhat liberal thresholds to continue the optimization for (ii), (iii), (iv), and (v), and an absolute bound on the maximal runtime for (i).

With many HPO tuners, it is possible to continue the optimization even after it has been terminated. One can even use (part of) the archive as an initialization of the next HPO run, e.g., as an initial design for BO or as the initial population for an ES.

Some methods, e.g., ES and under some circumstances BO, may get stuck in a subspace of the search space and fail to explore other regions. While some means of mitigation exist, such as interleaving proposed points with randomly generated points, there is always the possibility of re-starting the optimization from scratch multiple times, and selecting the best performance from the combined runs as the final result. The decision regarding the termination criterion is itself always a matter of cost-benefit calculation of potential increased performance against cost of additional evaluations. One should also consider the possibility that, for more local search strategies like ES, terminating early and restarting the optimization can be more cost-effective than letting a single optimization run continue for much longer. For more global samplers like BO and HB, it is much less clear whether and how such an efficient restart could be executed.

HPO may require substantial computational resources, as a large optimization search space may require many performance evaluations, and individual performance evaluations can also be very computationally expensive. One way to reduce the computational time are warm-starts, where information from previous experiments are used as a starting solution.

HP tuning is usually computationally expensive, and running optimization in parallel is one way to reduce the overall time needed for a satisfactory result. During parallelization, computational jobs (which can be arbitrary function calls) are mapped to a pool of workers (usually distributed among individual physical CPU cores, or remote computational nodes). In practice, a parallelization backend orchestrates starting up, shutting down, and communication with the workers. The achieved speedup can be proportional to the number of workers in ideal cases, but less than linear scaling is typically observed, depending on the algorithm used (Amdahl’s law, [116]). Besides parallelization, ensuring a fail-safe program flow is often also mandatory, and requires special care.

In practical model selection, it is often the case that practitioners care not only about predictive performance, but also about other aspects of the model such as its simplicity and interpretability. If a practitioner considers only a low number of different model classes, then a practical alternative to AutoML on a full pipeline space of different ML models is to construct a manual sequence of preferred models of increasing complexity, to separately tune them, and to compare their performance. An informed choice can then be made about how much additional complexity is acceptable for a certain amount of improved predictive performance. [57] suggest the ?one-standard-error rule? and use the simplest model within one standard error of the performance of the best model. The specific sequence of learners will differ for each use case but might look something like this: (1) linear model with only a single feature or tree stump, (2) L1/L2 regularized linear model or a CART tree, (3) component-wise boosting of a GAM, (4) random forest, (5) gradient tree boosting, and (6) full AutoML search space with pipelines. If even more focus is put on predictive performance (multi-level), stacking a diverse set of models can give a final (often rather small) boost [143].

While few well-designed large-scale benchmarks of HPO exist [81, 36, 35], there is still an ongoing endeavour of the community for better and more comprehensive studies. Usually, tuners are compared on fixed sets of data sets and search spaces or cheaper approximations via empirical performance models / surrogates [36]. It is still unclear how an ideal and representative benchmark set for HPO that also enables efficient benchmarking might look like.