#  HEBO: Pushing The Limits of Sample-Efficient Hyperparameter Optimisation

Alexander I. Cowen-Rivers<EMAIL_ADDRESS>Wenlong
Lyu11footnotemark: 1 22footnotemark: 2<EMAIL_ADDRESS>Rasul
Tutunov11footnotemark: 1 22footnotemark: 2<EMAIL_ADDRESS>
Zhi Wang 22footnotemark: 2<EMAIL_ADDRESS>Antoine Grosnit
22footnotemark: 2<EMAIL_ADDRESS>Ryan Rhys Griffiths
22footnotemark: 2<EMAIL_ADDRESS>Alexandre Max Maravel
22footnotemark: 2<EMAIL_ADDRESS>Hao Jianye 22footnotemark: 2
<EMAIL_ADDRESS>Jun Wang 22footnotemark: 2<EMAIL_ADDRESS>Jan Peters
33footnotemark: 3<EMAIL_ADDRESS>Haitham Bou-Ammar 22footnotemark:
2 66footnotemark: 6<EMAIL_ADDRESS>Equal contribution Huawei Noah’s
Ark Lab Technische Universität Darmstadt Corresponding Author University of
Cambridge University College London Honorary position.

###### Abstract

In this work we rigorously analyse assumptions inherent to black-box
optimisation hyper-parameter tuning tasks. Our results on the Bayesmark
benchmark indicate that heteroscedasticity and non-stationarity pose
significant challenges for black-box optimisers. Based on these findings, we
propose a Heteroscedastic and Evolutionary Bayesian Optimisation solver
(HEBO). HEBO performs non-linear input and output warping, admits exact
marginal log-likelihood optimisation and is robust to the values of learned
parameters. We demonstrate HEBO’s empirical efficacy on the NeurIPS 2020
Black-Box Optimisation challenge, where HEBO placed first. Upon further
analysis, we observe that HEBO significantly outperforms existing black-box
optimisers on 108 machine learning hyperparameter tuning tasks comprising the
Bayesmark benchmark. Our findings indicate that the majority of hyper-
parameter tuning tasks exhibit heteroscedasticity and non-stationarity, multi-
objective acquisition ensembles with Pareto front solutions improve queried
configurations, and robust acquisition maximisers afford empirical advantages
relative to their non-robust counterparts. We hope these findings may serve as
guiding principles for practitioners of Bayesian optimisation. All code is
made available at https://github.com/huawei-noah/HEBO.

##  1 Introduction

Although achieving significant success across numerous applications (?, ?, ?,
?, ?), the performance of machine learning models chiefly depends on the
correct setting of hyper-parameters. As models grow larger and more complex,
efficient and autonomous hyper-parameter tuning algorithms become crucial
determinants of performance. A variety of methods from black-box and multi-
fidelity optimisation (?, ?) have been adopted for hyperparameter tuning with
varying degrees of success. Techniques such as Bayesian optimisation (BO), for
example, enable sample efficiency (in terms of black-box evaluations) at the
expense of high computational demands, while “unguided” bandit-based
approaches can fail to converge (?). Identifying such failure modes, the
authors in (?) built on (?) and proposed a combination of bandits and BO that
achieves the best of both worlds; fast convergence and computational
scalability. More recently in the context of the 2020 NeurIPS competition on
Black-Box Optimisation, many BO variants have been convincingly demonstrated
to be superior to random search for the task of hyper-parameter tuning (?).
Though impressive, such successes of BO and alternative black-box optimisers,
belie a set of restrictive modelling and acquisition function assumptions. We
begin by describing these assumptions.

Modelling Assumptions: A core determinant of BO performance is the set of data
modelling assumptions required to specify an appropriate probabilistic model
of the black-box objective (e.g., the choice of validation loss in hyper-
parameter tuning tasks). The model should not only provide accurate point
estimates, but should also maintain calibrated uncertainty estimates to guide
exploration of the objective. Amongst many possible surrogates (?, ?),
Gaussian processes (?) (GPs) are the default choice due to their flexibility
and sample efficiency. Growing interest in applications of Bayesian
optimisation has catalysed engineering feats that enhance scalability and
training efficiency of GP surrogates by exploiting graphical processing units
(?, ?).

Similar to any other framework, the correct specification of a GP model is
dictated by the data modelling assumptions imposed by the user. For instance,
a homoscedastic GP suffers from misspecification when required to model data
with heteroscedastic noise whilst stationary GPs fail to track non-stationary
targets. The aforementioned shortcomings are not unnatural across a range of
real-world problems (?, ?, ?) and hyper-parameter tuning of machine learning
algorithms is no exception, as illustrated in our hypothesis tests of Section
3.2. Hence, even if one succeeds in improving computational efficiency,
frequently-made assumptions such as homoscedasticity and stationarity can
easily inhibit the performance of any BO-based hyper-parameter tuning
algorithm. Despite the importance of these assumptions in practice, GPs that
presume homoscedasticity and stationarity still constitute the most common
choice of surrogate.

Acquisition Function & Optimiser Assumptions: Modelling choices such as those
described above are not unique to the GP fitting procedure but rather
transcend to other steps in the BO algorithm. Precisely, given a model that
adheres to some (or all) assumptions mentioned above, the second step involves
maximising an acquisition function to query novel input locations that are
then evaluated. Hence, practitioners introduce additional constraints relating
to the category of optimisation variables and the choice of acquisition
function. When it comes to variable categories, mainstream implementations (?,
?) assume continuous domains and employ first and second-order optimisers such
as LBFGS (?) and ADAM (?) to propose query locations. Real-valued
configurations cover but a subset of possible machine learning hyper-
parameters rendering discrete variable categories out of scope, an example
being the hidden layer size in deep networks. Moreover, from the point of view
of acquisition functions, libraries tend to presuppose that one unique
acquisition performs best in a given task, while research has shown that
benefits that can arise from a combined solution (?, ?, ?) as we demonstrate
in Section 5.

Contributions: Having identified important modelling choices in BO, our goal
in this paper is to provide empirical insight into the impact of modelling
choice on empirical performance. As a case study, we consider best practices
for hyper-parameter tuning. We wish for our findings to be applicable across a
broad range of tasks and datasets, be attentive to the effect of random
initialisation on algorithmic performance, and naturally, be reproducible. As
such, we prefer to build on established benchmark packages, especially those
that facilitate fast and scalable evaluations with multi-seeding protocols. To
that end, we undertake our evaluation in 2140 experiments from 108 real-world
problems from the UCI repository (?), which was also the testbed of choice for
the NeurIPS 2020 Black-Box Optimisation challenge (?). Our findings point
towards the following conclusions:

  1. 1.

Hyper-parameter tuning tasks exhibit significant levels of heteroscedasticity
and non-stationarity.

  2. 2.

Input and output warping mitigate the effects of heteroscedasticity and non-
stationarity giving rise to better performing tuning algorithms with higher
mean and median performance across all 108 black-box functions under
examination.

  3. 3.

Individual acquisition functions tend to conflict in their solution (i.e., an
optimum for one acquisition function can be a sub-optimal point for another
and vice versa). Using a multi-objective formulation significantly improves
performance;.

To verify our principal conclusions, we conduct additional ablation studies on
our proposed solution method, Heteroscedastic and Evolutionary Bayesian
Optimisation (HEBO) which attempts to address the shortcomings identified in
our analysis and placed first in the 2020 NeurIPS Black-Box Optimisation
Challenge. We obtain a ranked order of importance for significant components
of HEBO, finding that output warping, multi-objective acquisitions and input
warping lead to the most significant improvements followed by robust
acquisition function formulations.

##  2 Standard Design Choices in BO

As discussed earlier, the problem of hyper-parameter tuning can be framed as
an instance of black-box optimisation:

$\arg\max_{\bm{x}\in\mathcal{X}}f(\bm{x}),$ (1)

with $\bm{x}$ denoting a configuration choice, $\mathcal{X}$ a (potentially)
mixed design space, and $f(\bm{x})$ a validation accuracy we wish to maximise.
In this paper, we focus on BO as a solution concept for black-box problems of
the form depicted in Equation 1. BO considers a sequential decision approach
to the global optimisation of a black-box function
$f:\mathcal{X}\rightarrow\mathbb{R}$ over a bounded input domain
$\mathcal{X}$. At each decision round, $i$, the algorithm selects a collection
of $q$ inputs $\bm{x}^{(\text{new})}_{1:q}\in\mathcal{X}^{q}$ and observes
values of the _black-box_ function
$\bm{y}^{(\text{new})}_{1:q}=f(\bm{x}^{(\text{new})}_{1:q})$. The goal is to
rapidly approach the maximum
$\bm{x}^{\star}=\arg\max_{\bm{x}\in\mathcal{X}}f(\bm{x})$. Since both
$f(\cdot)$ and $\bm{x}^{\star}$ are unknown, solvers need to trade off
exploitation and exploration during this search process.

To achieve this goal, BO algorithms operate in two steps. In the first, a
Bayesian model is learned, while in the second an acquisition function
determining new query locations is maximised. Next, we survey frequently-made
assumptions in mainstream BO implementations and contemplate their
implications for performance.

###  2.1 Modelling Assumptions

When black-boxes are real-valued, Gaussian processes (?) are effective
surrogates due to their flexibility and ability to maintain calibrated
uncertainty estimates. In established implementations of BO, designers place
GP priors on latent functions, $f(\cdot)$, which are fully specified through a
mean function, $m(\bm{x})$, and a covariance function or kernel
$k_{\bm{\theta}}(\bm{x},\bm{x}^{\prime})$ with $\bm{\theta}\in\mathbb{R}^{p}$
representing kernel hyper-parameters. The model specification is completed by
defining a likelihood. Here, practitioners typically assume that observations
$y_{l}$ adhere to a Gaussian noise model such that
$y_{l}=f(\bm{x}_{l})+\epsilon_{l}$ where
$\epsilon_{l}\sim\mathcal{N}(0,\sigma_{\text{noise}}^{2})$. This, in turn,
generates a Gaussian likelihood of the form
$y_{l}|\bm{x}_{l}\sim\mathcal{N}(f_{l},\sigma_{\text{noise}}^{2})$ where we
use $f_{l}$ to denote $f(\bm{x}_{l})$ with
$f(\bm{x})\sim\mathcal{G}\mathcal{P}(m(\bm{x}),k_{\bm{\theta}}(\bm{x},\bm{x}^{\prime}))$.
Additionally, a further design choice commonly made by practitioners is that
the GP kernel is stationary, depending only on the norm between $\bm{x}$ and
$\bm{x}^{\prime}$, $||\bm{x}-\bm{x}^{\prime}||$. From this exposition, we
conclude two important modelling assumptions stated as _data stationarity_ and
_homoscedasticity of the noise distribution_. Where homoscedasticity implies a
constant noise term $\sigma_{\text{noise}}^{2}$. Heteroscedasticity is usually
harder to model as implies $\sigma_{\text{noise}}^{2}$ is a function of the
input: i.e., depending on the data, the noise changes around the mean. Of
course, it is clear that there are significant differences between
homoscedastic functions and heteroscedastic functions, and later we show
indeed heteroscedastic functions require a different approach to optimise over
than the typical homoscedastic (synthetic) functions usually researched in
Bayesian Optimisation. If the true latent process does not adhere to these
assumptions, the resultant model will be a poor approximation to the black-
box. Realising the potential empirical implications of these modelling
choices, we identify the first two questions addressed by this paper:

Q.I. Are hyper-parameter tuning tasks stationary?

Q.II. Are hyper-parameter tuning tasks homoscedastic?

In Section 3.2, we show that even amongst the simplest hyper-parameter tuning
tasks, the null hypothesis may be rejected in the case of statistical
hypothesis tests for heteroscedasticity and non-stationarity.

###  2.2 Acquisition Function & Optimisation Assumptions

Acquisition functions trade off exploration and exploitation by utilising
statistics from the posterior $p_{\bm{\theta}}(f(\cdot)|\mathcal{D})$ with
$\mathcal{D}$ denoting the data (hyper-parameter configurations as inputs and
validation accuracy as outputs) collected so far. Under a GP surrogate with
Gaussian-corrupted observations
$y_{\ell}=f(\boldsymbol{x}_{\ell})+\epsilon_{\ell}$ where
$\epsilon_{\ell}\sim\mathcal{N}(0,\sigma^{2})$, and given a data set
$\mathcal{D}=\\{\boldsymbol{x},\boldsymbol{y}\\}$, the joint distribution of
$\mathcal{D}$ and an arbitrary set of input points $\boldsymbol{x}_{1:q}$ is
given by

$\displaystyle\left[\begin{array}[]{c}\bm{y}\\\
f(\bm{x}_{1:q})\end{array}\right]\Bigg{|}\
\bm{\theta}\sim\mathcal{N}\left(\left[\begin{array}[]{cc}m(\bm{x})\\\
m(\bm{x}_{1:q})\end{array}\right],\left[\begin{array}[]{cc}\bm{K}_{\bm{\theta}}+\sigma^{2}\bm{I}&\bm{k}_{\bm{\theta}}(\bm{x}_{1:q})\\\
\bm{k}^{\mathsf{T}}_{\bm{\theta}}(\bm{x}_{1:q})&\bm{k}_{\bm{\theta}}(\bm{x}_{1:q},\bm{x}_{1:q})\end{array}\right]\right),$

where $\bm{K}_{\bm{\theta}}=\bm{K}_{\bm{\theta}}(\bm{x},\bm{x})$ and
$\bm{k}_{\bm{\theta}}(\bm{x}_{1:q})=\bm{k}_{\bm{\theta}}(\bm{x},\bm{x}_{1:q})$.
From this joint distribution one can derive though marginalisation (?) the
posterior predictive
$p(f(\bm{x}_{1:q})|\mathcal{D})=\mathcal{N}(\bm{\mu}_{\bm{\theta}}(\bm{x}_{1:q}),\bm{\Sigma}_{\bm{\theta}}(\bm{x}_{1:q}))$
with:

$\displaystyle\bm{\mu}_{\bm{\theta}}(\bm{x}_{1:q})$
$\displaystyle=m(\bm{x}_{1:q})+\bm{k}_{\bm{\theta}}(\bm{x}_{1:q})^{\top}(\bm{K}_{\bm{\theta}}+\sigma^{2}\bm{I})^{-1}(\bm{y}-m(\bm{x}))$
$\displaystyle\bm{\Sigma}_{\bm{\theta}}(\bm{x}_{1:q})$
$\displaystyle=\bm{K}_{\bm{\theta}}(\bm{x}_{1:q},\bm{x}_{1:q})-\bm{k}_{\bm{\theta}}(\bm{x}_{1:q})^{\top}(\bm{K}_{\bm{\theta}}+\sigma^{2}\bm{I})^{-1}\bm{k}_{\bm{\theta}}(\bm{x}_{1:q}).$

As such we note that
$p(f(\bm{x}_{1:q})|\mathcal{D})=\mathcal{N}(\bm{\mu}_{\bm{\theta}}(\bm{x}_{1:q}),\bm{\Sigma}_{\bm{\theta}}(\bm{x}_{1:q}))$.
In this paper, we focus on three widely-used myopic acquisition functions
which in a reparameterised form can be written as (?):

Expected Improvement (EI):

$\displaystyle\alpha^{\bm{\theta}}_{\text{EI}}(\bm{x}_{1:q}|\mathcal{D})$
$\displaystyle=\mathbb{E}_{\text{post.}}\Bigg{[}\max_{j\in
1:q}\\{\text{ReLU}(f(\bm{x}_{j})-f(\bm{x}^{+}))\\}\Bigg{]},$

where the subscript ’post.’ is the predictive posterior of a GP (?),
$\bm{x}_{j}$ is the $j^{th}$ vector of $\bm{x}_{1:q}$, and $\bm{x}^{+}$ is the
best performing input in the data so far.

Probability of Improvement (PI):

$\displaystyle\alpha^{\bm{\theta}}_{\text{PI}}(\bm{x}_{1:q}|\mathcal{D})$
$\displaystyle=\mathbb{E}_{\text{post.}}\Bigg{[}\max_{j\in
1:q}\\{\mbox{1{1}}\\{{f}(\bm{x}_{j})-f(\bm{x}^{+})\\}\\}\Bigg{]},$

where $\mbox{1{1}}\\{\cdot\\}$ is the left-continuous Heaviside step function.

Upper Confidence Bound (UCB):

$\displaystyle\alpha^{\bm{\theta}}_{\text{UCB}}(\bm{x}_{j})$
$\displaystyle=\mathbb{E}_{\text{post.}}\Bigg{[}\max_{j\in
1:q}\Bigg{\\{}{\mu}_{\bm{\theta}}(\bm{x}_{j})+\sqrt{\nicefrac{{\beta\pi}}{{2}}}|{\gamma}_{\bm{\theta}}(\bm{x}_{j})|\Bigg{\\}}\Bigg{]},$

where ${\mu}_{\bm{\theta}}(\bm{x}_{j})$ is the posterior mean of the
predictive distribution and
${\gamma}_{\bm{\theta}}(\bm{x}_{j})={f}(\bm{x}_{j})-{\mu}_{\bm{\theta}}(\bm{x}_{j})$.
When it comes to practicality, generic BO implementations make additional
assumptions during the acquisition maximisation step. First, it is assumed
that one of the aforementioned acquisitions works best for a specific task,
and that the GP model is an accurate approximation to the black-box. However,
when it comes to real-world applications, both of these assumptions are
difficult to validate; the best-performing acquisition is challenging to
identify upfront and GP models may easily be misspecified. With this in mind,
we identify a third question that we wish to address:

Q.III. Can acquisition function solutions conflict in hyper-parameter tuning
tasks?

In the following section, we affirm that acquisitions can conflict even on the
simplest of hyper-parameter tuning tasks. Moreover, we show that a robust
formulation to tackle misspecification of acquisition maximisation can improve
overall performance (see Section 4.2.1).

##  3 Modelling Assumption Analysis

Before discussing the improvements afforded to BO via our solution method, we
detail analyses conducted to answer questions (Q.I., Q.II., and Q.III.) posed
in the previous section. Our analyses indicate:

A.I.: Even simple hyper-parameter tuning tasks exhibit significant
heteroscedasticity.
A.II.: Even simple hyper-parameter tasks exhibit significant non-stationarity.
A.III.: Acquisition functions conflict in their optima, occasionally leading
to opposing solutions.

Experiment Setting: We create a wide range of hyper-parameter tasks (108)
across a variety of classification and regression problems. We use nine
models, (e.g. multilayer perceptrons, support vector machines) and six
datasets (two regression and four classification) from the UCI repository, and
two metrics per dataset (such as negative log-likelihood or mean squared
error). Each model possesses tuneable hyper-parameters, e.g. the number of
hidden units of a neural network. The goal is to fit these hyper-parameters so
as to maximise/minimise one of the specified metrics. Values of the black-box
objective are stochastic with noise contributions originating from the train-
test splits used to compute the losses. Experimentation was facilitated by the
Bayesmark111https://github.com/uber/bayesmark package. Full hyper-parameter
search spaces are defined in Table 2 and Table 3. 222It is these search spaces
that are used by the random search baseline..

