Linear models for binomial outcomes

Linear models for binomial outcomes derive from the stats package which is part of the R core distribution (R Core Team 2019). Hyperparameters for these models include linkage functions and are shown in the table below.

Linear models for multinomial outcomes

Prior to version 1.3.0 the linear model for multinomial outcomes was implemented using the VGAM::vglm function (T. W. Yee and Wild 1996; T. Yee 2010; Thomas W. Yee 2015). For better scalability, this has been replaced by the nnet::multinom function (Venables and Ripley 2002). Hyperparameters for this model are shown in the table below.

Linear models for continuous and count-type outcomes

Linear models for continuous and count-type outcomes derive from the stats package of the R core distribution (R Core Team 2019). Hyperparameters for these models include linkage functions and are shown in the table below.

Linear models for survival outcomes

Linear models for survival outcomes are divided into semi-parametric and parametric models. The semi-parametric Cox proportional hazards model (Cox 1972) is based on the survival::coxph function (Therneau and Grambsch 2000). Tied survival times are resolved using the default method by Efron (1977).

Various fully parametric models are also available, which differ in the assumed distribution of the outcome, i.e. the linking function. The parametric models are all based on survival::survreg (Therneau and Grambsch 2000).

Hyperparameters for semi-parametric and parametric survival models are shown in the table below.

Predicting run time of model

Models take a certain time to train. Familiar is actively measuring this time during hyperparameter optimisation for two reasons: first, to optimise assignment of jobs to parallel nodes; and secondly, to prune potential challenger sets that produce models that tend to run longer than the best-known model. The strictness increases with increased exploration.

Let \(m\) be the number of bootstraps used to quantify the most visited hyperparameter sets, and \(M\) the total number of bootstrap samples that can be visited. Then let \(t_{\texttt{opt}}\) be the runtime of the best known model. The maximum time that a challenger hyperparameter set is allowed for training is then empirically set to:

\[t_{\texttt{max}} = \left(5 - 4 \frac{m}{M} \right) t_{\texttt{opt}} \]

Hence, \(t_{\texttt{max}}\) converges to \(t_{\texttt{opt}}\) for \(m \rightarrow M\).

If maximum runtime is relatively insignificant, i.e. \(t_{\texttt{max}} < 10.0\) seconds, a threshold of \(t_{\texttt{max}}= 10.0\) is used.

If the maximum runtime is not known, i.e. none of the hyperparameter sets evaluated so far produced a valid model, the maximum time threshold is set to infinite. In effect, no hyperparameter sets are pruned based on expected run time.

A random forest is trained to infer runtimes for (new) hyperparameter sets, based on the runtimes observed for visited hyperparameter sets. The random forest subsequently infers runtime for a challenger hyperparameter set. The runtime estimate is compared against \(t_{\texttt{max}}\), and if it exceeds the threshold, it is rejected and not evaluated.

Assessing goodness of hyperparameter sets

A summary score \(s\) determines how good a set of hyperparameters \(\mathbf{x}\) is. This score is computed in two steps. First an objective score is computed to assess model performance for a hyperparameter set in a bootstrap subsample. Subsequently an optimisation score is computed within the subsample. Finally a summary score is determined from the optimisation score over all subsamples.

An objective score \(s''\) is computed from the performance metric value for a specific hyperparameter set for in-bag (IB; \(s''_{\textrm{IB}}\)) and out-of-bag data (OOB; \(s''_{\textrm{OOB}}\)). The objective score always lies in the interval \([-1.0, 1.0]\). An objective score of \(1.0\) always indicates the best possible score. A score of \(0.0\) indicates that the hyperparameter set leads to the same degree of performance as the best educated guess, i.e. the majority class (for binomial and multinomial outcomes), the median value (for continuous and count outcomes), or tied risk or survival times (for survival outcomes). Objective scores that are either missing or below \(-1.0\) are truncated to the value of \(-1.0\). In case multiple metrics are used to assess model performance, the mean of the respective objective scores is used.