Statistical Hypothesis Testing for Heteroscedasticity and Non-Stationarity: We
describe here the statistical hypothesis tests we use to answer Q.II.. GP
regression typically considers a conditional normal distribution of the
observations $y|\cdot\sim\mathcal{N}(f(\cdot),\sigma^{2}(\cdot))$ and in most
cases $\sigma(\cdot)^{2}$ is assumed to be constant, in which case the GP is
termed homoscedastic. To assess whether the homoscedasticity assumption holds
for the tasks under examination, we make use of Levene’s test and the Fligner-
Killeen test. To give the reader intuition as to how we apply these tests,
Levene’s test asseses whether the variance is equal in two groups of data,
assuming the data is normally distributed. I.e for a given task, given
multiple evaluations of the black-box of two distinct hyperparameter sets, do
the share the same variance (Homoscedasticity), or do their variances differ
(Heteroscedasticity). Secondly, the Fligner-Killeen test is similarly a test
for Homoscedasticity, however it is particularly useful when the data is non-
normal. We refer the reader to the Appendix A for additional information
regarding the tests.

To run these tests on a given task, we evaluate $k=50$ distinct sets of
hyperparameters $\\{x_{i}\\}_{1\leq i\leq k}$ for $n=10$ times and obtain
scores $\\{Y_{ij}\\}_{1\leq i\leq k,1\leq j\leq n}$, where $Y_{ij}$ is the
$j^{\text{th}}$ score observed when evaluating the $i^{\text{th}}$
configuration. For $i=1,\dots,k$, let $\sigma_{i}^{2}$ denote the observed
variance of $y|x_{i}$, then both Levene’s test and the Fligner-Killeen test
share the same null hypothesis of homoscedasticity:

$H_{0}:\sigma_{1}^{2}=\dots=\sigma_{k}^{2}.$

In all 108 tests, we see a p-value significantly lower than $0.05$ in $72$
tasks using Levene’s test, and in $73$ tasks using Fligner-Killeen test. Such
results (shown in detail in Appendix C) imply that at least $66\%$ of the
experimental tasks exhibit heteroscedastic behaviour.

Table 1: Hypothesis Testing for 108 tasks with respect to GP fit. In the table below we show, out of all 108 tasks, whether the GP fit (marginal log-likelihoods) was improved (Better) when either the Output transform or Input warping was added into the surrogate model, or was worse. Furthermore, we include significantly testing using the one sided t-test and detail how many tasks the GP fit was significantly better or worse with these additional modelling components. We find that output transformations which tackle heteroscedasticity significantly improve GP modelling capabilities in general (improve marginal log-likelihoods). Similarly, input transformations which tackle non-stationarity significantly improve GP modelling capabilities in general. | Better | Sig. Better | Worse | Sig. Worse
---|---|---|---|---
Heteroscedasticity (Output Transform) | 70 (65%) | 58 (54%) | 38 (35%) | 25 (23%)
Non-Stationarity (Input Warping) | 106 (98%) | 79 (73%) | 2 (2%) | 0 (0%)

###  3.1 Answer A.I.: Simple Hyper-parameter Tuning Tasks are Non-Stationary

To assess the impact of the extent of non-stationarity on BO performance, we
conduct probabilistic regression experiments to gauge the predictive
performance of a stationary GP on the hyper-parameter tuning tasks with and
without input warping transformations which correct for non-stationarity. We
first run a one-sided t-test for each of the 108 tasks where the null
hypothesis is that the application of the input warping yields no difference
in the log probability metric. In Table 1 significance tests show that in
106/108 tasks, the log probability metric is more favourable when input
warping is applied. In 79/108 tasks, the gain is significant at the 95% level
of confidence (p-value $<0.025$). It is clear that tackling Non-stationarity
improves GP fit as shown in Table 1 and improves BO performance, as shown by
the algorithm BO Base w Non-stationarity in Figure 4b. We thus conclude that
non-stationarity is an important consideration for BO performance due to the
observed effect on the log probability metric.

###  3.2 Answers A.II.: Simple Hyper-parameter Tuning Tasks are
Heteroscedastic

We perform an analogous hypothesis test as in Section 3.1, assessing a vanilla
GP’s performance with and without output transformations (Box-Cox/ Yeo-
Johnson). We run a two-sided paired t-test for each of the 108 tasks where the
null hypothesis is that the application of the output transform yields no
difference in the log probability metric. In Table 1 significance tests show
that in 70/108 tasks, the log probability metric is more favourable when
output transformations are applied. In 58/108 tasks, the gain is significant
at the 95% level of confidence (p-value $<0.025$). It is clear that tackling
Heteroscedasticity improves GP fit as shown in Table 1 and improves BO
performance, as shown by the algorithm BO Base w Heteroscedasticity in Figure
4b. We thus conclude that heteroscedasticity is an important consideration for
BO performance due to it isimpact on the log probability metric.

Furthermore, to gauge the level of heteroscedasticity in the underlying data,
we use the Fligner-Killeen (?) and Levene (?) tests. For both tests, the null
hypothesis is that the underlying black-box function noise process is
homoscedastic. In all 108 tests, we see a p-value significantly lower than
$0.05$ in $72$ tasks using Levene’s test, and in $73$ tasks using Fligner-
Killeen. Such results (shown in detail in Appendix C) imply that at least
$66\%$ of the experimental tasks exhibit heteroscedastic behaviour.

###  3.3 Answer A.III.: No Clear Winner

Figure 1: Examples depicting conflicting acquisitions across data sets (Wine,
Boston Housing, and Iris) and models (AdaBoost, Multilayer perceptron,
K-Nearest neighbours, and support vector machines). The y-axis shows the
acquisition value, and x-axis a given configuration of hyperparameters.
Clearly, in these examples, not only different acquisitions lead to different
optima, but it can be seen that such solutions might conflict (minimum value
for one acquisition function is a maximum value for another acquisition
function).

It has previously been observed that acquisition functions can conflict in
their optima (?). To provide further support for the answer to Q.III., we
collect 128 samples from each task by evaluating various hyper-parameter
configurations across metrics. We then assemble a data set
$\mathcal{D}=\\{\textbf{hyper-param}_{i},y_{i}\\}_{i=1}^{32}$, where
$\textbf{hyper-param}_{i}$ is a vector with dimensionality dependent on the
number of hyper-parameters in a given model, and $y_{i}$ is an evaluation
metric, (e.g., mean squared error) We subsequently fit a GP surrogate model
and consider each of the three acquisition functions from Section 2.2. Given
the difficulty involved in the graphical depiction of an acquisition function
conflict in more than two dimensions, we examine a simple, two-dimensional
illustrative example. From Figure 1, it is apparent that even in the simplest
2D case, many examples of conflicting acquisitions exist. Thus, in higher
dimensions this behaviour will also occur.

##  4 Optimising Bayesian Optimisation

In this section we describe the component design choices that may mitigate for
heteroscedastic and non-stationary aspects of commonly-encountered BO
problems. Input and output transformations as well as multi-objective
acquisition functions have been introduced in isolation previously, whilst
acquisition function robustness is unique to this work. The overall design
choices produce the method which we name Heteroscedastic and Evolutionary
Bayesian Optimisation (HEBO).

###  4.1 Tackling Heteroscedasticity and Non-Stationarity

To parsimoniously handle heteroscedasticity and non-stationarity, we leverage
ideas from the warped GP literature (?) where output transformations
facilitate the modelling of complex noise processes. We observe that the well-
known Box-Cox (?) and Yeo-Jonhson (?) output transformations in conjunction
with the Kumaraswamy (?) input transformation, offer a balance between
simplicity of implementation and empirical performance. In our ablation study
(Section 5), we demonstrate that the addition of these two modelling
components alone yields large performance gains. Note, we refit the parameters
for the output transformation before we refit the GP after receiving a new
samples.

Output Transformation for Heteroscedasticity: We consider the Box-Cox
transformation most frequently used as a corrective mapping for non-Gaussian
data. The transform depends on a tuneable parameter $\zeta$ and applies the
following map to each of the labels:
$\text{T}_{\zeta}(y_{l})=\nicefrac{{{y}_{l}^{\zeta}-1}}{{\zeta}}$ for
$\zeta\neq 0$ and $\text{T}_{\zeta}(y_{l})=\log y_{l}$ if $\zeta=0$, where in
our case $y_{l}$ denotes the validation accuracy of the $l^{th}$ hyper-
parameter configuration. $\zeta$ must be fit based on the observed data such
that the distribution of the transformed labels closely resembles a Gaussian
distribution. This is achieved by minimising the negative Box-Cox likelihood
function:

$\displaystyle\log\left[\sum_{l=1}^{n}\frac{(\text{T}_{\zeta}(y_{l})-\overline{\text{T}}_{\zeta}(\bm{y}))^{2}}{n}\right]^{\frac{n}{2}}+\sum_{l=1}^{n}\log\left[\text{T}_{\zeta}(y_{l})\right]^{(1-\zeta)},$

where $n$ is the number of datapoints and
$\overline{\text{T}_{\zeta}(\bm{y})}$ is the sample mean of the transformed
labels. Box-Cox transforms only consider strictly positive (or strictly
negative) labels $y_{l}$.

When labels take on arbitrary values, we use the Yeo-Johnson transform in
place of the Box-Cox transform. The Yeo-Johnson transform is defined as
follows:

$\displaystyle\texttt{Y.J.}_{\zeta}(y_{l})=\left\\{\begin{array}[]{lr}\frac{(y_{l}+1)^{\zeta}-1}{\zeta},&\text{if
$\zeta\neq 0$, $y_{l}\geq 0$}\\\ \log(y_{l}+1),&\text{if $\zeta=0$, $y_{l}\geq
0$}\\\ \frac{(1-y_{l})^{2-\zeta}-1}{\zeta-2}&\text{if $\zeta\neq 2$,
$y_{l}<0$}\\\ -\log(1-y_{l})&\text{if $\zeta=2$,
$y_{l}<0$.}\end{array}\right.$

In an analogous fashion to the Box-Cox transform, the Yeo-Johnson’s parameter
is fit based on the observed data through solving the following 1-dimensional
optimisation problem:

$\displaystyle\max_{\zeta}$
$\displaystyle-\frac{n}{2}\log\left[\frac{\sum_{j=1}^{n}(\texttt{Y.J.}_{\zeta}({y}_{l})-\overline{\texttt{Y.J.}_{\zeta}(\bm{y})})^{2}}{n-1}\right]+(\zeta-1)\sum_{i=1}^{n}\left[\text{sign}({y}_{l})\log(|{y}_{l}|+1)\right],$

with $\overline{\texttt{Y.J.}_{\zeta}(\bm{y})}$ the sample average computed
after applying the Yeo-Johnson transformation.

Input Transformations for Non-Stationarity: As a general solution concept for
correcting for non-stationarity, we consider input warping see (?). Input
warping performs a (usually non-linear and learnable) transformation to the
input variables $(\bm{x}_{l})$. It was proven in (?) that Input warping also
helps tackle non-stationary functions. We rely on the Kumaraswamy input
warping transform as used in (?), which operates as follows for each input
dimension:

$[\texttt{Kumaraswamy}_{\bm{\gamma}}(\bm{x}_{l})]_{k}=1-\left(1-[\bm{x}_{l}]_{k}^{a_{k}}\right)^{b_{k}}\
\forall k\in[1:d],$

where $d$ is the dimensionality of the decision variable (i.e. the number of
free hyper-parameters), $a_{k}$ and $b_{k}$ are tuneable warping parameters
for each of the dimensions, and $\bm{\gamma}$ is a vector concatenating all
free parameters, i.e., $\bm{\gamma}=[a_{1:d},b_{1:d}]^{\mathsf{T}}$.
$\bm{\gamma}$ is fit based on the observed data. Similar to (?), we optimise
$\bm{\gamma}$ under the marginal likelihood objective used to fit the GP
surrogate.

All Modelling Improvements Together: Combining the above corrective measures
for heteroscedasticity and non-stationarity leads us to an improved GP
surrogate with more flexible modelling capabilities. The implementation of
such a model is relatively simple and involves maximising a new marginal
likelihood which may be written as:

$\displaystyle\max_{\bm{\theta},\bm{\gamma}}$
$\displaystyle-\frac{1}{2}\text{T}_{\zeta^{\star}}(\bm{y})^{\mathsf{T}}(\bm{K}^{\bm{\gamma}}_{\bm{\theta}}+\sigma_{\text{noise}}\bm{I})^{-1}\text{T}_{\zeta^{\star}}(\bm{y})-\frac{1}{2}|\bm{K}^{\bm{\gamma}}_{\bm{\theta}}+\sigma_{\text{noise}}^{2}\bm{I}|-\text{const,}$

where $\bm{\theta}$ are GP hyper-parameters, $\bm{\gamma}$ indicates the use
of non-stationary transformations, and $\zeta^{\star}$ denotes the solution to
a Box-Cox likelihood objective. It is worth noting that we use Box-Cox as a
running example but as mentioned previously we interchange Box-Cox with Yeo-
Johnson transforms based on the properties of the label $y_{l}$. We use
$\bm{K}^{\bm{\gamma}}_{\bm{\theta}}\in\mathbb{R}^{n\times n}$ to represent a
matrix such that each entry depends on both $\bm{\theta}$ and $\bm{\gamma}$,
where
$k_{\bm{\theta}}^{\gamma}(\bm{x},\bm{x}^{\prime})=k_{\bm{\theta}}(\texttt{Kumaraswamy}_{\bm{\gamma}}(\bm{x}),\texttt{Kumaraswamy}_{\bm{\gamma}}(\bm{x}^{\prime}))$.

###  4.2 Tackling Acquisition Conflict & Robustness

Having proposed modifications to the surrogate model component of the Bayesian
optimisation scheme, we now turn our attention to the acquisition maximisation
step. In particular, we focus on two considerations, the first related to the
assumption of a perfect GP surrogate, and the second centred on conflicting
acquisitions.

####  4.2.1 A Robust Acquisition Objective

As mentioned in Section 2.2, the acquisition maximisation step assumes that an
adequate surrogate model is readily available. During early rounds of training
especially, where data is scarce, such a property is often violated, leading
to potentially severe model misspecification. One way to tackle such model
misspecification is to adopt a robust formulation (?, ?) which attempts to
identify the best-performing query location under the worst-case GP model,
i.e., solving
$\max_{\bm{x}}\min_{\bm{\theta}}\alpha^{\bm{\theta}}(\bm{x}|\mathcal{D})$.
Though such a formulation admits a solution $\bm{x}^{\star}$ that is robust to
worst-case misspecification in $\bm{\theta}$, having a $\max\min$ acquisition
is problematic for several reasons. From a conceptual perspective $\max\min$
formulations are known to lead to very conservative solutions if not correctly
constrained or regularised since the optimiser possesses the power to impair
the GP fit while updating $\bm{\theta}$333One can make a case for augmenting
the objective with a constraint such that updates for $\bm{\theta}$ remain
close to $\bm{\theta}^{\star}$ of the marginal likelihood. The ideal enforced
proximity value however remains unclear in the robust acquisition literature
to date (?, ?).. From the perspective of implementation, one encounters two
further issues. First, no global convergence guarantees are known for the non-
convex, non-concave case (?), and second, ensuring gradients can propagate
through the computation graph restricts surrogates and acquisition functions
to be within the same programming framework.

To avoid worst-case solutions and engender independence between acquisition
functions and surrogate models, given a set of parameters from a trained GP
$\bm{\theta}$, we leverage ideas from domain randomisation (?) and consider an
expected formulation instead over these parameters:
$\max_{\bm{x}}\alpha^{\bm{\theta}}_{\text{rob.}}(\bm{x}|\mathcal{D})\equiv\max_{\bm{x}}\mathbb{E}_{\epsilon\sim\mathcal{N}(\bm{0},\sigma_{\epsilon}^{2}\bm{I})}\left[\alpha^{\bm{\theta}+\epsilon}(\bm{x}|\mathcal{D})\right]$.
Importantly, this problem seeks to find new query locations that perform well
on average over a distribution of surrogate models in favour of assuming a
perfect surrogate. Despite on an intractable nature of
$\alpha^{\bm{\theta}}_{\text{rob.}}(\cdot|\mathcal{D})$, in HEBO we show (the
rigorous representation of this result is presented in Appendix B) that it can
be approximated with any arbitrary precision and high confidence with
$\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})=\alpha^{\bm{\theta}}(\bm{x}|\mathcal{D})+\mathcal{N}(0,\sigma^{2}_{n})$
by properly choosing parameters $\sigma_{\epsilon}$ and $\sigma_{n}$:

Theorem: ( Informal ) Let us consider the stochastic version of the
acquisition function utilised in HEBO and given by
$\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})=\alpha^{\bm{\theta}}(\bm{x}|\mathcal{D})+\mathcal{N}(0,\sigma^{2}_{n})$
and Let
$\alpha^{\bm{\theta}}_{\text{rob.}}(\bm{x}|\mathcal{D})\equiv\mathbb{E}_{\epsilon\sim\mathcal{N}(\bm{0},\sigma_{\epsilon}^{2}\bm{I})}\left[\alpha^{\bm{\theta}+\epsilon}(\bm{x}|\mathcal{D})\right]$
be the robust form of the standard acquisition function given as expectation
over random perturbation of parameter $\bm{\theta}$. Then, with proper choice
of parameters $\sigma_{n}$ and $\sigma_{\epsilon}$, HEBO acquisition function
$\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})$ accurately approximates
the robust acquisition function
$\alpha^{\bm{\theta}}_{\text{rob.}}(\bm{x}|\mathcal{D})$ with high probability
444Here we use the common approach for proving stochastic expressions with
high probability (see (?), (?)). Specifically, we show that for any confidence
parameter $\delta\in(0,1)$ the stochastic expression under consideration is
valid with probability at least $1-\delta$. and for any $\bm{\theta},\bm{x}$.

####  4.2.2 Multi-Objective Acquisition functions

As a final component of our general framework, we propose the use of multi-
objective acquisitions seeking a Pareto-front solution. This formulation
facilitates the process of “hedging” between different acquisitions such that
no single acquisition dominates the solution (?). Formally, we solve

$\max_{\bm{x}}\left(\overline{\alpha}^{\bm{\theta}}_{\text{EI}}(\bm{x}|\mathcal{D}),\overline{\alpha}^{\bm{\theta}}_{\text{PI}}(\bm{x}|\mathcal{D}),\overline{\alpha}^{\bm{\theta}}_{\text{UCB}}(\bm{x}|\mathcal{D})\right),$
(2)

where $\overline{\alpha}^{\bm{\theta}}_{\text{type}}(\bm{x}|\mathcal{D})$ is a
robust acquisition of $\text{type}\in\\{\text{EI},\text{PI},\text{UCB}\\}$ as
introduced in the previous section. We also note that our formulation is
designed to admit the use of a robust objective value of
$\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})=\alpha^{\bm{\theta}}(\bm{x}|\mathcal{D})+\eta_{n}$
with $\eta_{n}$ being a sample from $\mathcal{N}(0,\sigma_{n})$ at each
iteration of the evolutionary solver.

Although solving the problem in Equation 2 is a formidable challenge, we note
the existence of many mature multi-objective optimisation algorithms. These
range from first-order (?) to zero-order (?, ?) and evolutionary methods (?,
?). Due to the discrete nature of hyper-parameters in machine learning tasks,
we advocate the use of evolutionary solvers that naturally handle categorical
and integer-valued variables. In our experiments, we employ the non-dominated
sorting genetic algorithm II (NSGA-II) which allows for mixed variable
crossover and mutation to optimise real-valued and integer-valued inputs (?).
We use the implementation of NSGA-II found in the Pymoo (?) library.
Alternatively, one may use the GP Hedge acquisition as used in Dragonfly (?)
in (?) or in SkOpt to select between acquisitions. We however, observed this
formulation to perform poorly when compared against individual acquisitions.

##  5 Experiments and Results

Figure 2: HEBO compared against all baselines for 16 iterations and a batch
size of 8 query points per iteration. Each experiment is repeated with 20
random seeds. We average each seed over both metrics in all tasks and display
a subset of 36 summary plots for the 108 black-box functions. HEBO achieves
the highest normalised mean score in $68.5\%$ of the 108 black-box functions.
Full results for the 108 tasks are presented in Appendix D in tabular format.

In this section, we continue our empirical evaluation and validate gains (if
any) that arise from the improvements proposed in Section 4. The experimental
setup remains as described in Section 3. To assess performance, we use the
normalised task score555Note, we don’t report the time to compute query points
per algorithm as this was under 20 seconds per query point batch.. We run
experiments on either 16 iterations with a batch of 8 query points per
iteration or 100 iterations with 1 query point. Each experiment is repeated
for 20 random seeds. We baseline against a wide range of solvers that either
rely on BO-strategies or follow zero-order techniques such as differential
evolution or particle swarms. These include SkOpt (?)
666https://github.com/scikit-optimize/scikit-optimize pySOT
777https://github.com/dme65/pySOT a parallel global optimisation package (?),
HyperOpt (?) 888https://github.com/hyperopt/hyperopt, OpenTuner
999https://github.com/jansel/opentuner a package for ensembling methods (?),
NeverGrad (?) 101010https://github.com/facebookresearch/nevergrad a gradient-
free optimisation toolbox (where we use the One Plus One optimiser with the
associated label NeverGrad (1+1)), BOHB (?) 111111
https://github.com/automl/HpBandSter and Dragonfly (?) 121212
https://github.com/dragonfly/dragonfly. Additionally, we carried our modelling
improvements to TuRBO 131313
https://github.com/rdturnermtl/bbo_challenge_starter_kit/ (?), augmenting the
standard GP with mitigation strategies from Section 4 producing a new baseline
that we entitle TuRBO+. Finally, we introduce Heteroscedastic Evolutionary
Bayesian Optimisation (HEBO), in which we construct an optimiser with the
improvements introduced in Section 4.

Implementation Details for BOHB: BOHB is a scalable hyper-parameter tuning
algorithm introduced in (?) mixing bandits and BO approaches to achieve both
competitive anytime and final performance. Contrary to the other solvers
considered in this paper, BOHB is specifically designed to tackle multi-
fidelity optimisation and uses the Hyperband (?) routine to define the
fidelity levels under which points are asynchronously evaluated. The selection
of points follows a BO strategy based on the Tree Parzen Estimator (TPE)
method. Given a data set $\mathcal{D}$ of observed data points and a threshold
$\alpha\in\mathbb{R}$, the TPE models $p(\bm{x}|y)$, using kernel density
estimates of

$\displaystyle\ell(\bm{x})$ $\displaystyle=p(y<\alpha|\bm{x},\mathcal{D})$
$\displaystyle g(\bm{x})$ $\displaystyle=p(y\geq\alpha|\bm{x},\mathcal{D}).$

In the TPE algorithm, maximising the expected improvement criterion

$\alpha_{\text{EI}}(\boldsymbol{x})=\int\max(0,\alpha-p(y|\boldsymbol{x}))p(y|\boldsymbol{x})dy$

is equivalent to maximising the ratio $r(x)=\tfrac{\ell(\bm{x})}{g(\bm{x})}$
which is carried out to select a single new candidate point at a time.

In the absence of a multi-fidelity setup in our experiments, we run a modified
version of the BOHB algorithm implemented in the HpBandSter package. We leave
the TPE method for modelling unchanged but ignore the fidelity level
assignment from Hyperband. Moreover, as our experimental setup involves batch
acquisitions, we tested two alternatives to the standard BOHB acquisition
procedure to support synchronous suggestion of multiple points. In the first
approach, we run $q$ independent maximisation processes of $r(\bm{x})$ from
random starting points and recover a single candidate from each process to
form the $q$-batch suggestion. In the second approach, we obtain one point as
a result of a single maximisation of $r(\bm{x})$ and we sample $q-1$ random
points to complete the $q$-batch suggestion. As the latter method yields
better overall performance, the results reported under the BOHB-BB label are
obtained using the second approach.

###  5.1 Black-Box Functions

As discussed in Section 3, we evaluate black-box optimisation solvers on a
large set of tasks from the Bayesmark package. Each task involves optimising
the hyper-parameters of a machine learning algorithm to minimise the cross
validation loss incurred when this model is applied in a regression (reg) or a
classification (clf) setting for a given data set. Thus, a task is
characterised by a model, a data set and a loss function (metric) quantifying
the quality of the regression or classification performance. In total, $108$
distinct tasks can be defined from the valid combinations of the nine models
specified in Table 2, the following six real-world UCI datasets (?), Boston
(reg), Breast Cancer (clf), Diabetes (reg), Digits (clf), Iris (clf) and Wine
(clf); the following two regression metrics, negative mean-squared error
(MSE), negative mean absolute error (MAE), and two classification metrics,
negative log-likelihood (NLL) and negative accuracy (ACC). The results
reported in Figures 3 and 4 have been obtained by applying each black-box
optimisation method using $16$ iterations of $8$-batch acquisition steps on
all of the $108$ tasks. In order to provide a reliable evaluation of the
different solvers, we repeated each run with $20$ random seeds and considered
the normalised score given by:

$\textbf{Normalised
Score}=100\times\frac{\mathcal{L}-\mathcal{L}^{*}}{\mathcal{L}^{\text{rand}}-\mathcal{L}^{*}}$
(3)

where $\mathcal{L}$ is the best-achieved cross validation loss at the end of
the $16$ acquisition steps, $\mathcal{L}^{*}$ is the estimated optimal loss
for the task and $\mathcal{L}^{\text{rand}}$ is the mean loss (across multiple
runs) obtained using random search with the same number of acquisition steps.
The normalisation procedure permits aggregation of the scores across tasks
despite the different cross-validation loss functions used.

###  5.2 Black-Box Optimisation Input Variables

We provide in Table 2 and Table 3 the list of the hyper-parameters controlling
the behaviour of each model along with their optimisation domains, which can
differ depending on whether the model is used for a classification or a
regression task. The search domain may include a mix of continuous and
integer-valued variables (e.g. the MLP-SGD hyper-parameter set includes an
integer-valued hidden layer size, and a continuous-valued initial learning
rate that can take on values between $10^{-5}$ and $10^{-1}$). The
dimensionality of the input space, i.e. the number of hyper-parameters to
tune, ranges from $2$ to $9$. We specify in the final column of the tables
whether the search domain is modified through a standard transformation (
$\log$ or $\operatorname{logit}$) in order to facilitate optimisation.

Table 2: Search spaces for hyper-parameter tuning on classification tasks. We specify the variable type of each hyper-parameter (with $\mathbb{R}$ for real-valued and $\mathbb{Z}$ for integer- valued) as well as the search domain. We specify $\log-\mathcal{U}$ (resp. $\text{logit}-\mathcal{U}$) to indicate that a $\log$ (resp. logit) transformation is applied to the optimisation domain. Model | Parameter | Type | Domain
---|---|---|---
kNN | n_neighbors | $\mathbb{Z}$ | $\mathcal{U}(1,25)$
| p | $\mathbb{Z}$ | $\mathcal{U}(1,4)$
Support Vector Machine | C | $\mathbb{R}$ | $\log-\mathcal{U}(1,10^{3})$
| gamma | $\mathbb{R}$ | $\log-\mathcal{U}(10^{-4},10^{-3})$
| tol | $\mathbb{R}$ | $\log-\mathcal{U}(10^{-5},10^{-1})$
Decision Tree | max_depth | $\mathbb{Z}$ | $\mathcal{U}(1,15)$
| min_samples_split | $\mathbb{R}$ | $\text{logit}-\mathcal{U}(0.01,0.99)$
| min_samples_leaf | $\mathbb{R}$ | $\text{logit}-\mathcal{U}(0.01,0.49)$
| min_weight_fraction_leaf | $\mathbb{R}$ | $\text{logit}-\mathcal{U}(0.01,0.49)$
| max_features | $\mathbb{R}$ | $\text{logit}-\mathcal{U}(0.01,0.99)$
| min_impurity_decrease | $\mathbb{R}$ | $\mathcal{U}(0,0.5)$
Random Forest | max_depth | $\mathbb{Z}$ | $\mathcal{U}(1,15)$
| max_features | $\mathbb{R}$ | $\text{logit}-\mathcal{U}(0.01,0.99)$
| min_samples_split | $\mathbb{R}$ | $\text{logit}-\mathcal{U}(0.01,0.99)$
| min_samples_leaf | $\mathbb{R}$ | $\text{logit}-\mathcal{U}(0.01,0.49)$
| min_weight_fraction_leaf | $\mathbb{R}$ | $\text{logit}-\mathcal{U}(0.01,0.49)$
| min_impurity_decrease | $\mathbb{R}$ | $\mathcal{U}(0,0.5)$
MLP-Adam | hidden_layer_sizes | $\mathbb{Z}$ | $\mathcal{U}(50,200)$
| alpha | $\mathbb{R}$ | $\log-\mathcal{U}(10^{-5},10^{1})$
| batch_size | $\mathbb{Z}$ | $\mathcal{U}(10,250)$
| learning_rate_init | $\mathbb{R}$ | $\log-\mathcal{U}(10^{-5},10^{-1})$
| tol | $\mathbb{R}$ | $\log-\mathcal{U}(10^{-5},10^{-1})$
| validation_fraction | $\mathbb{R}$ | $\text{logit}-\mathcal{U}(0.1,0.9)$
| beta_1 | $\mathbb{R}$ | $\text{logit}-\mathcal{U}(0.5,0.99)$
| beta_2 | $\mathbb{R}$ | $\text{logit}-\mathcal{U}(0.9,1-10^{-6})$
| epsilon | $\mathbb{R}$ | $\log-\mathcal{U}(10^{-9},10^{-6})$
MLP-SGD | hidden_layer_sizes | $\mathbb{Z}$ | $\mathcal{U}(50,200)$
| alpha | $\mathbb{R}$ | $\log-\mathcal{U}(10^{-5},10^{1})$
| batch_size | $\mathbb{Z}$ | $\mathcal{U}(10,250)$
| learning_rate_init | $\mathbb{R}$ | $\log-\mathcal{U}(10^{-5},10^{-1})$
| power_t | $\mathbb{R}$ | $\text{logit}-\mathcal{U}(0.1,0.9)$
| tol | $\mathbb{R}$ | $\log-\mathcal{U}(10^{-5},10^{-1})$
| momentum | $\mathbb{R}$ | $\text{logit}-\mathcal{U}(0.001,0.999)$
| validation_fraction | $\mathbb{R}$ | $\text{logit}-\mathcal{U}(0.1,0.9)$
AdaBoost | n_estimators | $\mathbb{Z}$ | $\mathcal{U}(10,100)$
| learning_rate | $\mathbb{R}$ | $\log-\mathcal{U}(10^{-4},10^{1})$
Lasso | C | $\mathbb{R}$ | $\log-\mathcal{U}(10^{-2},10^{2})$
| intercept_scaling | $\mathbb{R}$ | $\log-\mathcal{U}(10^{-2},10^{2})$
Linear | C | $\mathbb{R}$ | $\log-\mathcal{U}(10^{-2},10^{2})$
| intercept_scaling | $\mathbb{R}$ | $\log-\mathcal{U}(10^{-2},10^{2})$
Table 3: Models and search spaces for hyper-parameter tuning on regression tasks. Models having the same search spaces for classification and regression tasks are omitted (cf. Table 2). Model | Parameter | Type | Domain
---|---|---|---
AdaBoost | n_estimators | $\mathbb{Z}$ | $\mathcal{U}(10,100)$
| learning_rate | $\mathbb{R}$ | $\log-\mathcal{U}(10^{-4},10^{1})$
Lasso | alpha | $\mathbb{R}$ | $\log-\mathcal{U}(10^{-2},10^{2})$
| fit_intercept | $\mathbb{Z}$ | $\mathcal{U}(0,1)$
| normalize | $\mathbb{Z}$ | $\mathcal{U}(0,1)$
| max_iter | $\mathbb{Z}$ | $\log-\mathcal{U}(10,5000)$
| tol | $\mathbb{R}$ | $\log-\mathcal{U}(10^{-5},10^{-1})$
| positive | $\mathbb{Z}$ | $\mathcal{U}(0,1)$
Linear | alpha | $\mathbb{R}$ | $\log-\mathcal{U}(10^{-2},10^{2})$
| fit_intercept | $\mathbb{Z}$ | $\mathcal{U}(0,1)$
| normalize | $\mathbb{Z}$ | $\mathcal{U}(0,1)$
| max_iter | $\mathbb{Z}$ | $\log-\mathcal{U}(10,5000)$
| tol | $\mathbb{R}$ | $\log-\mathcal{U}(10^{-4},10^{-1})$

Table 4 synthesises the performance achieved on the $108$ tasks by the black-
box optimisation solvers considered in our experiments. We note that the
distribution of the scores attained by HEBO has the largest mean and the
smallest standard deviation, indicating that HEBO significantly outperforms
competitor algorithms.

Algorithm | Mean | Std | Median | 40th Centile | 30th Centile | 20th Centile | 5th Centile
---|---|---|---|---|---|---|---
HEBO | $\bm{100.12}$ | $\bm{8.70}$ | $\bm{100.01}$ | $\bm{100.00}$ | $\bm{99.88}$ | $\bm{98.64}$ | $\bm{85.71}$
PySOt | 98.18 | 9.03 | 100.00 | 99.81 | 98.60 | 95.36 | 80.00
TuRBO | 97.95 | 10.80 | 100.00 | 99.88 | 98.75 | 95.26 | 78.63
HyperOpt | 96.37 | 8.79 | 99.31 | 98.16 | 95.94 | 92.38 | 78.52
SkOpt | 96.18 | 11.51 | 99.78 | 98.66 | 96.73 | 91.62 | 74.77
TuRBO+ | 95.29 | 10.93 | 98.97 | 97.60 | 95.27 | 90.92 | 74.77
OpenTuner | 94.32 | 14.18 | 98.44 | 96.93 | 93.84 | 89.97 | 68.96
Nevergrad (1+1) | 93.20 | 17.52 | 99.65 | 97.84 | 94.57 | 88.28 | 55.34
BOHB | 92.03 | 11.16 | 96.02 | 93.55 | 90.14 | 85.71 | 67.82
Random-Search | 92.00 | 11.71 | 96.18 | 93.55 | 90.05 | 85.16 | 69.55
Table 4: Mean and n-th percentile normalised scores over $108$ black-box
functions, each repeated with 20 random seeds. We observe significant mean
improvements from HEBO compared to all competitor algorithms.  (a) Empirical
Performance Gain Figure 3: Analysis of the results on 108 tuning tasks. (Left)
Normalised score comparison demonstrating that HEBO (i.e., BO with
improvements from Section 4) outperforms competitor algorithms. We observe a
5% relative improvement to SOTA optimisers such as TuRBO. (Right) HEBO yields
an 8% improvement compared to random search.

Figure 2 and Figure 3 demonstrates gains from adopting the general HEBO
framework. We note that due to optimising over numerous regression and
classification metrics, we show that irrelevant of the validation score HEBO
performs better than other optimisers. In Figure 4, we compare HEBO against
baselines and report up to an $8\%$ performance gain relative to a random
search strategy. It is also worth noting that TuRBO+ tends to underperform
141414We believe this due to the trust region not being modelled correctly
with input warping., achieving ca. $4\%$ improvement relative to random
search. We believe such a result is related to the interplay between our
approach’s capabilities to address heteroscedasticity and non-stationarity as
well as the size of the trust regions; an interesting avenue that we plan to
explore in future work, as well as experimenting with deeper neural networks
as well as other architectures such as convolutional/ recurrent neural
networks. Overall, HEBO achieves the highest normalised mean scores on 74 of
the 108 datasets. Complete results on all tasks may be found in Appendix D.

Comparison to Asynchronous BO Algorithms: We perform a comparison to black-box
optimisers, such as Dragonfly and BOHB, which operate in the asynchronous
setting. We run each method for 100 iterations of data collection with a
single query location per iteration. We label the asynchronous algorithms
without their multi-fidelity components with an addition BB for black-box
optimiser (Dragonfly-BB and BOHB-BB) to assess black-box optimisation
performance only. The results of Figure 4a show that in the asynchronous
setting, both Dragonfly-BB and BOHB-BB under-perform relative to other black-
box optimisers, with HEBO performing best. However, this result is not
surprising as asynchronous methods trade off sample efficiency with speed.
Nevertheless, this experiment reveals a large gap in suggestion power between
SOTA asynchronous and synchronous methods.

(a) Competitor Comparison.

(b) Ablation Study

Figure 4: (a) We compare HEBO against several popular hyper-parameter tuning
approaches including BOHB-BB and Dragonfly-BB, running all methods for 100
iterations with a batch size of 1 (i.e. one set of hyper-parameters queried
per iteration). BOHB-BB and Dragonfly-BB feature asynchronous queries,
suggesting a batch of one set of hyper-parameters at each iteration. We remove
the multi-fidelity components from BOHB and Dragonfly to assess Black-Box
optimisation alone, hence the additional BB appended to their label. (b)
Ablation study where X denotes a general component of HEBO. HEBO w/o X takes
one component X out at a time and BO Base w X adds one component X in at a
time. We show TuRBO as a baseline and refer to HEBO with all significant
components removed as BO Base. The ablation demonstrates that the corrections
for each misspecified modelling assumption yield a tangible gain in empirical
performance.

###  5.3 Ablation Results

To better understand the relative importance of each component of the HEBO
algorithm, we conduct an ablation study by first removing each component of
HEBO and testing the remaining components and second, by starting with basic
BO and sequentially adding and testing each component of HEBO. The components
comprise the consideration of heteroscedasticity, non-stationarity and
robustness, as well as the use of a multiobjective acquisition function. We
report average normalised scores in Figure 4b. The precedence order observed
is: heteroscedasticity, multi-objective acquisition functions, non-
stationarity and robustness.

##  6 Related Work

We introduce work on the following topics relating to modelling, acquisition
and optimisers in Bayesian optimisation:

Heteroscedasticity with output transforms: Among various approaches to
handling heteroscedasticity (?, ?, ?, ?, ?), transforming the output variables
is a straightforward option giving rise to warped Gaussian processes (?). More
recently, output transformations have been extended to compositions of
elementary functions (?) and normalising flows (?, ?). Output transformations
have not featured prominently in the Bayesian optimisation literature, perhaps
due to the commonly-held opinion that warped GPs require more data relative to
standard GPs in order to function as effective surrogates (?). Rather than
introduce additional hyper-parameters to the GP, we enable efficient output
warping through methods that only require pre-training. Recent work (?) has
also investigated Gaussian copula transforms which may prove to be
particularly effective in situations where there are outliers.

Non-stationarity with input warping: Many surrogate models with input warping
exist for optimising non-stationary black-box objectives (?, ?, ?) and have
enjoyed particular success in hyper-parameter tuning where the natural scale
of parameters is often logarithmic. Traditionally, a Beta cumulative
distribution function is used. In this paper, we adopt the Kumaraswamy warping
which is another instance of the generalised Beta class of distributions which
we have observed to achieve superior performance (?) 151515For clarity we note
that the input warping function used in (?) is the same one used in this
work.; confirming results reported in (?).

Multi-objective acquisition ensembles: Multi-objective acquisition ensembles
were first proposed in (?) and are closely related to portfolios of
acquisition functions (?, ?, ?). In this form, the optimisation problem
involves at least two conflicting and expensive black-box objectives and as
such, solutions are located along the Pareto-efficient frontier. The multi-
objective acquisition ensemble employs these ideas to find a Pareto-efficient
solution amongst multiple acquisition functions. Although we utilised the
multi-objective acquisition ensemble, we note that our framework is solver
agnostic in so far as any multi-objective optimiser (?) may be applied.

Robustness of Acquisitions: Methods achieving robustness with respect to
either surrogates (?) or the optimisation process (?, ?) have been previously
proposed. Most relevant to our setting, is the approach of (?) that introduces
robustness to BO by solving a $\max\min$ objective to determine optimal input
perturbations. Their method, however, relies on gradient ascent-descent-type
algorithms that require real-valued variables and are not guaranteed to
converge in the general non-convex, non-concave setting (?). On the other
hand, our solution possesses two advantages: 1) simplicity of implementation
as we merely require random perturbations of acquisition functions to
guarantee robustness, and 2) support for mixed variable solutions through the
use of evolutionary solvers.

##  7 Conclusion & Future Work

In this paper, we presented an in-depth empirical study of Bayesian
optimisation for hyper-parameter tuning tasks. We demonstrated that even the
simplest among machine learning problems can exhibit heteroscedasticity and
non-stationarity. We also reflected on the affects of misspecified models and
conflicting acquisition functions. We augmented BO algorithms with various
enhancements and revealed that with a revised set of assumptions BO can in
fact act as a competitive baseline in hyper-parameter tuning. We highlight the
large discrepancy between suggestion power of synchronous and asynchronous
methods. We hope for future work to focus on integrating the best of
asynchronous and synchronous methods for optimal performance. We hope this
paper’s findings can guide the community when employing black-box and Bayesian
optimisation in practice.

##  A Addition Detail Of Hypothesis Tests

Levene’s Test Levene’s test statistic is defined as

$W=\frac{N-k}{k-1}\cdot\frac{\sum_{i=1}^{k}n(\bar{Z}_{i\cdot}-\bar{Z}_{\cdot\cdot})^{2}}{\sum_{i=1}^{k}\sum_{j=1}^{n}(Z_{ij}-\bar{Z}_{i\cdot})^{2}},$

where $N=k\times n$, $Z_{ij}=|Y_{ij}-\tfrac{1}{n}\sum_{j=1}^{n}Y_{ij}|$,
$\bar{Z}_{i\cdot}=\tfrac{1}{n}\sum_{j=1}^{n}Z_{ij}$ and
$\bar{Z}_{\cdot\cdot}=\tfrac{1}{k}\sum_{i=1}^{k}\bar{Z}_{i\cdot}$, for all
$i=1,\dots,k$, $j=1,\dots,n$. Levene’s test rejects the null hypothesis of
homoscedasticity $H_{0}$ if

$W>F_{\alpha,k-1,N-k},$

where $F_{\alpha,k-1,N-k}$ is the upper critical value at a significance level
$\alpha$ of the $F$ distribution with $k-1$ and $N-k$ degrees of freedom. The
Fligner-Killeen test is an alternative to Levene’s test that is particularly
robust to outliers.

Fligner-Killeen Test: Computation of the Fligner-Killeen test involves ranking
all the absolute values $\\{|Y_{ij}-\tilde{Y}_{i}|\\}_{1\leq i\leq k,1\leq
j\leq n}$, where $\tilde{Y}_{i}$ is the median of $\\{Y_{ij}\\}_{1\leq j\leq
n}$. Increasing scores
$a_{N,r}=\Phi^{-1}\left(\tfrac{1+\tfrac{r}{N+1}}{2}\right)$ are associated
with each rank $r=1,\dots,N$, where $N=kn$ and $\Phi(\cdot)$ is the cumulative
distribution function for a standard normal random variable. We denote the
rank score associated with $Y_{ij}$ as $r_{ij}$. The Fligner-Killeen test
statistic is given by

$\chi_{o}^{2}=\frac{\sum_{i=1}^{k}n\left(\bar{A}_{i}-\bar{a}\right)^{2}}{V^{2}},$

where $\bar{A}_{i}=\tfrac{1}{n}\sum_{j=1}^{n}a_{N,r_{i}j}$,
$\bar{a}=\tfrac{1}{N}\sum_{r=1}^{N}a_{N,r}$ and
$V^{2}=\tfrac{1}{N-1}\sum_{r=1}^{N}(a_{N,r}-\bar{a})^{2}$. As $\chi_{0}$ has
an asymptotic $\mathcal{X}^{2}$ distribution with $(k-1)$ degrees of freedom,
the test rejects the null hypothesis of homoscedasticity $H_{0}$ if

$\chi_{0}>\mathcal{X}^{2}_{\alpha,k-1}$

where $\mathcal{X}^{2}_{\alpha,k-1}$ is the upper critical value at a
significance level $\alpha$ of the $\mathcal{X}^{2}$ distribution with $k-1$
degrees of freedom.

##  B Details of Robust Acquisition

Though appealing, our formulation still assumes access to the GP hyper-
parameters which complicates the implementation by restricting models and
optimisers to the same programming paradigm. Ideally, we would wish to illicit
robustness through only the GP predictive mean and predictive variance.
Fortunately, we are able to show that upon a simple acquisition perturbation
it becomes possible to approximate $\alpha_{\text{rob}}(\cdot)$ above. As
such, we demonstrate that robust acquisition formulations are achievable using
only the GP predictive mean and variance.

Theorem: Let us consider the stochastic version of the acquisition function
utilised in HEBO and given by
$\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})=\alpha^{\bm{\theta}}(\bm{x}|\mathcal{D})+\eta\sigma_{n}$
with $\eta\sim\mathcal{N}(0,1)$ and standard deviation parameter
$\sigma_{n}>0$ 161616We note that gradient-based algorithms remain applicable
upon addition of the $\eta\sigma_{n}$ term. In our formulation however, we use
an evolutionary method which utilises acquisition function values.
Consequently, the path followed by the optimiser will be altered based on
$\eta$ samples leading to more robust query locations.. Let
$\alpha^{\bm{\theta}}_{\text{rob.}}(\bm{x}|\mathcal{D})\equiv\mathbb{E}_{\epsilon\sim\mathcal{N}(\bm{0},\sigma_{\epsilon}^{2}\bm{I})}\left[\alpha^{\bm{\theta}+\epsilon}(\bm{x}|\mathcal{D})\right]$
be the robust form of the standard acquisition function given as expectation
over random perturbation of parameter $\bm{\theta}$. Then, by properly
choosing parameters $\sigma_{n}$ and $\sigma_{\epsilon}$ with high
probability171717Here we use the common approach for proving stochastic
expressions with high probability (see (?), (?)). Specifically, we show that
for any confidence parameter $\delta\in(0,1)$ the stochastic expression under
consideration is valid with probability at least $1-\delta$., HEBO acquisition
function $\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})$ accurately
approximates the robust acquisition function
$\alpha^{\bm{\theta}}_{\text{rob.}}(\bm{x}|\mathcal{D})$ for any
$\bm{\theta},\bm{x}$. More formally, for any $\rho\in(0,1)$ and
$\delta\in(0,1)$, there are parameters $\sigma_{n}=\sigma_{n}(\rho,\delta)$
and $\sigma_{\epsilon}=\sigma_{\epsilon}(\rho,\delta)$ such that:

$\left|\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})-\alpha^{\bm{\theta}}_{\text{rob.}}(\bm{x}|\mathcal{D})\right|\leq\rho,\
\ \ \ \ \ \forall\bm{\theta,x}$

with probability at least $1-\delta$.

####  B.0.1 Proof of the Robustness Bound

Let $\delta\in(0,1)$ be the desired probability threshold, and $\rho\in(0,1)$
be a desired accuracy parameter. Consider the GP with mean function
$m(\bm{x})$ and covariance function $k_{\bm{\theta}}(\bm{x},\bm{x}^{\prime})$
such that $\forall\bm{x},\bm{x}^{\prime}\in\mathcal{X}$,
$\bm{\theta}\in\mathbb{R}^{p}$:

$\displaystyle|k_{\bm{\theta}}(\bm{x},\bm{x})|\geq M_{0},\ \
|k_{\bm{\theta}}(\bm{x},\bm{x}^{\prime})|\leq M_{1},$ (4)
$\displaystyle||\nabla_{\boldsymbol{\theta}}k_{\bm{\theta}}(\bm{x},\bm{x}^{\prime})||_{2}\leq
M_{2},\ \ |m(\bm{x})|\leq M_{4}.$

Moreover, assume that observations $y\in\mathcal{D}$ are bounded, i.e.
$|y|\leq C$ and let
$\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})=\alpha^{\bm{\theta}}(\bm{x}|\mathcal{D})+\eta\sigma_{n}$
with $\eta$ a standard normal random variable. Then, we are going to show that
there are values $c_{1}=c_{1}(\rho,\delta)$ and $c_{2}=c_{2}(\rho,\delta)$,
such that choosing $\sigma_{n}\leq c_{1}$ and $\sigma_{\epsilon}\leq c_{2}$:

$\left|\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})-\mathbb{E}_{\epsilon\sim\mathcal{N}(\bm{0},\sigma_{\epsilon}^{2}\bm{I})}\left[\alpha^{\bm{\theta}+\epsilon}(\bm{x}|\mathcal{D})\right]\right|\leq\rho.$

with probability at least $1-\delta$. Note, the robust form of the acquisition
function given as
$\alpha^{\bm{\theta}}_{\text{rob.}}(\bm{x}|\mathcal{D})\equiv\mathbb{E}_{\epsilon\sim\mathcal{N}(\bm{0},\sigma_{\epsilon}^{2}\bm{I})}\left[\alpha^{\bm{\theta}+\epsilon}(\bm{x}|\mathcal{D})\right]$
constitutes an intractable integral. Therefore, in order to be optimised
during the course of Bayesian optimisation, the intractable integral must be
replaced by an accurate approximation. Without loss of generality we choose
the UCB acquisition function
$\alpha^{\bm{\theta}}(\bm{x}|\mathcal{D})=\alpha_{\text{UCB}}^{\bm{\theta}}(\bm{x}|\mathcal{D})$
and to avoid technical complications relating to multivariate calculus we
consider a batch size $q=1$. In this case, the UCB acquisition function can be
written as
$\alpha_{\text{UCB}}^{\bm{\theta}}(\bm{x}|\mathcal{D})=\mu_{\bm{\theta}}(\bm{x}|\mathcal{D})+\sqrt{\frac{\beta\pi}{2}}\sigma_{\boldsymbol{\theta}}(\bm{x}|\mathcal{D})$,
where $\mu_{\bm{\theta}}(\bm{x}|\mathcal{D})$ and
$\sigma_{\boldsymbol{\theta}}(\bm{x}|\mathcal{D})$ are the posterior mean and
posterior standard deviation respectively. Consider a Monte-Carlo estimation
of
$\alpha^{\bm{\theta}}_{\text{rob.}}(\bm{x}|\mathcal{D})\equiv\mathbb{E}_{\epsilon\sim\mathcal{N}(\bm{0},\sigma_{\epsilon}^{2}\bm{I})}\left[\alpha^{\bm{\theta}+\epsilon}(\bm{x}|\mathcal{D})\right]$:

$\displaystyle\hat{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})=\frac{1}{N_{\epsilon}}\sum_{j=1}^{N_{\epsilon}}\alpha^{\bm{\theta}+\bm{\epsilon}_{j}}(\bm{x}|\mathcal{D})$

where $\bm{\epsilon}_{j}$ are i.i.d samples drawn from
$\mathcal{N}(\bm{0},\sigma^{2}_{\epsilon}\bm{I})$. Then, adding and
subtracting $\hat{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})$ gives:

$\displaystyle\left|\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})-\mathbb{E}_{\epsilon\sim\mathcal{N}(\bm{0},\sigma_{\epsilon}^{2}\bm{I})}\left[\alpha^{\bm{\theta}+\epsilon}(\bm{x}|\mathcal{D})\right]\right|\leq$
$\displaystyle\left|\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})-\hat{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})\right|+\left|\hat{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})-\mathbb{E}_{\epsilon\sim\mathcal{N}(\bm{0},\sigma_{\epsilon}^{2}\bm{I})}\left[\alpha^{\bm{\theta}+\epsilon}(\bm{x}|\mathcal{D})\right]\right|.$

Using the definition of $\hat{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})$ in
the above result yields:

$\displaystyle\left|\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})-\mathbb{E}_{\epsilon\sim\mathcal{N}(\bm{0},\sigma_{\epsilon}^{2}\bm{I})}\left[\alpha^{\bm{\theta}+\epsilon}(\bm{x}|\mathcal{D})\right]\right|\leq$
(5)
$\displaystyle\frac{1}{N_{\epsilon}}\sum_{j=1}^{N_{\epsilon}}\left|\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})-\alpha^{\bm{\theta}+\bm{\epsilon}_{j}}(\bm{x}|\mathcal{D})\right|+\left|\frac{1}{N_{\epsilon}}\sum_{j=1}^{N_{\epsilon}}\alpha^{\bm{\theta}+\bm{\epsilon}_{j}}(\bm{x}|\mathcal{D})-\mathbb{E}_{\epsilon\sim\mathcal{N}(\bm{0},\sigma_{\epsilon}^{2}\bm{I})}\left[\alpha^{\bm{\theta}+\epsilon}(\bm{x}|\mathcal{D})\right]\right|$

Let us now study separately each term in the above result. Applying the
Chebyshev inequality for the second term in the above expression, we have that
with probability at least
$p_{1}=1-\frac{8\left[\mathbb{E}_{\bm{\epsilon}}\left[\mu^{2}_{\bm{\theta}+\bm{\epsilon}}(\bm{x}|\mathcal{D})\right]+\frac{\beta\pi}{2}\mathbb{E}_{\bm{\epsilon}}\left[\sigma^{2}_{\bm{\theta}+\bm{\epsilon}}(\bm{x}|\mathcal{D})\right]\right]}{N_{\epsilon}\rho^{2}}$:

$\displaystyle\left|\frac{1}{N_{\epsilon}}\sum_{j=1}^{N_{\epsilon}}\alpha^{\bm{\theta}+\bm{\epsilon}_{j}}(\bm{x}|\mathcal{D})-\mathbb{E}_{\epsilon\sim\mathcal{N}(\bm{0},\sigma_{\epsilon}^{2}\bm{I})}\left[\alpha^{\bm{\theta}+\epsilon}(\bm{x}|\mathcal{D})\right]\right|\leq\frac{\rho}{2}.$
(6)

In order to ensure that $p_{1}=1-\frac{\delta}{2}$, the sample number
$\bm{\epsilon}_{j}$ should be taken:

$\displaystyle
N_{\epsilon}=\left\lceil\frac{16\left[\mathbb{E}_{\bm{\epsilon}}\left[\mu^{2}_{\bm{\theta}+\bm{\epsilon}}(\bm{x}|\mathcal{D})\right]+\frac{\beta\pi}{2}\mathbb{E}_{\bm{\epsilon}}\left[\sigma^{2}_{\bm{\theta}+\bm{\epsilon}}(\bm{x}|\mathcal{D})\right]\right]}{\delta\rho^{2}}\right\rceil.$

We will later simplify this expression using the bounds in (4). For now, we
restrict our focus on the second term in (5). To bound it, we will establish a
bound on
$|\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})-\alpha^{\bm{\theta}+\bm{\epsilon}_{j}}(\bm{x}|\mathcal{D})$.
For a small random perturbation $\bm{\epsilon}_{j}$ we have (with probability
1):

$\displaystyle\alpha^{\bm{\theta}+\bm{\epsilon}_{j}}(\bm{x}|\mathcal{D})=\alpha^{\bm{\theta}}(\bm{x}|\mathcal{D})+\bm{\epsilon}^{\mathsf{T}}_{j}\nabla_{\bm{\theta}}\alpha^{\bm{\theta}}(\bm{x}|\mathcal{D})+o(||\bm{\epsilon_{j}}||)=$
$\displaystyle\alpha^{\bm{\theta}}(\bm{x}|\mathcal{D})+\bm{\epsilon}^{\mathsf{T}}_{j}\nabla_{\bm{\theta}}\left[\mu_{\bm{\theta}}(\bm{x}|\mathcal{D})+\sqrt{\frac{\beta\pi}{2}}\sigma_{\bm{\theta}}(\bm{x}|\mathcal{D})\right]+o(||\bm{\epsilon_{j}}||_{2}).$

Let us define

$\bm{h}_{\bm{\theta}}(\bm{x}|\mathcal{D})=\nabla_{\bm{\theta}}\left[\mu_{\bm{\theta}}(\bm{x}|\mathcal{D})+\sqrt{\frac{\beta\pi}{2}}\sigma_{\bm{\theta}}(\bm{x}|\mathcal{D})\right],$

then, using the Cauchy–Schwarz inequality we have:

$\displaystyle\left|\alpha^{\bm{\theta}+\bm{\epsilon}_{j}}(\bm{x}|\mathcal{D})-\alpha^{\bm{\theta}}(\bm{x}|\mathcal{D})\right|$
$\displaystyle\leq||\bm{\epsilon_{j}}||_{2}||[\bm{h}_{\bm{\theta}}(\bm{x}|\mathcal{D})||_{2}+o(1)]$

Since $\bm{\epsilon}_{j}\sim\mathcal{N}(0,1)$, then with probability at least
$1-\frac{\delta}{4N_{\epsilon}}$:

$\displaystyle||\bm{\epsilon}_{j}||_{2}\leq
4\sigma_{\epsilon}\sqrt{p}+2\sigma_{\epsilon}\sqrt{\log\frac{4N_{\epsilon}}{\delta}}$