The optimisation scores \(s'\) and summary score \(s\) of each hyperparameter set are subsequently computed using one of the following functions, which can be set using the optimisation_function parameter:

1. validation, max_validation (default): The optimisation score for each bootstrap subsample is \(s'=s''_{\textrm{OOB}}\). The summary score \(s\) is then the average of optimisation scores \(s'\). This is a commonly used criterion that tries to maximise the score on the OOB validation data.

2. balanced: The optimisation score for each bootstrap subsample is \(s'=s''_{\textrm{OOB}} - \left| s''_{\textrm{OOB}} - s''_{\textrm{IB}} \right|\). The summary score \(s\) is then the average of optimisation scores \(s'\). A variation on max_validation with a penalty for differences in performance between the IB and OOB data. The underlying idea is that a good set of hyperparameters should lead to models that perform well on both development and validation data.

3. strong_balance: The optimisation score for each bootstrap subsample is \(s'=s''_{\textrm{OOB}} - 2\left| s''_{\textrm{OOB}}-s''_{\textrm{IB}} \right|\). The summary score \(s\) is then the average of optimisation scores \(s'\). A variant of balanced with a stronger penalty term.

4. validation_minus_sd: The optimisation score for each bootstrap subsample is \(s'=s''_{\textrm{OOB}}\). The summary score \(s\) is then the average of optimisation scores \(s'\) minus the standard deviation of \(s'\). This penalises hyperparameter sets that lead to wide variance on OOB data.

5. validation_25th_percentile: The optimisation score for each bootstrap subsample is \(s'=s''_{\textrm{OOB}}\). The summary score \(s\) is the 25th percentile of optimisation scores \(s'\). Like validation_minus_sd this penalises hyperparameter sets that lead to wide variance on OOB data.

6. model_estimate: The optimisation score for each bootstrap subsample is \(s'=s''_{\textrm{OOB}}\). The model used to infer utility of new hyperparameter sets (see the next section Predicting optimisation score for new hyperparameter sets) is used to predict the expected optimisation score \(\mu(\mathbf{x})\). This score is used as summary score \(s\). Not available for random search.

7. model_estimate_minus_sd: The optimisation score for each bootstrap subsample is \(s'=s''_{\textrm{OOB}}\). The model used to infer utility of new hyperparameter sets (see the next section Predicting optimisation score for new hyperparameter sets) is used to predict the expected optimisation score \(\mu(\mathbf{x})\) and its standard deviation \(\sigma(\mathbf{x})\). Then \(s = \mu(\mathbf{x}) - \sigma(\mathbf{x})\). This penalises hyperparameter sets for which a large variance is expected. Not available for random search.

8. model_balanced_estimate: The optimisation score for each bootstrap subsample is \(s'=s''_{\textrm{OOB}} - \left| s''_{\textrm{OOB}} - s''_{\textrm{IB}} \right|\). The model used to infer utility of new hyperparameter sets (see the next section Predicting optimisation score for new hyperparameter sets) is used to predict the expected optimisation score \(\mu(\mathbf{x})\). This score is used as summary score \(s\). The method is similar to model_estimate, but predicts the balanced score instead of the validation score. Not available for random search.

9. model_balanced_estimate_minus_sd: The optimisation score for each bootstrap subsample is \(s'=s''_{\textrm{OOB}} - \left| s''_{\textrm{OOB}} - s''_{\textrm{IB}} \right|\). The model used to infer utility of new hyperparameter sets (see the next section Predicting optimisation score for new hyperparameter sets) is used to predict the expected optimisation score \(\mu(\mathbf{x})\) and its standard deviation \(\sigma(\mathbf{x})\). Then \(s = \mu(\mathbf{x}) - \sigma(\mathbf{x})\). This penalises hyperparameter sets for which a large variance is expected. The method is similar to model_estimate_minus_sd, but predicts the balanced score and its variance instead. Not available for random search.

The summary score is then used to select the best known hyperparameter set, i.e. the set that maximises the summary score. Moreover, the optimisation scores are used to identify candidate hyperparameter sets as described in the following section. The optimisation function is set using the optimisation_function argument.

Predicting optimisation score for new hyperparameter sets

Any point in the hyperparameter space has a single, scalar, optimisation score value that is a priori unknown. During the optimisation process, the algorithm samples from the hyperparameter space by selecting hyperparameter sets and computing the optimisation score for one or more bootstraps. For each hyperparameter set the resulting scores are distributed around the actual value.

A key point of Bayesian optimisation is the ability to estimate the usefulness or utility of new hyperparameter sets. This ability is facilitated by modelling the optimisation score of new hyperparameter sets using the optimisation scores of observed hyperparameter sets.

The following models can be used for this purpose:

• gaussian_process (default): Creates a localised approximate Gaussian deterministic Gaussian Process. This is implemented using the laGP package (Gramacy 2016).

• bayesian_additive_regression_trees or bart: Uses Bayesian Additive Regression Trees for inference. Unlike standard random forests, BART allows for estimating posterior distributions directly and can extrapolate. BART is implemented using the BART package (Sparapani, Spanbauer, and McCulloch 2021).

• random_forest: Creates a random forest for inference. A random forest was originally used in the SMBO algorithm by Hutter et al. (Hutter, Hoos, and Leyton-Brown 2011). A weakness of random forests is their lack of extrapolation beyond observed values, which limits their usefulness in exploiting promising areas of hyperparameter space somewhat.

In addition, familiar can perform random search (random or random_search). This forgoes the use of models to steer optimisation. Instead, the hyperparameter space is sampled at random.

The learner used to predict optimisation scores can be specified using the hyperparameter_learner parameter.

Acquisition functions for utility of hyperparameter sets

The expected values and the posterior values of optimisation scores from the models are used to compute utility of new hyperparameter sets using an acquisition function. The following acquisition functions are available in familiar (Shahriari et al. 2016). Let \(\nu=f(\mathbf{x})\) be the posterior distribution of hyperparameter set \(\mathbf{x}\), and \(\tau\) the best observed optimisation score. Let also \(\mu(\mathbf{x})\) and \(\sigma(\mathbf{x})\) be the sample mean and sample standard deviation for set \(\mathbf{x}\) (implicitly for round \(m\)):

• improvement_probability: The probability of improvement quantifies the probability that the expected optimisation score for a set \(\mathbf{x}\) is better than the best observed optimisation score \(\tau\):

  \[\alpha(\mathbf{x}) = P(\nu > \tau) = \Phi \left(\frac{\mu(\mathbf{x}) - \tau} {\sigma(\mathbf{x})}\right)\]

  Here \(\Phi\) is the cumulative distribution function of a normal distribution. Note that this acquisition function is typically prone to convergence to local minima, although in our algorithm this may be mitigated by the alternating random search strategy.

• improvement_empirical_probability: Similar to improvement_probability, but based directly on optimisation scores predicted by the individual decision trees:

  \[\alpha(\mathbf{x}) = P(\nu > \tau) = \frac{1}{N} \sum_{i=1}^N \left[\nu_i > \tau\right]\]

  with \(\left[\ldots\right]\) denoting an Iverson bracket, which takes the value 1 if the condition specified by \(\ldots\) is true, and 0 otherwise.

• expected_improvement (default): Expected improvement is based on an improvement function, and is be computed as follows:

  \[\alpha(\mathbf{x})=(\mu(\mathbf{x})-\tau) \Phi \left(\frac{\mu(\mathbf{x}) - \tau}{\sigma(\mathbf{x})}\right) + \sigma(\mathbf{x}) \phi\left(\frac{\mu(\mathbf{x}) - \tau}{\sigma(\mathbf{x})}\right)\]

  with \(\phi\) denoting the normal probability density function. If \(\sigma(\mathbf{x}) = 0\), \(\alpha(\mathbf{x})=0\).

• upper_confidence_bound: This acquisition function is based on the upper confidence bound of the distribution, and is defined as (Srinivas et al. 2012):

  \[\alpha(\mathbf{x}, t)=\mu(\mathbf{x}) + \beta_t \sigma(\mathbf{x})\]

  Here the parameter \(\beta_t\) depends on the current round \(t\), which in our implementation would be roughly equivalent to the maximum number of bootstraps performed for any of the hyperparameter sets, minus \(1\). The \(\beta_t\) parameter is then defined as \(\beta_t= 2 \log(|D|(t+1)^2\pi^2/6\delta)\) (Srinivas et al. 2012). Here, \(\delta \in [0,1]\) and \(|D|\) the dimensionality of the hyperparameter space. For their experiments Srinivas et al. used \(\delta=0.1\), and noted that scaling down \(\beta_t\) by a factor \(5\) improves the algorithm. Hence, we use \(\beta_t= \frac{2}{5} \log\left(\frac{5}{3} |D|(t+1)^2\pi^2\right)\).

• bayes_upper_confidence_bound: Proposed by Kaufmann et al. (Kaufmann, Cappé, and Garivier 2012), the Bayesian upper confidence bound is defined as:

  \[\alpha(\mathbf{x}, t)=Q \left(1-\frac{1}{(t + 1) \log(n)^c}, \nu_t\right)\]

  with \(Q\) being a quantile function. Note that we here deviate from the original definition by substituting the original \(t\) by \(t+1\), and not constraining the posterior distribution \(\nu\) to a known function. In the original, \(n\) defines the number of rounds, with \(t\) the current round. In our algorithm, this is interpreted as \(n\) being the total number of bootstraps, and \(t\) the maximum number of bootstraps already performed for any of the hyperparameter sets. Hence our implementation is the same, aside from the lack of constraint on the posterior distribution. Based on the simulations and recommendations by Kaufmann et al. we choose parameter \(c=0\).

The acquisition function can be specified using the acquisition_function parameter.

Exploring challenger sets

After selecting challenger hyperparameter sets, the optimisation score of models trained using these hyperparameter sets are compared against the best-known score in a run-off intensify step. During this step, models will be trained on new bootstrap data, compared, and then trained on another set of bootstraps, and so on. Pruning the set of challenger hyperparameters after each iteration reduces computational load. Familiar implements the following methods:

• single_shot (default): The set of alternative parameter sets is not pruned, and each intensify iteration contains only a single (initial) intensification step that only uses a single bootstrap. This is the exploration method that under default settings does the least in-depth exploration, with the advantage of being fast. Extensive internal tests show that this exploration method yields similarly performant models as more extensive exploration methods.

• successive_halving: The set of alternative parameter sets is pruned by removing the worst performing half of the sets after each step (Jamieson and Talwalkar 2016). The set of investigated parameter sets gets progressively smaller.

• stochastic_reject: The set of alternative parameter sets is pruned by comparing the performance of each parameter set with that of the incumbent best parameter set using a paired Wilcoxon test.

• none: The set of alternative parameter sets is not pruned. Compared to single_shot, hyperparameter sets are assessed using multiple bootstraps (smbo_step_bootstraps * smbo_intensify_steps).

The method used to steer exploration can be set using the exploration_method parameter.

Providing hyperparameters manually

It is possible to set hyperparameters manually. This can be used to change hyperparameters that are fixed by default, to set a fixed value for randomised hyperparameters, or to provide a different search range for randomised hyperparameters.

Hyperparameters can be provided using the hyperparameter tag. For the glm_logistic learner an example tag may for example look as follows:

<hyperparameter>
  <glm_logistic>
    <sign_size>5</sign_size>
  </glm_logistic>
</hyperparameter>

Or as a nested list passed as the hyperparameter argument to summon_familiar:

hyperparameter = list("glm_logistic"=list("sign_size"=5))

More than one value can be provided. The behaviour changes depending on whether the hyperparameter is categorical or numeric variable. In case of categorical hyperparameters, the provided values define the search range. For numerical hyperparameters, providing two values sets the bounds of the search range, whereas providing more than two values will define the search range itself.

Configuration options for hyperparameter optimisation

Hyperparameter optimization may be configured using the tags/arguments in the table below.