Let us assume (and later we will prove the existence of such a bound) that
$||\bm{h}_{\bm{\theta}}(\bm{x}|\mathcal{D})||_{2}\leq A_{1}$. Then, with
probability at least $1-\frac{\delta}{4N_{\epsilon}}$:

$\displaystyle\left|\alpha^{\bm{\theta}+\bm{\epsilon}_{j}}(\bm{x}|\mathcal{D})-\alpha^{\bm{\theta}}(\bm{x}|\mathcal{D})\right|\leq\left[4\sigma_{\epsilon}\sqrt{p}+2\sigma_{\epsilon}\sqrt{\log\frac{4N_{\epsilon}}{\delta}}\right]\left[A_{1}+o(1)\right]$

On the other hand, for
$\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})=\alpha^{\bm{\theta}}(\bm{x}|\mathcal{D})+\eta\sigma_{\eta}$
with probability at least $1-\frac{\delta}{4N_{\epsilon}}$ we have:

$\displaystyle\left|\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})-\alpha^{\bm{\theta}}(\bm{x}|\mathcal{D})\right|\leq\Phi^{-1}\left(1-\frac{\delta}{8N_{\epsilon}}\right)\sigma_{n}.$

where $\Phi(\cdot)$ is the cumulative distribution function for a standard
Gaussian random variable. Hence, by choosing
$\sigma_{\bm{\epsilon}}=\min\left\\{1,\frac{\Phi^{-1}\left(1-\frac{\delta}{8N_{\epsilon}}\right)\sigma_{n}}{\left[4\sqrt{p}+2\sqrt{\log\frac{4N_{\epsilon}}{\delta}}\right]\left[A_{1}+o(1)\right]}\right\\}$
with probability at least $1-\frac{\delta}{2N_{\epsilon}}$ we have that both
$\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})$ and
$\alpha^{\bm{\theta}+\bm{\epsilon}_{j}}(\bm{x}|\mathcal{D})$ belong to the
interval centred at $\alpha^{\bm{\theta}}(\bm{x}|\mathcal{D})$ of size
$\Phi^{-1}\left(1-\frac{\delta}{8N_{\epsilon}}\right)\sigma_{n}$. Therefore,
with probability at least $1-\frac{\delta}{2N_{\epsilon}}$:

$\displaystyle\left|\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})-\alpha^{\bm{\theta}+\bm{\epsilon}_{j}}(\bm{x}|\mathcal{D})\right|\leq
2\Phi^{-1}\left(1-\frac{\delta}{8N_{\epsilon}}\right)\sigma_{n}$

Hence, by choosing
$\sigma_{n}=\frac{\rho}{4\Phi^{-1}\left(1-\frac{\delta}{8N_{\epsilon}}\right)}$
we arrive at:

$\displaystyle\left|\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})-\alpha^{\bm{\theta}+\bm{\epsilon}_{j}}(\bm{x}|\mathcal{D})\right|\leq\frac{\rho}{2}$

and, therefore, for the first term in (5) with probability at least
$1-\frac{\delta}{2}$ we have:

$\displaystyle\frac{1}{N_{\epsilon}}\sum_{j=1}^{N_{\epsilon}}\left|\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})-\alpha^{\bm{\theta}+\bm{\epsilon}_{j}}(\bm{x}|\mathcal{D})\right|\leq\frac{\rho}{2}$

Combining this result with (6) gives, that with probability at least
$1-\delta$:

$\displaystyle\left|\overline{\alpha}^{\bm{\theta}}(\bm{x}|\mathcal{D})-\mathbb{E}_{\epsilon\sim\mathcal{N}(\bm{0},\sigma_{\epsilon}^{2}\bm{I})}\left[\alpha^{\bm{\theta}+\epsilon}(\bm{x}|\mathcal{D})\right]\right|\leq\rho,\
\ \ \ \ \ \ \ \forall\bm{\theta,x},$

upon defining:

$\displaystyle\sigma_{n}=\frac{\rho}{4\Phi^{-1}\left(1-\frac{\delta}{8N_{\epsilon}}\right)},\
\sigma_{\bm{\epsilon}}=\min\left\\{1,\frac{\rho}{8\left[2\sqrt{p}+\sqrt{\log\frac{4N_{\epsilon}}{\delta}}\right]\left[A_{1}+o(1)\right]}\right\\},$
(7)

with

$\displaystyle
N_{\epsilon}=\lceil\frac{16\left[\mathbb{E}_{\bm{\epsilon}}\left[\mu^{2}_{\bm{\theta}+\bm{\epsilon}}(\bm{x}|\mathcal{D})\right]+\frac{\beta\pi}{2}\mathbb{E}_{\bm{\epsilon}}\left[\sigma^{2}_{\bm{\theta}+\bm{\epsilon}}(\bm{x}|\mathcal{D})\right]\right]}{\delta\rho^{2}}\rceil.$

Our last step is to prove the existence of a constant $A_{1}$ such that
$||\bm{h}_{\bm{\theta}}(\bm{x})||_{2}\leq A_{1}$ and also to simplify these
expressions by deriving bounds on
$\mathbb{E}_{\bm{\epsilon}}\left[\mu_{\bm{\theta}+\bm{\epsilon}}(\bm{x}|\mathcal{D})\right]$
and
$\mathbb{E}_{\bm{\epsilon}}\left[\sigma^{2}_{\bm{\theta}+\bm{\epsilon}}(\bm{x}|\mathcal{D})\right]$.
This will be provided as a separate Claim:

Claim: Let the bounds in (4) hold, then there are positive constants
$A_{1},A_{2}$ and $A_{3}$, such that

$||\bm{h}_{\bm{\theta}}(\bm{x})||_{2}\leq A_{1},\ \
\mathbb{E}_{\bm{\epsilon}}\left[\mu_{\bm{\theta}+\bm{\epsilon}}(\bm{x}|\mathcal{D})\right]\leq
A_{2},\ \ \
\mathbb{E}_{\bm{\epsilon}}\left[\sigma^{2}_{\bm{\theta}+\bm{\epsilon}}(\bm{x}|\mathcal{D})\right]\leq
A_{3}.$ (8)

###### Proof.

We start with the bound on $||\bm{h}_{\bm{\theta}}(\bm{x})||_{2}$. Let us
denote for simplicity
$\bm{a}_{\bm{\theta}}=[k_{\bm{\theta}}(\bm{x},\bm{x}_{i})]_{\bm{x}_{i}\in\mathcal{D}}$,
$\bm{B}_{\bm{\theta}}=\left[\left[k_{\bm{\theta}}(\bm{x},\bm{x}^{\prime})\right]_{\bm{x}\in\mathcal{D},\bm{x}^{\prime}\in\mathcal{D}}+\bm{I}\right]^{-1}$,
$\bm{y}=[y(\bm{x})]_{\bm{x\in\mathcal{D}}}$,
$\bm{m}_{\mathcal{D}}=[m(\bm{x})]_{\bm{x}\in\mathcal{D}}$, $m=m(\bm{x})$, and
$k_{\bm{\theta}}=k_{\bm{\theta}}(\bm{x},\bm{x})$, then

$\displaystyle\mu_{\bm{\theta}}(\bm{x}|\mathcal{D})=\bm{a}^{\mathsf{T}}_{\bm{\theta}}\bm{B}_{\bm{\theta}}[\bm{y}-\bm{m}_{\mathcal{D}}]+m,\
\
\sigma^{2}_{\bm{\theta}}(\bm{x}|\mathcal{D})=\bm{a}_{\bm{\theta}}^{\mathsf{T}}\bm{B}_{\bm{\theta}}\bm{a}_{\bm{\theta}}+k_{\bm{\theta}}$

Let us also denote the size of $\mathcal{D}$ as $N$, then we have:

$\displaystyle\nabla_{\boldsymbol{\theta}}\mu_{\bm{\theta}}(\bm{x}|\mathcal{D})=\sum_{i=1}^{N}\sum_{j=1}^{N}\nabla_{\bm{\theta}}\left[[y_{j}-m_{j}][\bm{a}_{\bm{\theta}}]_{i}[\boldsymbol{B}_{\bm{\theta}}]_{ij}\right]=$
$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}\left[[y_{j}-m_{j}][\bm{B}_{\bm{\theta}}]_{ij}\nabla_{\bm{\theta}}[\bm{a}_{\bm{\theta}}]_{i}\right]+\sum_{i=1}^{N}\sum_{j=1}^{N}\left[[y_{j}-m_{j}][\bm{a}_{\bm{\theta}}]_{i}\nabla_{\bm{\theta}}[[\bm{B}_{\bm{\theta}}]_{ij}]\right].$

Consider each term in this expression separately:

$\displaystyle\left|\left|\sum_{i=1}^{N}\sum_{j=1}^{N}[y_{j}-m_{j}][\bm{B}_{\bm{\theta}}]_{ij}\nabla_{\bm{\theta}}[\bm{a}_{\bm{\theta}}]_{i}\right|\right|_{2}=\left|\left|\sum_{i=1}^{N}[\bm{B}_{\bm{\theta}}[\bm{y}-\bm{m}_{\mathcal{D}}]]_{i}\nabla_{\bm{\theta}}[\bm{a}_{\bm{\theta}}]_{i}\right|\right|_{2}\leq$
$\displaystyle\sum_{i=1}^{N}||\bm{B}_{\bm{\theta}}(i,:)||_{2}||\bm{y}-\bm{m}_{\mathcal{D}}||_{2}\left|\left|\nabla_{\bm{\theta}}[\bm{a}_{\bm{\theta}}]_{i}\right|\right|_{2}.$

Using $|y|\leq C$ and $|m(\boldsymbol{x})|\leq M_{4}$, we have
$||\bm{y}-\bm{m}_{\mathcal{D}}||_{2}\leq(C+M_{4})\sqrt{N}$ and
$||\nabla_{\boldsymbol{\theta}}[\boldsymbol{a}(\boldsymbol{\theta})]_{i}||_{2}\leq
M_{2}$ and as such:

$\displaystyle\left|\left|\sum_{i=1}^{N}\sum_{j=1}^{N}[y_{j}-m_{j}][\bm{B}_{\bm{\theta}}]_{ij}\nabla_{\bm{\theta}}[\bm{a}_{\bm{\theta}}]_{i}\right|\right|_{2}\leq(C+M_{4})\sqrt{N}\sum_{i=1}^{N}||\bm{B}(\bm{\theta})(i,:)||_{2}\left|\left|\nabla_{\bm{\theta}}[\bm{a}_{\bm{\theta}}]_{i}\right|\right|_{2}\leq$
(9)
$\displaystyle(C+M_{4})\sqrt{N}||\bm{B}_{\bm{\theta}}||_{F}\sum_{i=1}^{N}||\nabla_{\bm{\theta}}[\boldsymbol{a}_{\bm{\theta}}]_{i}||_{2}\leq(C+M_{4})\sqrt{N}\sqrt{\text{rank}(\bm{B}_{\bm{\theta}}}||\bm{B}_{\bm{\theta}}||_{2}\sum_{i=1}^{N}||\nabla_{\bm{\theta}}[\bm{a}_{\bm{\theta}}]_{i}||_{2}\leq$
$\displaystyle\frac{(C+M_{4})N}{\sigma^{2}_{n}}\sum_{i=1}^{N}||\nabla_{\boldsymbol{\theta}}[\boldsymbol{a}(\boldsymbol{\theta})]_{i}||_{2}=\frac{(C+M_{4})N^{2}M_{2}}{\sigma^{2}_{n}}$

Now, let us consider the second term in the expression for the posterior mean:

$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}[y_{j}-m_{j}][\bm{a}_{\bm{\theta}}]_{i}\nabla_{\bm{\theta}}\left[[\bm{B}_{\bm{\theta}}]_{ij}\right]=\sum_{i=1}^{N}[\bm{a}_{\bm{\theta}}]_{i}\left[\sum_{j=1}[y_{j}-m_{j}]\nabla_{\boldsymbol{\theta}}\left[[\bm{B}_{\bm{\theta}}]_{ij}\right]\right].$

Notice, that the gradient expression above is presented in the form of a
vector:

$\displaystyle\nabla_{\bm{\theta}}\left[[\bm{B}_{\bm{\theta}}]_{ij}\right]=\left[\begin{array}[]{c}\frac{\partial}{\partial\theta_{1}}\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]^{-1}_{ij},\\\
\vdots\\\
\frac{\partial}{\partial\theta_{p}}\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]^{-1}_{ij}\\\
\end{array},\right]$

where we use the notation
$\bm{K}_{\bm{\theta}}=[k_{\bm{\theta}}(\bm{x}_{i},\bm{x}_{j})]^{N,N}_{i=1,j=1}$.
For the $r^{th}$ component we have

$\frac{\partial}{\partial\theta_{r}}\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]^{-1}_{ij}=\left[-\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]^{-1}\frac{\partial}{\partial\theta_{r}}\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]^{-1}\right]_{ij}$
(10)

Now we can study the gradient of the second term in the posterior mean
expression,

$\displaystyle\left|\left|\sum_{i=1}^{N}[\bm{a}_{\bm{\theta}}]_{i}\left[\sum_{j=1}[y_{j}-m_{j}]\nabla_{\bm{\theta}}\left[[\bm{B}_{\bm{\theta}}]_{ij}\right]\right]\right|\right|_{2}\leq\sum_{i=1}^{N}\left|[\bm{a}_{\bm{\theta}}]_{i}\right|\left[\sum_{j=1}^{N}||\bm{y}-\bm{m}_{\mathcal{D}}||_{2}\left|\left|\nabla_{\bm{\theta}}\left[[\bm{B}_{\bm{\theta}}]_{ij}\right]\right|\right|_{2}\right]\leq$
$\displaystyle(C+M_{4})\sqrt{N}M_{1}\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{r=1}^{p}\left|\frac{\partial}{\partial\theta_{r}}\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]^{-1}_{ij}\right|.$

Using result (10) in the above expression we now have

$\displaystyle\left|\left|\sum_{i=1}^{N}[\bm{a}_{\bm{\theta}}]_{i}\left[\sum_{j=1}[y_{j}-m_{j}]\nabla_{\bm{\theta}}\left[[\bm{B}_{\bm{\theta}}]_{ij}\right]\right]\right|\right|_{2}\leq$
$\displaystyle(C+M_{4})\sqrt{N}M_{1}\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{r=1}^{p}\left|\frac{\partial}{\partial\theta_{r}}\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]^{-1}_{ij}\right|\leq$
$\displaystyle(C+M_{4})\sqrt{N}M_{1}\sum_{r=1}^{p}\sum_{i=1}^{N}\sum_{j=1}^{N}\left|\frac{\partial}{\partial\theta_{r}}\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]^{-1}_{ij}\right|\leq$
$\displaystyle(C+M_{4})N\sqrt{N}M_{1}\times\sum_{r=1}^{p}\left|\left|\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]^{-1}\frac{\partial}{\partial\theta_{r}}\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]^{-1}\right|\right|_{F},$

where we used that $\sum_{i=1}^{N}\sum_{j=1}^{N}|\boldsymbol{C}_{ij}|\leq
N||\boldsymbol{C}||_{F}$ for any arbitrary matrix
$\boldsymbol{C}\in\mathbb{R}^{N\times N}$. Because
$\frac{\partial}{\partial\theta_{r}}[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}]=\frac{\partial}{\partial\theta_{r}}\bm{K}_{\boldsymbol{\theta}}$

$\displaystyle\frac{\left|\left|\sum_{i=1}^{N}[\bm{a}_{\bm{\theta}}]_{i}\left[\sum_{j=1}[y_{j}-m_{j}]\nabla_{\bm{\theta}}\left[[\bm{B}_{\bm{\theta}}]_{ij}\right]\right]\right|\right|_{2}}{(C+M_{4})N\sqrt{N}M_{1}}\leq$
$\displaystyle\sum_{r=1}^{p}\left|\left|\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]^{-1}\frac{\partial}{\partial\theta_{r}}\bm{K}_{\boldsymbol{\theta}}\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]^{-1}\right|\right|_{F}\leq$
$\displaystyle\sqrt{N}\sum_{r=1}^{p}\left|\left|\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]^{-1}\frac{\partial}{\partial\theta_{r}}\bm{K}_{\boldsymbol{\theta}}\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]^{-1}\right|\right|_{2}$

and using the properties of the matrix 2-norm $||\cdot||_{2}$

$\displaystyle\left|\left|\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]^{-1}\right|\right|_{2}\leq\frac{1}{\sigma^{2}_{n}}$

Hence,

$\displaystyle\frac{\left|\left|\sum_{i=1}^{N}[\bm{a}_{\bm{\theta}}]_{i}\left[\sum_{j=1}[y_{j}-m_{j}]\nabla_{\bm{\theta}}\left[[\bm{B}_{\bm{\theta}}]_{ij}\right]\right]\right|\right|_{2}}{(C+M_{4})N^{2}M_{1}}\leq$
(11)
$\displaystyle\sum_{r=1}^{p}\left|\left|\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]^{-1}\right|\right|_{2}\left|\left|\frac{\partial}{\partial\theta_{r}}\bm{K}_{\boldsymbol{\theta}}\right|\right|_{2}\left|\left|\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]^{-1}\right|\right|_{2}\leq\frac{1}{\sigma^{4}_{n}}\sum_{r=1}^{p}\left|\left|\frac{\partial}{\partial\theta_{r}}\bm{K}_{\boldsymbol{\theta}}\right|\right|_{2}.$

Let us study the last term in the expression above. Using
$\sqrt{c^{2}_{1}+\ldots+c^{2}_{R}}\leq|c_{1}|+\ldots+|c_{R}|$ for any set of
real numbers $c_{1},\ldots,c_{R}\in\mathbb{R}$ we have

$\displaystyle\left|\left|\frac{\partial}{\partial\theta_{r}}\bm{K}_{\boldsymbol{\theta}}\right|\right|_{2}=\left|\left|\left[\begin{array}[]{ccc}\frac{\partial}{\partial\theta_{r}}k_{\boldsymbol{\theta}}(\boldsymbol{x}_{1},\boldsymbol{x}_{1}),&\ldots&\frac{\partial}{\partial\theta_{r}}k_{\boldsymbol{\theta}}(\boldsymbol{x}_{1},\boldsymbol{x}_{N})\\\
\vdots&\vdots&\vdots\\\
\frac{\partial}{\partial\theta_{r}}k_{\boldsymbol{\theta}}(\boldsymbol{x}_{N},\boldsymbol{x}_{1}),&\ldots&\frac{\partial}{\partial\theta_{r}}k_{\boldsymbol{\theta}}(\boldsymbol{x}_{N},\boldsymbol{x}_{N})\\\
\end{array}l\right]\right|\right|_{2}\leq$
$\displaystyle\left|\left|\left[\begin{array}[]{ccc}\frac{\partial}{\partial\theta_{r}}k_{\boldsymbol{\theta}}(\boldsymbol{x}_{1},\boldsymbol{x}_{1}),&\ldots&\frac{\partial}{\partial\theta_{r}}k_{\boldsymbol{\theta}}(\boldsymbol{x}_{1},\boldsymbol{x}_{N})\\\
\vdots&\vdots&\vdots\\\
\frac{\partial}{\partial\theta_{r}}k_{\boldsymbol{\theta}}(\boldsymbol{x}_{N},\boldsymbol{x}_{1}),&\ldots&\frac{\partial}{\partial\theta_{r}}k_{\boldsymbol{\theta}}(\boldsymbol{x}_{N},\boldsymbol{x}_{N})\\\
\end{array}\right]\right|\right|_{F}=$
$\displaystyle\sqrt{\sum_{i=1}^{N}\sum_{j=1}^{N}\left[\frac{\partial}{\partial\theta_{r}}k_{\boldsymbol{\theta}}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})\right]^{2}}\leq\sum_{i=1}^{N}\sum_{j=1}^{N}\left|\frac{\partial}{\partial\theta_{r}}k_{\boldsymbol{\theta}}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})\right|.$

Substituting this expression in (11) gives us

$\displaystyle\frac{\left|\left|\sum_{i=1}^{N}[\bm{a}_{\bm{\theta}}]_{i}\left[\sum_{j=1}[y_{j}-m_{j}]\nabla_{\bm{\theta}}\left[[\bm{B}_{\bm{\theta}}]_{ij}\right]\right]\right|\right|_{2}}{(C+M_{4})N^{2}M_{1}}\leq\frac{1}{\sigma^{4}_{n}}\sum_{r=1}^{d}\sum_{i=1}^{N}\sum_{j=1}^{N}\left|\frac{\partial}{\partial\theta_{r}}k_{\boldsymbol{\theta}}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})\right|\leq$
(12)
$\displaystyle\frac{\sqrt{p}}{\sigma^{4}_{n}}\sum_{i=1}^{N}\sum_{j=1}^{N}\left|\left|\nabla_{\boldsymbol{\theta}}k_{\boldsymbol{\theta}}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})\right|\right|_{2}\leq\frac{N^{2}\sqrt{p}M_{2}}{\sigma^{4}_{n}}.$

Hence, combining results (9) and (12) we have

$||\nabla_{\boldsymbol{\theta}}\mu_{\bm{\theta}}(\bm{x}|\mathcal{D})||_{2}\leq\frac{(C+M_{4})N^{2}M_{2}}{\sigma^{2}_{n}}\left[1+\frac{N^{2}M_{1}\sqrt{p}}{\sigma^{2}_{n}}\right].$
(13)

Now, let us focus on the gradient of the posterior standard deviation,

$\displaystyle\nabla_{\bm{\theta}}\sigma_{\bm{\theta}}(\bm{x}|\mathcal{D})=\nabla_{\bm{\theta}}\left[\sqrt{k_{\bm{\theta}}(\boldsymbol{x},\bm{x})-\bm{a}_{\bm{\theta}}^{\mathsf{T}}[\bm{K}_{\bm{\theta}}+\sigma^{2}_{n}\textbf{I}]^{-1}\bm{a}_{\bm{\theta}}}\right]=$
$\displaystyle\frac{1}{2\sigma_{\bm{\theta}}(\bm{x}|\mathcal{D})}\nabla_{\bm{\theta}}\left[k_{\bm{\theta}}(\bm{x},\boldsymbol{x})-\bm{a}_{\bm{\theta}}^{\mathsf{T}}[\bm{K}_{\bm{\theta}}+\sigma^{2}_{n}\textbf{I}]^{-1}\bm{a}_{\bm{\theta}}\right]=\frac{1}{2\sigma_{\bm{\theta}}(\bm{x})}\left[\nabla_{\bm{\theta}}k_{\bm{\theta}}(\bm{x},\bm{x})-\nabla_{\bm{\theta}}\left[\bm{a}_{\bm{\theta}}^{\mathsf{T}}[\bm{K}_{\bm{\theta}}+\sigma^{2}_{n}\textbf{I}]^{-1}\bm{a}_{\bm{\theta}}\right]\right].$

Let us study the second gradient expression. Using our notation we have

$\displaystyle\bm{a}_{\bm{\theta}}^{\mathsf{T}}[\bm{K}_{\bm{\theta}}+\sigma^{2}_{n}\textbf{I}]^{-1}\bm{a}_{\bm{\theta}}=\bm{a}^{\mathsf{T}}_{\bm{\theta}}\bm{B}_{\bm{\theta}}\bm{a}_{\bm{\theta}}=\sum_{i=1}^{N}\sum_{j=1}^{N}[\bm{a}_{\bm{\theta}}]_{i}[\bm{a}_{\bm{\theta}}]_{j}\left[\bm{B}_{\bm{\theta}}\right]_{ij}.$

Hence, for the gradient,

$\displaystyle\nabla_{\bm{\theta}}\left[\bm{a}^{\mathsf{T}}_{\bm{\theta}}\bm{B}_{\bm{\theta}}\bm{a}_{\bm{\theta}}\right]=\sum_{i=1}^{N}\sum_{j=1}^{N}\nabla_{\bm{\theta}}\left[[\bm{a}_{\bm{\theta}}]_{i}[\bm{a}_{\bm{\theta}}]_{j}\left[\bm{B}_{\bm{\theta}}\right]_{ij}\right]=$
$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}\nabla_{\bm{\theta}}\left[[\bm{a}_{\bm{\theta}}]_{i}\right][\bm{a}_{\bm{\theta}}]_{j}\left[\bm{B}_{\bm{\theta}}\right]_{ij}+\sum_{i=1}^{N}\sum_{j=1}^{N}\nabla_{\bm{\theta}}\left[[\bm{a}_{\bm{\theta}}]_{j}\right][\bm{a}_{\bm{\theta}}]_{i}\left[\bm{B}_{\bm{\theta}}\right]_{ij}+$
$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}\nabla_{\bm{\theta}}\left[\left[\bm{B}_{\bm{\theta}}\right]_{ij}\right][\bm{a}_{\bm{\theta}}]_{i}[\bm{a}_{\bm{\theta}}]_{j}.$

and for the norm of the above expression we have

$\displaystyle\left|\left|\nabla_{\bm{\theta}}\left[\bm{a}^{\mathsf{T}}_{\bm{\theta}}\bm{B}_{\bm{\theta}}\bm{a}_{\bm{\theta}}\right]\right|\right|_{2}=\sum_{i=1}^{N}\sum_{j=1}^{N}\nabla_{\bm{\theta}}\left[[\bm{a}_{\bm{\theta}}]_{i}[\bm{a}_{\bm{\theta}}]_{j}\left[\bm{B}_{\bm{\theta}}\right]_{ij}\right]=$
$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}\left|\left|\nabla_{\bm{\theta}}\left[\left[\bm{B}_{\bm{\theta}}\right]_{ij}\right]\right|\right|_{2}\left|[\bm{a}_{\bm{\theta}}]_{i}[\bm{a}_{\bm{\theta}}]_{j}\right|+\sum_{i=1}^{N}\sum_{j=1}^{N}\left|\left|\nabla_{\bm{\theta}}\right|\right|_{2}\left|\left[[\bm{a}_{\bm{\theta}}]_{i}\right][\bm{a}_{\bm{\theta}}]_{j}\left[\bm{B}_{\bm{\theta}}\right]_{ij}\right|+$
$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}\left|\left|\nabla_{\bm{\theta}}\left[[\bm{a}_{\bm{\theta}}]_{j}\right]\right|\right|_{2}\left|[\bm{a}_{\bm{\theta}}]_{i}\left[\bm{B}_{\bm{\theta}}\right]_{ij}\right|.$

Let us now bound each term in this expression:

  1. 1.

The first term:

$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}\left|\left|\nabla_{\bm{\theta}}\left[\left[\bm{B}_{\bm{\theta}}\right]_{ij}\right]\right|\right|_{2}\left|[\bm{a}_{\bm{\theta}}]_{i}[\bm{a}_{\bm{\theta}}]_{j}\right|\leq$
$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}\left|\left|\nabla_{\bm{\theta}}\left[\left[\bm{B}_{\bm{\theta}}\right]_{ij}\right]\right|\right|_{2}\left|\left|\bm{a}_{\bm{\theta}}\right|\right|_{2}\left|\left|\bm{a}_{\bm{\theta}}\right|\right|_{2}\leq
M^{2}_{1}\sum_{i=1}^{N}\sum_{j=1}^{N}\left|\left|\nabla_{\bm{\theta}}\left[\left[\bm{B}_{\bm{\theta}}\right]_{ij}\right]\right|\right|_{2}$

Using the previous bound for
$\left|\left|\nabla_{\bm{\theta}}\left[\left[\bm{B}_{\bm{\theta}}\right]_{ij}\right]\right|\right|_{2}$
we have:

$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}\left|\left|\nabla_{\bm{\theta}}\left[\left[\bm{B}_{\bm{\theta}}\right]_{ij}\right]\right|\right|_{2}\left|[\bm{a}_{\bm{\theta}}]_{i}[\bm{a}_{\bm{\theta}}]_{j}\right|\leq
M^{2}_{1}\sum_{i=1}^{N}\sum_{j=1}^{N}\sum_{r=1}^{p}\left|\frac{\partial}{\partial\theta_{r}}\left[\bm{K}_{\bm{\theta}}+\sigma^{2}_{n}\textbf{I}\right]^{-1}_{ij}\right|=$
$\displaystyle
NM^{2}_{1}\sum_{r=1}^{p}\left|\left|\left[\bm{K}_{\bm{\theta}}+\sigma_{n}\bm{I}\right]^{-1}\frac{\partial}{\partial\theta_{r}}\bm{K}_{\bm{\theta}}\left[\bm{K}_{\bm{\theta}}+\sigma_{n}\bm{I}\right]^{-1}\right|\right|_{F}\leq$
$\displaystyle
N^{\frac{3}{2}}M^{2}_{1}\sum_{r=1}^{p}\left|\left|\left[\bm{K}_{\bm{\theta}}+\sigma_{n}\bm{I}\right]^{-1}\frac{\partial}{\partial\theta_{r}}\bm{K}_{\boldsymbol{\theta}}\left[\bm{K}_{\bm{\theta}}+\sigma_{n}\bm{I}\right]^{-1}\right|\right|_{2}$

Since
$\left|\left|\left[\bm{K}_{\boldsymbol{\theta}}+\sigma_{n}\boldsymbol{I}\right]^{-1}\right|\right|_{2}\leq\frac{1}{\sigma^{2}_{n}}$
we have:

$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}\left|\left|\nabla_{\bm{\theta}}\left[\left[\bm{B}_{\bm{\theta}}\right]_{ij}\right]\right|\right|_{2}\left|[\bm{a}_{\bm{\theta}}]_{i}[\bm{a}_{\bm{\theta}}]_{j}\right|\leq\frac{N\sqrt{N}M^{2}_{1}}{\sigma^{4}_{n}}\sum_{r=1}^{p}\sum_{i=1}^{N}\sum_{j=1}^{N}\left|\frac{\partial}{\partial\theta_{r}}k_{\boldsymbol{\theta}}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})\right|.$

Using
$\sum_{r=1}^{p}\sum_{i=1}^{N}\sum_{j=1}^{N}\left|\frac{\partial}{\partial\theta_{r}}k_{\bm{\theta}}(\bm{x}_{i},\bm{x}_{j})\right|=\sqrt{p}\sum_{i=1}^{N}\sum_{j=1}^{N}\left|\left|\nabla_{\boldsymbol{\theta}}k_{\boldsymbol{\theta}}(\boldsymbol{x}_{i},\boldsymbol{x}_{j})\right|\right|_{2}\leq
N^{2}\sqrt{p}M_{2}$, we have:

$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}\left|\left|\nabla_{\bm{\theta}}\left[\left[\bm{B}_{\bm{\theta}}\right]_{ij}\right]\right|\right|_{2}\left|[\bm{a}_{\bm{\theta}}]_{i}[\bm{a}_{\bm{\theta}}]_{j}\right|\leq\frac{N^{\frac{7}{2}}\sqrt{p}M^{2}_{1}M_{2}}{\sigma^{4}_{n}}$

  2. 2.

The second and the third terms are identical with respect to the bounding
strategy,

$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}||\nabla_{\bm{\theta}}\left[[\bm{a}_{\bm{\theta}}]_{i}\right]||_{2}\left|[\bm{a}_{\bm{\theta}}]_{j}\left[\bm{B}_{\bm{\theta}}\right]_{ij}\right|=\sum_{i=1}^{N}\left|\bm{B}_{\bm{\theta}}(i,:)\bm{a}_{\bm{\theta}}\right|\left|\left|\nabla_{\bm{\theta}}\left[[\bm{a}_{\bm{\theta}}]_{i}\right]\right|\right|_{2}\leq$
$\displaystyle\sum_{i=1}^{N}\left|\left|\bm{B}_{\bm{\theta}}(i,:)\right|\right|_{2}\left|\left|\bm{a}_{\bm{\theta}}\right|\right|_{2}\left|\left|\nabla_{\bm{\theta}}\left[[\bm{a}_{\bm{\theta}}]_{i}\right]\right|\right|_{2}\leq||\bm{B}_{\bm{\theta}}||_{F}||\bm{a}_{\bm{\theta}}||_{2}\sum_{i=1}^{N}\left|\left|\nabla_{\bm{\theta}}\left[[\bm{a}_{\bm{\theta}}]_{i}\right]\right|\right|_{2},$

since
$||\bm{B}_{\bm{\theta}}||_{F}\leq\sqrt{\text{rank}(\bm{B}_{\bm{\theta}})}||\bm{B}_{\bm{\theta}}||_{2}\leq\frac{\sqrt{N}}{\sigma^{2}_{n}}$.
Hence,

$\displaystyle\sum_{i=1}^{N}\sum_{j=1}^{N}||\nabla_{\bm{\theta}}\left[[\bm{a}_{\bm{\theta}}]_{i}\right]||_{2}\left|[\bm{a}_{\bm{\theta}}]_{j}\left[\bm{B}_{\bm{\theta}}\right]_{ij}\right|\leq\frac{N\sqrt{N}M_{1}M_{2}}{\sigma^{2}_{n}}.$

Combining these results and using
$||\nabla_{\bm{\theta}}k_{\bm{\theta}}(\bm{x},\bm{x})||\leq M_{2}$,
$\left|\sigma_{\bm{\theta}}(\bm{x}|\mathcal{D})\right|\geq
k_{\bm{\theta}}(\bm{x},\bm{x})\geq M_{0}$, we have

$\left|\left|\nabla_{\bm{\theta}}\left[\sigma_{\bm{\theta}}(\bm{x}|\mathcal{D})\right]\right|\right|_{2}\leq\frac{N\sqrt{N}M_{1}M_{2}}{2\sigma^{2}_{n}M_{0}}\left[\frac{N^{2}\sqrt{p}M_{1}}{\sigma^{2}_{n}}+2\right]$
(14)

Hence, combining (13) and (14) we have

$\displaystyle||\bm{h}_{\bm{\theta}}(\bm{x}|\mathcal{D})||_{2}\leq||\nabla_{\boldsymbol{\theta}}\mu_{\bm{\theta}}(\bm{x}|\mathcal{D})||_{2}+\sqrt{\frac{\beta\pi}{2}}\left|\left|\nabla_{\bm{\theta}}\left[\sigma_{\bm{\theta}}(\bm{x}|\mathcal{D})\right]\right|\right|_{2}\leq$
$\displaystyle\frac{(C+M_{4})N^{2}M_{2}}{\sigma^{2}_{n}}\left[1+\frac{N^{2}M_{1}\sqrt{p}}{\sigma^{2}_{n}}\right]+\sqrt{\frac{\beta\pi}{2}}\frac{N\sqrt{N}M_{1}M_{2}}{2\sigma^{2}_{n}M_{0}}\left[\frac{N^{2}\sqrt{p}M_{1}}{\sigma^{2}_{n}}+2\right]\triangleq
A_{1}.$

Now, we are ready to bound the other two terms in the claim:

$\displaystyle\mu^{2}_{\bm{\theta}+\bm{\epsilon}}(\bm{x}|\mathcal{D})\leq
2\left[\bm{a}^{\mathsf{T}}_{\bm{\theta}+\bm{\epsilon}}\bm{B}_{\bm{\theta}+\bm{\epsilon}}[\bm{y}-\bm{m}_{\mathcal{D}}]\right]^{2}+2|m|^{2}\leq
2\frac{(C+M_{4})^{2}M^{2}_{1}}{\sigma^{4}_{n}}+2M^{2}_{4}$

Therefore, for
$\mathbb{E}_{\bm{\epsilon}}\left[\mu^{2}_{\bm{\theta}+\bm{\epsilon}}(\bm{x}|\mathcal{D})\right]$
we have

$\displaystyle\mathbb{E}_{\bm{\epsilon}}\left[\mu^{2}_{\bm{\theta}+\bm{\epsilon}}(\bm{x}|\mathcal{D})\right]\leq
2\frac{(C+M_{4})^{2}M^{2}_{1}}{\sigma^{4}_{n}}+2M^{2}_{4}\triangleq A_{2}.$

Finally, for the posterior mean

$\displaystyle\sigma^{2}_{\bm{\theta}+\bm{\epsilon}}(\bm{x}|\mathcal{D})\leq
k_{\bm{\theta}}(\bm{x},\bm{x})+\bm{a}^{\mathsf{T}}_{\bm{\theta}+\bm{\epsilon}}\bm{B}_{\bm{\theta}+\bm{\epsilon}}\bm{a}_{\bm{\theta}+\bm{\epsilon}}\leq
M_{1}+\frac{M^{2}_{1}}{\sigma^{2}_{n}}$

Therefore, for
$\mathbb{E}_{\bm{\epsilon}}\left[\sigma^{2}_{\bm{\theta}+\bm{\epsilon}}(\bm{x}|\mathcal{D})\right]$
we have

$\displaystyle\mathbb{E}_{\bm{\epsilon}}\left[\sigma^{2}_{\bm{\theta}+\bm{\epsilon}}(\bm{x}|\mathcal{D})\right]\leq
M_{1}+\frac{M^{2}_{1}}{\sigma^{2}_{n}}\triangleq A_{3}.$

This finishes the proof of the claim. ∎

Equipped with these results, we can further simplify the expressions (7):

$\displaystyle\sigma_{n}=\frac{\rho}{4\Phi^{-1}\left(1-\frac{\delta}{8N_{\epsilon}}\right)},\
\
\sigma_{\bm{\epsilon}}=\min\left\\{1,\frac{\rho}{8\left[2\sqrt{p}+\sqrt{\log\frac{4N_{\epsilon}}{\delta}}\right]\left[A_{1}+o(1)\right]}\right\\},$

with

$\displaystyle
N_{\epsilon}=\left\lceil\frac{16\left[A_{2}+\frac{\beta\pi}{2}A_{3}\right]}{\delta\rho^{2}}\right\rceil.$

This finishes the proof of the lemma.

As such, we may now implement robust formulations of acquisition functions
using only the GP predictive mean and variance.

##  C Statistical Hypothesis Tests for Heteroscedasticity

In this section we present the full results for the statistical hypothesis
testing using Levene’s test and the Fligner-Killeen test in Table 5, Table 6,
Table 7, Table 8, Table 9 and Table 10 for the Boston, breast cancer,
diabetes, digits, iris and wine datasets respectively.

Table 5: Heteroscedasticity tests on tasks involving Boston data set. Data set | Model | Metric | Fligner Statistic | p-value | Levene Statistic | p-value
---|---|---|---|---|---|---
Boston | DT | mae | 73.51 | 0.01327 | 1.752 | 1.900e-03
| MLP-adam | mae | 336.3 | 1.737e-44 | 14.4 | 3.611e-65
| MLP-SGD | mae | 272.6 | 8.694e-33 | 6.561 | 1.480e-29
| RF | mae | 28.79 | 0.9906 | 0.6768 | 0.9537
| SVM | mae | 48.08 | 0.5106 | 0.9612 | 0.5508
| ada | mae | 218.7 | 2.692e-23 | 13.59 | 5.542e-62
| kNN | mae | 33.15 | 0.9597 | 0.619 | 0.98
| lasso | mae | 30.4 | 0.983 | 0.6091 | 0.983
| linear | mae | 16.17 | 1 | 0.251 | 1
| DT | mse | 60.75 | 0.1211 | 1.33 | 0.07387
| MLP-adam | mse | 387 | 4.504e-54 | 15.32 | 1.147e-68
| MLP-SGD | mse | 353.2 | 1.185e-47 | 8.239 | 3.548e-38
| RF | mse | 35.59 | 0.9242 | 0.8985 | 0.6692
| SVM | mse | 25.01 | 0.9983 | 0.4491 | 0.9996
| ada | mse | 249.1 | 1.398e-28 | 14.4 | 3.682e-65
| kNN | mse | 27.75 | 0.9938 | 0.8247 | 0.7951
| lasso | mse | 31.38 | 0.9764 | 0.5397 | 0.9955
| linear | mse | 16.67 | 1 | 0.1726 | 1
Table 6: Heteroscedasticity tests on tasks involving Breast cancer (BC) data set. Data set | Model | Metric | Fligner Statistic | p-value | Levene Statistic | p-value
---|---|---|---|---|---|---
BC | DT | acc | 97.79 | 4.302e-05 | 4.62 | 6.650e-19
| MLP-adam | acc | 133 | 1.113e-09 | 2.939 | 1.923e-09
| MLP-SGD | acc | 116.8 | 1.854e-07 | 2.469 | 6.495e-07
| RF | acc | 154.9 | 6.469e-13 | 6.661 | 4.353e-30
| SVM | acc | 20.7 | 0.9999 | 0.3995 | 0.9999
| ada | acc | 272.5 | 9.178e-33 | 13.57 | 6.582e-62
| kNN | acc | 33.16 | 0.9596 | 0.5519 | 0.9941
| lasso | acc | 20.78 | 0.9999 | 0.4291 | 0.9998
| linear | acc | 21.15 | 0.9998 | 0.4545 | 0.9995
| DT | nll | 260.5 | 1.280e-30 | 9.52 | 2.294e-44
| MLP-adam | nll | 166.6 | 1.008e-14 | 3.643 | 2.247e-13
| MLP-SGD | nll | 141.2 | 7.115e-11 | 2.669 | 5.661e-08
| RF | nll | 185.8 | 8.495e-18 | 7.553 | 1.013e-34
| SVM | nll | 76.98 | 6.526e-03 | 1.707 | 2.970e-03
| ada | nll | 142 | 5.458e-11 | 4.283 | 5.274e-17
| kNN | nll | 125.7 | 1.155e-08 | 4.337 | 2.635e-17
| lasso | nll | 71.41 | 0.02 | 1.011 | 0.4565
| linear | nll | 18.55 | 1 | 0.2714 | 1
Table 7: Heteroscedasticity tests on tasks involving diabetes data set. Dataset | Model | Metric | Fligner Statistic | p-value | Levene Statistic | p-value
---|---|---|---|---|---|---
Diabetes | DT | mae | 56.52 | 0.2146 | 1.131 | 0.2601
| MLP-adam | mae | 74.64 | 0.01059 | 2.573 | 1.847e-07
| MLP-SGD | mae | 191.3 | 1.062e-18 | 17.87 | 8.498e-78
| RF | mae | 79.38 | 3.898e-03 | 1.558 | 0.01174
| SVM | mae | 2.436 | 1 | 1.810e-04 | 1
| ada | mae | 179.8 | 7.883e-17 | 7.542 | 1.154e-34
| kNN | mae | 67.48 | 0.04106 | 2.101 | 4.747e-05
| lasso | mae | 176.2 | 2.950e-16 | 4.75 | 1.225e-19
| linear | mae | 206 | 3.792e-21 | 5.714 | 5.490e-25
| DT | mse | 44.52 | 0.6551 | 0.8264 | 0.7925
| MLP-adam | mse | 100.4 | 2.109e-05 | 3.582 | 4.951e-13
| MLP-SGD | mse | 202.9 | 1.257e-20 | 14.31 | 7.960e-65
| RF | mse | 37.1 | 0.8938 | 0.8063 | 0.8224
| SVM | mse | 4.004 | 1 | 4.740e-04 | 1
| ada | mse | 189 | 2.510e-18 | 7.348 | 1.138e-33
| kNN | mse | 88.62 | 4.545e-04 | 2.964 | 1.407e-09
| lasso | mse | 257.6 | 4.341e-30 | 10.86 | 1.637e-50
| linear | mse | 278.2 | 8.540e-34 | 10.01 | 1.216e-46
Table 8: Heteroscedasticity tests on tasks involving digits data set. Data set | Model | Metric | Fligner Statistic | p-value | Levene Statistic | p-value
---|---|---|---|---|---|---
Digits | DT | acc | 205 | 5.670e-21 | 14.29 | 9.219e-65
| MLP-adam | acc | 256.7 | 6.239e-30 | 7.342 | 1.219e-33
| MLP-SGD | acc | 210 | 8.188e-22 | 6.53 | 2.167e-29
| RF | acc | 184.3 | 1.458e-17 | 15.61 | 9.379e-70
| SVM | acc | 91.72 | 2.093e-04 | 2.187 | 1.790e-05
| ada | acc | 99.34 | 2.832e-05 | 2.305 | 4.601e-06
| kNN | acc | 35 | 0.9343 | 0.7042 | 0.9349
| lasso | acc | 22.97 | 0.9994 | 0.4292 | 0.9998
| linear | acc | 17.3 | 1 | 0.2963 | 1
| DT | nll | 249.6 | 1.140e-28 | 15.71 | 3.892e-70
| MLP-adam | nll | 339.8 | 3.816e-45 | 6.882 | 3.012e-31
| MLP-SGD | nll | 244.8 | 7.740e-28 | 6.104 | 4.129e-27
| RF | nll | 144 | 2.791e-11 | 7.435 | 4.059e-34
| SVM | nll | 4.373 | 1 | 0.06091 | 1
| ada | nll | 135.1 | 5.444e-10 | 3.294 | 2.061e-11
| kNN | nll | 108.2 | 2.326e-06 | 3.059 | 4.211e-10
| lasso | nll | 88.4 | 4.799e-04 | 2.116 | 3.995e-05
| linear | nll | 103 | 1.024e-05 | 3.328 | 1.335e-11
Table 9: Heteroscedasticity tests on tasks involving iris data set. Data set | Model | Metric | Fligner Statistic | p-value | Levene Statistic | p-value
---|---|---|---|---|---|---
Iris | DT | acc | 207.1 | 2.440e-21 | 6.523 | 2.355e-29
| MLP-adam | acc | 83.81 | 1.436e-03 | 1.838 | 7.989e-04
| MLP-SGD | acc | 68.52 | 0.03413 | 1.409 | 0.04082
| RF | acc | 155.5 | 5.311e-13 | 6.138 | 2.726e-27
| SVM | acc | 198.4 | 6.990e-20 | 3.345 | 1.065e-11
| ada | acc | 155.7 | 4.788e-13 | 5.018 | 3.858e-21
| kNN | acc | 55.68 | 0.2378 | 1.124 | 0.2701
| lasso | acc | 19.72 | 0.9999 | 0.4045 | 0.9999
| linear | acc | 106.4 | 3.965e-06 | 2.959 | 1.502e-09
| DT | nll | 322.2 | 7.375e-42 | 6.118 | 3.506e-27
| MLP-adam | nll | 106.3 | 4.070e-06 | 3.123 | 1.869e-10
| MLP-SGD | nll | 155.6 | 4.966e-13 | 6.386 | 1.264e-28
| RF | nll | 321.3 | 1.066e-41 | 8.339 | 1.136e-38
| SVM | nll | 188.4 | 3.217e-18 | 4.736 | 1.470e-19
| ada | nll | 74.04 | 0.01194 | 1.414 | 0.03938
| kNN | nll | 212.6 | 2.863e-22 | 8.838 | 4.118e-41
| lasso | nll | 45.45 | 0.6177 | 0.5045 | 0.998
| linear | nll | 36.64 | 0.9037 | 0.733 | 0.9101
Table 10: Heteroscedasticity tests on tasks involving the wine data set. Data set | Model | Metric | Fligner Statistic | p-value | Levene Statistic | p-value
---|---|---|---|---|---|---
Wine | DT | acc | 127.3 | 6.912e-09 | 3.553 | 7.195e-13
| MLP-adam | acc | 85.37 | 9.945e-04 | 1.874 | 5.544e-04
| MLP-SGD | acc | 109 | 1.845e-06 | 2.48 | 5.701e-07
| RF | acc | 128.5 | 4.717e-09 | 5.069 | 2.014e-21
| SVM | acc | 28.73 | 0.9908 | 0.5136 | 0.9975
| ada | acc | 156.6 | 3.527e-13 | 3.968 | 3.215e-15
| kNN | acc | 37.67 | 0.8807 | 0.6869 | 0.9473
| lasso | acc | 29.8 | 0.9862 | 0.5981 | 0.9859
| linear | acc | 21.28 | 0.9998 | 0.3839 | 1
| DT | nll | 349.2 | 6.614e-47 | 10.46 | 1.115e-48
| MLP-adam | nll | 57.19 | 0.1971 | 1.21 | 0.1646
| MLP-SGD | nll | 110.1 | 1.362e-06 | 2.597 | 1.380e-07
| RF | nll | 258 | 3.660e-30 | 6.468 | 4.597e-29
| SVM | nll | 57.18 | 0.1975 | 1.006 | 0.4663
| ada | nll | 152.8 | 1.323e-12 | 3.072 | 3.555e-10
| kNN | nll | 178.2 | 1.410e-16 | 5.446 | 1.635e-23
| lasso | nll | 83.94 | 1.394e-03 | 1.782 | 1.416e-03
| linear | nll | 185.8 | 8.404e-18 | 5.01 | 4.312e-21

##  D Task-Level Results Breakdown

In this section we present the full task-level breakdown of the results with
each metric, data set and model combination for each black-box optimiser
summarised with the mean and variance achieved across 20 seeds. We show a
summary plot in Table 11.

Table 11: Number of tasks for which each optimiser performed best.

HEBO | TuRBO | PySOT | Skopt | Nevergrad (1+1) | BOHB-BB | Opentuner | Hyperopt | TuRBO+
---|---|---|---|---|---|---|---|---
71 (65.7%) | 14 (13.0%) | 7 (6.5 %) | 5 (4.6 %) | 4 (3.7 %) | 3 (2.8%) | 2 (1.9%) | 1 (0.9%) | 1 (0.9%)

We now present sequentially the full results for each of the 6 datasets:
Boston (Section D.1), Breast cancer (Section D.2), Diabetes (Section D.3),
Digits (Section D.4), Iris (Section D.5), Wine (Section D.6). For each
optimiser we give the mean and variance of the performance metric across all
18 tasks (2 metrics x 9 models) for a given data set.

###  D.1 Boston Data Set

Table 12: Boston with MAE loss for tuning DT model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Boston | DT | MAE | 106.658094 | 3.514569
PySOT | Boston | DT | MAE | 104.811444 | 4.220520
Skopt | Boston | DT | MAE | 104.303237 | 4.311549
TuRBO | Boston | DT | MAE | 103.863640 | 6.184650
Nevergrad (1+1) | Boston | DT | MAE | 102.938925 | 17.667195
Hyperopt | Boston | DT | MAE | 99.231898 | 5.681478
TuRBO+ | Boston | DT | MAE | 98.464758 | 35.071496
Random-search | Boston | DT | MAE | 95.273667 | 21.448475
BOHB-BB | Boston | DT | MAE | 93.360381 | 37.490514
Opentuner | Boston | DT | MAE | 86.421660 | 313.239599
Table 13: Boston with MAE loss for tuning MLP-adam model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is TuRBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
TuRBO | Boston | MLP-adam | MAE | 104.423509 | 17.216156
PySOT | Boston | MLP-adam | MAE | 103.369026 | 15.037764
HEBO | Boston | MLP-adam | MAE | 101.377994 | 31.981786
TuRBO+ | Boston | MLP-adam | MAE | 98.759781 | 46.065189
Nevergrad (1+1) | Boston | MLP-adam | MAE | 96.279956 | 45.560329
Hyperopt | Boston | MLP-adam | MAE | 95.255250 | 28.728828
Random-search | Boston | MLP-adam | MAE | 92.866456 | 23.342156
BOHB-BB | Boston | MLP-adam | MAE | 92.001680 | 17.523725
Opentuner | Boston | MLP-adam | MAE | 89.595926 | 47.161525
Skopt | Boston | MLP-adam | MAE | 86.257482 | 24.245633
Table 14: Boston with MAE loss for tuning MLP-SGD model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is PySOT. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
PySOT | Boston | MLP-SGD | MAE | 96.325329 | 34.592455
HEBO | Boston | MLP-SGD | MAE | 95.885011 | 4.520840
TuRBO | Boston | MLP-SGD | MAE | 94.936543 | 21.232934
Hyperopt | Boston | MLP-SGD | MAE | 93.605852 | 4.184636
Skopt | Boston | MLP-SGD | MAE | 91.084570 | 8.855305
Random-search | Boston | MLP-SGD | MAE | 89.854633 | 9.131124
BOHB-BB | Boston | MLP-SGD | MAE | 89.622364 | 20.485508
TuRBO+ | Boston | MLP-SGD | MAE | 89.288302 | 18.617165
Nevergrad (1+1) | Boston | MLP-SGD | MAE | 87.149541 | 115.487770
Opentuner | Boston | MLP-SGD | MAE | 85.130342 | 26.218217
Table 15: Boston with MAE loss for tuning RF model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Boston | RF | MAE | 101.726637 | 0.282878
TuRBO | Boston | RF | MAE | 100.095822 | 2.339170
PySOT | Boston | RF | MAE | 99.773979 | 1.147489
Skopt | Boston | RF | MAE | 99.352605 | 0.805112
Hyperopt | Boston | RF | MAE | 98.142703 | 1.649156
Nevergrad (1+1) | Boston | RF | MAE | 96.790676 | 20.349353
TuRBO+ | Boston | RF | MAE | 95.817383 | 7.005061
Opentuner | Boston | RF | MAE | 94.012783 | 32.962245
BOHB-BB | Boston | RF | MAE | 93.869871 | 6.026597
Random-search | Boston | RF | MAE | 93.223774 | 8.626423
Table 16: Boston with MAE loss for tuning SVM model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Boston | SVM | MAE | 102.469764 | 0.015274
Nevergrad (1+1) | Boston | SVM | MAE | 102.236442 | 0.135306
TuRBO | Boston | SVM | MAE | 101.826793 | 0.169932
PySOT | Boston | SVM | MAE | 101.262466 | 0.082329
Skopt | Boston | SVM | MAE | 100.914194 | 0.281085
Hyperopt | Boston | SVM | MAE | 98.943800 | 9.092074
Opentuner | Boston | SVM | MAE | 97.906378 | 12.428820
TuRBO+ | Boston | SVM | MAE | 90.935974 | 182.645481
Random-search | Boston | SVM | MAE | 88.908494 | 80.396927
BOHB-BB | Boston | SVM | MAE | 86.714806 | 46.398841
Table 17: Boston with MAE loss for tuning Ada model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Boston | Ada | MAE | 103.445454 | 19.376575
Opentuner | Boston | Ada | MAE | 101.293346 | 54.800061
TuRBO | Boston | Ada | MAE | 99.774003 | 84.126539
Nevergrad (1+1) | Boston | Ada | MAE | 97.261376 | 394.018233
Skopt | Boston | Ada | MAE | 97.141270 | 48.739222
TuRBO+ | Boston | Ada | MAE | 97.035052 | 178.220595
Hyperopt | Boston | Ada | MAE | 93.394811 | 59.032286
PySOT | Boston | Ada | MAE | 91.618187 | 25.548574
Random-search | Boston | Ada | MAE | 85.889844 | 39.688171
BOHB-BB | Boston | Ada | MAE | 84.104672 | 99.526333
Table 18: Boston with MAE loss for tuning Knn model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Boston | Knn | MAE | 100.000000 | 0.000000
Opentuner | Boston | Knn | MAE | 100.000000 | 0.000000
Skopt | Boston | Knn | MAE | 100.000000 | 0.000000
TuRBO | Boston | Knn | MAE | 100.000000 | 0.000000
TuRBO+ | Boston | Knn | MAE | 99.689888 | 1.923386
PySOT | Boston | Knn | MAE | 99.069665 | 5.162772
Hyperopt | Boston | Knn | MAE | 98.023553 | 9.829689
Nevergrad (1+1) | Boston | Knn | MAE | 97.379434 | 28.707287
BOHB-BB | Boston | Knn | MAE | 97.134442 | 18.411664
Random-search | Boston | Knn | MAE | 96.297696 | 45.566838
Table 19: Boston with MAE loss for tuning Lasso model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Boston | Lasso | MAE | 100.031658 | 0.000141
PySOT | Boston | Lasso | MAE | 100.004491 | 0.001119
Nevergrad (1+1) | Boston | Lasso | MAE | 99.989103 | 0.005199
Skopt | Boston | Lasso | MAE | 99.979143 | 0.001668
TuRBO | Boston | Lasso | MAE | 99.977645 | 0.002141
Hyperopt | Boston | Lasso | MAE | 99.967041 | 0.002263
BOHB-BB | Boston | Lasso | MAE | 99.926445 | 0.006747
TuRBO+ | Boston | Lasso | MAE | 99.921337 | 0.031836
Random-search | Boston | Lasso | MAE | 99.917670 | 0.002902
Opentuner | Boston | Lasso | MAE | 99.306403 | 0.956136
Table 20: Boston with MAE loss for tuning Linear model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Boston | Linear | MAE | 99.964657 | 0.002542
Hyperopt | Boston | Linear | MAE | 99.922142 | 0.006819
TuRBO | Boston | Linear | MAE | 99.695751 | 1.213154
PySOT | Boston | Linear | MAE | 99.394823 | 0.799652
BOHB-BB | Boston | Linear | MAE | 98.547907 | 5.197760
Skopt | Boston | Linear | MAE | 98.327627 | 6.150494
Random-search | Boston | Linear | MAE | 97.843101 | 5.252852
TuRBO+ | Boston | Linear | MAE | 95.874716 | 153.782738
Nevergrad (1+1) | Boston | Linear | MAE | 80.996813 | 1523.062377
Opentuner | Boston | Linear | MAE | 45.221227 | 2080.382131
Table 21: Boston with MSE loss for tuning DT model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is Skopt. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
Skopt | Boston | DT | MSE | 105.587343 | 12.872983
PySOT | Boston | DT | MSE | 105.072075 | 4.900371
HEBO | Boston | DT | MSE | 104.838117 | 6.188087
TuRBO | Boston | DT | MSE | 104.030707 | 10.799242
Hyperopt | Boston | DT | MSE | 102.229770 | 11.149492
TuRBO+ | Boston | DT | MSE | 100.450728 | 12.162198
Nevergrad (1+1) | Boston | DT | MSE | 97.629442 | 407.356319
Random-search | Boston | DT | MSE | 95.359894 | 28.675740
BOHB-BB | Boston | DT | MSE | 94.520125 | 20.675992
Opentuner | Boston | DT | MSE | 87.826986 | 205.761429
Table 22: Boston with MSE loss for tuning MLP-adam model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Boston | MLP-adam | MSE | 101.227304 | 15.854379
TuRBO | Boston | MLP-adam | MSE | 100.906498 | 4.709893
PySOT | Boston | MLP-adam | MSE | 98.878397 | 7.952874
TuRBO+ | Boston | MLP-adam | MSE | 98.090398 | 17.299141
Nevergrad (1+1) | Boston | MLP-adam | MSE | 97.631386 | 45.709714
Hyperopt | Boston | MLP-adam | MSE | 95.833644 | 15.829315
BOHB-BB | Boston | MLP-adam | MSE | 95.267050 | 21.466780
Random-search | Boston | MLP-adam | MSE | 93.691552 | 18.033615
Opentuner | Boston | MLP-adam | MSE | 92.489905 | 52.692896
Skopt | Boston | MLP-adam | MSE | 87.610428 | 15.511744
Table 23: Boston with MSE loss for tuning MLP-SGD model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is TuRBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
TuRBO | Boston | MLP-SGD | MSE | 104.817056 | 17.856834
HEBO | Boston | MLP-SGD | MSE | 103.285788 | 4.467401
Hyperopt | Boston | MLP-SGD | MSE | 102.049189 | 6.848005
PySOT | Boston | MLP-SGD | MSE | 99.992157 | 17.117778
TuRBO+ | Boston | MLP-SGD | MSE | 98.973335 | 45.162334
Skopt | Boston | MLP-SGD | MSE | 98.293397 | 6.746134
BOHB-BB | Boston | MLP-SGD | MSE | 94.731025 | 17.820451
Random-search | Boston | MLP-SGD | MSE | 94.524667 | 28.931228
Nevergrad (1+1) | Boston | MLP-SGD | MSE | 94.267722 | 207.308935
Opentuner | Boston | MLP-SGD | MSE | 88.380591 | 58.962581
Table 24: Boston with MSE loss for tuning RF model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Boston | RF | MSE | 103.461058 | 0.440697
PySOT | Boston | RF | MSE | 101.853905 | 0.492860
TuRBO | Boston | RF | MSE | 101.839078 | 0.800142
Skopt | Boston | RF | MSE | 101.472976 | 0.606579
Nevergrad (1+1) | Boston | RF | MSE | 100.814244 | 2.514980
Hyperopt | Boston | RF | MSE | 100.547643 | 0.964402
TuRBO+ | Boston | RF | MSE | 98.915015 | 3.362617
Random-search | Boston | RF | MSE | 98.167459 | 4.278077
Opentuner | Boston | RF | MSE | 98.049981 | 11.073665
BOHB-BB | Boston | RF | MSE | 97.839608 | 6.645583
Table 25: Boston with MSE loss for tuning SVM model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Boston | SVM | MSE | 104.852372 | 0.005667
Nevergrad (1+1) | Boston | SVM | MSE | 104.578053 | 0.288069
TuRBO | Boston | SVM | MSE | 103.821417 | 1.334006
PySOT | Boston | SVM | MSE | 103.398778 | 0.101558
Hyperopt | Boston | SVM | MSE | 101.271084 | 3.903429
Skopt | Boston | SVM | MSE | 100.753282 | 76.344135
TuRBO+ | Boston | SVM | MSE | 97.946333 | 50.658268
Opentuner | Boston | SVM | MSE | 91.009871 | 467.963686
BOHB-BB | Boston | SVM | MSE | 90.350500 | 127.546283
Random-search | Boston | SVM | MSE | 83.204838 | 135.308788
Table 26: Boston with MSE loss for tuning Ada model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Boston | Ada | MSE | 96.942494 | 18.658049
Opentuner | Boston | Ada | MSE | 94.274205 | 34.779775
TuRBO | Boston | Ada | MSE | 94.131115 | 43.476954
Nevergrad (1+1) | Boston | Ada | MSE | 92.333218 | 117.061847
TuRBO+ | Boston | Ada | MSE | 91.017786 | 86.796910
Hyperopt | Boston | Ada | MSE | 89.721591 | 41.637778
PySOT | Boston | Ada | MSE | 88.996278 | 29.924305
Skopt | Boston | Ada | MSE | 87.835681 | 55.855654
BOHB-BB | Boston | Ada | MSE | 82.947338 | 67.514235
Random-search | Boston | Ada | MSE | 80.119051 | 17.165168
Table 27: Boston with MSE loss for tuning Knn model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Boston | Knn | MSE | 100.000000 | 0.000000
Opentuner | Boston | Knn | MSE | 100.000000 | 0.000000
PySOT | Boston | Knn | MSE | 100.000000 | 0.000000
TuRBO | Boston | Knn | MSE | 100.000000 | 0.000000
TuRBO+ | Boston | Knn | MSE | 99.704490 | 1.253051
Hyperopt | Boston | Knn | MSE | 99.319467 | 6.373596
Random-search | Boston | Knn | MSE | 99.068697 | 7.272026
Skopt | Boston | Knn | MSE | 99.014267 | 19.433403
BOHB-BB | Boston | Knn | MSE | 96.129465 | 62.658638
Nevergrad (1+1) | Boston | Knn | MSE | 88.882774 | 227.331571
Table 28: Boston with MSE loss for tuning Lasso model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Boston | Lasso | MSE | 100.037610 | 0.000010
PySOT | Boston | Lasso | MSE | 100.008571 | 0.000269
Hyperopt | Boston | Lasso | MSE | 100.007350 | 0.000123
TuRBO | Boston | Lasso | MSE | 99.982366 | 0.021257
Skopt | Boston | Lasso | MSE | 99.806671 | 0.433679
BOHB-BB | Boston | Lasso | MSE | 99.679375 | 0.206443
Nevergrad (1+1) | Boston | Lasso | MSE | 99.372815 | 4.020588
Random-search | Boston | Lasso | MSE | 99.325799 | 0.769118
TuRBO+ | Boston | Lasso | MSE | 98.971017 | 2.235443
Opentuner | Boston | Lasso | MSE | 97.388970 | 6.008643
Table 29: Boston with MSE loss for tuning Linear model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Boston | Linear | MSE | 100.000001 | 3.207627e-10
Hyperopt | Boston | Linear | MSE | 99.999348 | 1.340314e-06
PySOT | Boston | Linear | MSE | 99.994196 | 1.378338e-04
Skopt | Boston | Linear | MSE | 99.986643 | 6.872575e-04
Nevergrad (1+1) | Boston | Linear | MSE | 99.976999 | 5.017296e-03
TuRBO | Boston | Linear | MSE | 99.889263 | 2.451020e-01
BOHB-BB | Boston | Linear | MSE | 99.837413 | 2.455636e-01
Random-search | Boston | Linear | MSE | 99.802681 | 2.731228e-01
TuRBO+ | Boston | Linear | MSE | 99.332841 | 1.085780e+00
Opentuner | Boston | Linear | MSE | 98.545654 | 1.344348e+00

###  D.2 Breast Cancer Data set

Table 30: Breast cancer dataset with ACC loss for tuning DT model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Breast cancer dataset | DT | ACC | 106.694915 | 52.457703
PySOT | Breast cancer dataset | DT | ACC | 103.220339 | 32.930646
TuRBO | Breast cancer dataset | DT | ACC | 99.830508 | 32.930646
Nevergrad (1+1) | Breast cancer dataset | DT | ACC | 99.661017 | 47.052420
Skopt | Breast cancer dataset | DT | ACC | 99.322034 | 27.033974
Opentuner | Breast cancer dataset | DT | ACC | 97.457627 | 19.201984
Hyperopt | Breast cancer dataset | DT | ACC | 97.288136 | 17.357384
TuRBO+ | Breast cancer dataset | DT | ACC | 96.610169 | 26.913017
Random-search | Breast cancer dataset | DT | ACC | 93.813559 | 29.415322
BOHB-BB | Breast cancer dataset | DT | ACC | 93.305085 | 22.823145
Table 31: Breast cancer dataset with ACC loss for tuning MLP-adam model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is Hyperopt. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
Hyperopt | Breast cancer dataset | MLP-adam | ACC | 99.078947 | 14.779851
HEBO | Breast cancer dataset | MLP-adam | ACC | 98.684211 | 20.046654
Nevergrad (1+1) | Breast cancer dataset | MLP-adam | ACC | 96.973684 | 53.415221
PySOT | Breast cancer dataset | MLP-adam | ACC | 96.578947 | 11.736405
TuRBO | Breast cancer dataset | MLP-adam | ACC | 96.447368 | 17.695728
TuRBO+ | Breast cancer dataset | MLP-adam | ACC | 96.447368 | 18.424697
Skopt | Breast cancer dataset | MLP-adam | ACC | 95.921053 | 20.757399
Opentuner | Breast cancer dataset | MLP-adam | ACC | 95.000000 | 21.067211
Random-search | Breast cancer dataset | MLP-adam | ACC | 94.078947 | 13.759294
BOHB-BB | Breast cancer dataset | MLP-adam | ACC | 93.684211 | 11.955095
Table 32: Breast cancer dataset with ACC loss for tuning MLP-SGD model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Breast cancer dataset | MLP-SGD | ACC | 100.847458 | 1.209574
TuRBO | Breast cancer dataset | MLP-SGD | ACC | 99.491525 | 1.391010
Skopt | Breast cancer dataset | MLP-SGD | ACC | 99.237288 | 0.824022
PySOT | Breast cancer dataset | MLP-SGD | ACC | 99.067797 | 1.504407
Hyperopt | Breast cancer dataset | MLP-SGD | ACC | 98.771186 | 1.396680
Nevergrad (1+1) | Breast cancer dataset | MLP-SGD | ACC | 98.771186 | 4.647409
TuRBO+ | Breast cancer dataset | MLP-SGD | ACC | 98.559322 | 2.502306
Random-search | Breast cancer dataset | MLP-SGD | ACC | 96.864407 | 1.973117
BOHB-BB | Breast cancer dataset | MLP-SGD | ACC | 96.355932 | 3.787478
Opentuner | Breast cancer dataset | MLP-SGD | ACC | 84.703390 | 190.241386
Table 33: Breast cancer dataset with ACC loss for tuning RF model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is Skopt. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
Skopt | Breast cancer dataset | RF | ACC | 98.846154 | 5.985365
HEBO | Breast cancer dataset | RF | ACC | 98.356643 | 8.250307
PySOT | Breast cancer dataset | RF | ACC | 97.377622 | 5.964775
Hyperopt | Breast cancer dataset | RF | ACC | 96.888112 | 3.988099
TuRBO | Breast cancer dataset | RF | ACC | 96.853147 | 8.416317
Nevergrad (1+1) | Breast cancer dataset | RF | ACC | 96.083916 | 13.249908
TuRBO+ | Breast cancer dataset | RF | ACC | 95.419580 | 6.047137
BOHB-BB | Breast cancer dataset | RF | ACC | 94.335664 | 2.980457
Opentuner | Breast cancer dataset | RF | ACC | 93.636364 | 3.289313
Random-search | Breast cancer dataset | RF | ACC | 93.321678 | 2.495296
Table 34: Breast cancer dataset with ACC loss for tuning SVM model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is TuRBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
TuRBO | Breast cancer dataset | SVM | ACC | 89.285714 | 40.279270
Hyperopt | Breast cancer dataset | SVM | ACC | 86.428571 | 53.168636
TuRBO+ | Breast cancer dataset | SVM | ACC | 86.428571 | 96.133190
PySOT | Breast cancer dataset | SVM | ACC | 85.714286 | 42.964554
HEBO | Breast cancer dataset | SVM | ACC | 84.285714 | 40.816327
BOHB-BB | Breast cancer dataset | SVM | ACC | 80.714286 | 113.319012
Opentuner | Breast cancer dataset | SVM | ACC | 78.571429 | 182.599356
Nevergrad (1+1) | Breast cancer dataset | SVM | ACC | 77.142857 | 481.203007
Random-search | Breast cancer dataset | SVM | ACC | 76.428571 | 48.872180
Skopt | Breast cancer dataset | SVM | ACC | 76.428571 | 349.624060
Table 35: Breast cancer dataset with ACC loss for tuning Ada model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Breast cancer dataset | Ada | ACC | 99.06250 | 4.214638
Opentuner | Breast cancer dataset | Ada | ACC | 98.59375 | 6.656044
TuRBO | Breast cancer dataset | Ada | ACC | 98.43750 | 3.597862
TuRBO+ | Breast cancer dataset | Ada | ACC | 97.96875 | 8.506373
Nevergrad (1+1) | Breast cancer dataset | Ada | ACC | 97.65625 | 12.207031
PySOT | Breast cancer dataset | Ada | ACC | 97.03125 | 5.628084
Hyperopt | Breast cancer dataset | Ada | ACC | 96.87500 | 8.223684
Random-search | Breast cancer dataset | Ada | ACC | 95.31250 | 6.681743
BOHB-BB | Breast cancer dataset | Ada | ACC | 94.53125 | 3.983347
Skopt | Breast cancer dataset | Ada | ACC | 93.12500 | 9.868421
Table 36: Breast cancer dataset with ACC loss for tuning Knn model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Breast cancer dataset | Knn | ACC | 100.000 | 0.000000
Opentuner | Breast cancer dataset | Knn | ACC | 100.000 | 0.000000
Hyperopt | Breast cancer dataset | Knn | ACC | 97.500 | 26.315789
Nevergrad (1+1) | Breast cancer dataset | Knn | ACC | 95.625 | 53.865132
Random-search | Breast cancer dataset | Knn | ACC | 95.625 | 53.865132
Skopt | Breast cancer dataset | Knn | ACC | 95.000 | 55.921053
BOHB-BB | Breast cancer dataset | Knn | ACC | 95.000 | 39.473684
TuRBO+ | Breast cancer dataset | Knn | ACC | 92.500 | 55.921053
TuRBO | Breast cancer dataset | Knn | ACC | 91.875 | 37.417763
PySOT | Breast cancer dataset | Knn | ACC | 91.250 | 34.539474
Table 37: Breast cancer dataset with ACC loss for tuning Lasso model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is PySOT. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
PySOT | Breast cancer dataset | Lasso | ACC | 95.000000 | 32.163743
TuRBO | Breast cancer dataset | Lasso | ACC | 93.888889 | 32.163743
Opentuner | Breast cancer dataset | Lasso | ACC | 93.888889 | 32.163743
HEBO | Breast cancer dataset | Lasso | ACC | 93.333333 | 31.189084
TuRBO+ | Breast cancer dataset | Lasso | ACC | 92.222222 | 27.290448
Hyperopt | Breast cancer dataset | Lasso | ACC | 90.555556 | 16.569201
Random-search | Breast cancer dataset | Lasso | ACC | 90.555556 | 16.569201
BOHB-BB | Breast cancer dataset | Lasso | ACC | 90.000000 | 11.695906
Nevergrad (1+1) | Breast cancer dataset | Lasso | ACC | 88.333333 | 58.154646
Skopt | Breast cancer dataset | Lasso | ACC | 87.777778 | 180.636777
Table 38: Breast cancer dataset with ACC loss for tuning Linear model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Breast cancer dataset | Linear | ACC | 107.0 | 95.789474
Nevergrad (1+1) | Breast cancer dataset | Linear | ACC | 103.0 | 390.526316
TuRBO | Breast cancer dataset | Linear | ACC | 99.0 | 188.421053
PySOT | Breast cancer dataset | Linear | ACC | 98.0 | 164.210526
TuRBO+ | Breast cancer dataset | Linear | ACC | 98.0 | 164.210526
Opentuner | Breast cancer dataset | Linear | ACC | 98.0 | 290.526316
Hyperopt | Breast cancer dataset | Linear | ACC | 96.0 | 109.473684
Skopt | Breast cancer dataset | Linear | ACC | 95.0 | 205.263158
Random-search | Breast cancer dataset | Linear | ACC | 89.0 | 146.315789
BOHB-BB | Breast cancer dataset | Linear | ACC | 85.0 | 78.947368
Table 39: Breast cancer dataset with NLL loss for tuning DT model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Breast cancer dataset | DT | NLL | 111.047413 | 51.627976
TuRBO+ | Breast cancer dataset | DT | NLL | 98.912425 | 121.134669
PySOT | Breast cancer dataset | DT | NLL | 98.696805 | 115.036349
Hyperopt | Breast cancer dataset | DT | NLL | 98.201008 | 75.868808
Skopt | Breast cancer dataset | DT | NLL | 97.647135 | 165.617171
TuRBO | Breast cancer dataset | DT | NLL | 95.406797 | 230.498379
Opentuner | Breast cancer dataset | DT | NLL | 91.076386 | 92.789945
Nevergrad (1+1) | Breast cancer dataset | DT | NLL | 90.579282 | 92.638289
BOHB-BB | Breast cancer dataset | DT | NLL | 90.418909 | 43.016544
Random-search | Breast cancer dataset | DT | NLL | 87.752435 | 64.953742
Table 40: Breast cancer dataset with NLL loss for tuning MLP-adam model, averaged over 20 seeds. Optimiser for this task with highest mean normalised score is HEBO. Algorithm | Dataset | Model | Metric | Normalised Score | Variance
---|---|---|---|---|---
HEBO | Breast cancer dataset | MLP-adam | NLL | 100.303210 | 0.252546
TuRBO | Breast cancer dataset | MLP-adam | NLL | 100.260440 | 0.386150
PySOT | Breast cancer dataset | MLP-adam | NLL | 99.965212 | 0.410841
Hyperopt | Breast cancer dataset | MLP-adam | NLL | 99.778153 | 0.295595
Opentuner | Breast cancer dataset | MLP-adam | NLL | 99.614006 | 0.559840
Skopt | Breast cancer dataset | MLP-adam | NLL | 99.422393 | 0.444736
TuRBO+ | Breast cancer dataset | MLP-adam | NLL | 99.298544 | 1.362840
Nevergrad (1+1) | Breast cancer dataset | MLP-adam | NLL | 98.865542 | 5.064746
Random-search | Breast cancer dataset | MLP-adam | NLL | 98.538089 | 1.227034
BOHB-BB | Breast cancer dataset | MLP-adam | NLL | 98.396155 | 0.577588