Modulating Surrogates for Bayesian Optimization

TL;DR

This paper introduces surrogate models that emphasize well-behaved parts of the objective function in Bayesian optimization, improving robustness and efficiency on complex real-world problems.

Contribution

It proposes a novel approach using latent Gaussian processes to focus on informative structures, ignoring challenging details that hinder optimization progress.

Findings

Latent Gaussian process surrogates outperform standard models.

The approach enhances robustness on difficult objective functions.

Experimental results show improved optimization reliability.

Abstract

Bayesian optimization (BO) methods often rely on the assumption that the objective function is well-behaved, but in practice, this is seldom true for real-world objectives even if noise-free observations can be collected. Common approaches, which try to model the objective as precisely as possible, often fail to make progress by spending too many evaluations modeling irrelevant details. We address this issue by proposing surrogate models that focus on the well-behaved structure in the objective function, which is informative for search, while ignoring detrimental structure that is challenging to model from few observations. First, we demonstrate that surrogate models with appropriate noise distributions can absorb challenging structures in the objective function by treating them as irreducible uncertainty. Secondly, we show that a latent Gaussian process is an excellent surrogate for this purpose, comparing with Gaussian processes with standard noise distributions. We perform numerous experiments on a range of BO benchmarks and find that our approach improves reliability and performance when faced with challenging objective functions.

Tables

Table 1. Table 1: Mean gap performance for various test functions; higher is better. The upper table shows the results after 50 objective function evaluations and the lower table after 100 evaluations. Due to computational cost, Warped GP results are only reported for 50 evaluations. Methods not significantly different from the best performing method with respect by a two-sided paired Wilcoxon signed-rank test at a 5% significance level over 20 repetitions are shown in bold  (Malkomes & Garnett, 2018 ) . For results in terms of regret, see the supplement.

Benchmark	Evals	Dim	Properties	GP	Warped GP	Homosced GP	Heterosced GP	LGP
Hartmann	50	6	boring	0.959	0.537	0.881	0.973	0.937
Griewank	50	2	oscillatory	0.914	0.493	0.752	0.913	0.897
Shubert	50	2	oscillatory	0.378	0.158	0.378	0.480	0.593
Ackley ${[-10,30]}^d$	50	2	complicated, oscillatory	0.924	0.274	0.892	0.912	0.927
Cross In Tray	50	2	complicated, oscillatory	0.954	0.385	0.929	0.977	0.945
Holder table	50	2	complicated, oscillatory	0.939	0.896	0.900	0.931	0.993
Corrupted Holder Table	50	2	complicated, oscillatory	0.741	0.798	0.826	0.729	0.896
Branin01	100	2	none	1.000		1.000	1.000	1.000
Branin02	100	2	none	0.991		0.964	0.990	0.981
Beale	100	2	boring	0.987		0.982	0.987	0.988
Hartmann	100	6	boring	0.987		0.947	0.984	0.979
Griewank	100	2	oscillatory	0.967		0.875	0.969	0.946
Levy	100	2	oscillatory	0.997		0.999	0.998	0.998
Deflected Corrugated Spring	100	10	oscillatory	0.347		0.840	0.406	0.697
Shubert ${[-10,10]}^d$	100	2	oscillatory	0.510		0.511	0.672	0.877
Weierstrass	100	8	complicated	0.600		0.704	0.577	0.625
Cross In Tray	100	2	complicated, oscillatory	1.000		0.995	1.000	1.000
Holder Table	100	2	complicated, oscillatory	0.971		0.963	0.964	1.000
Ackley ${[-10,30]}^d$	100	2	complicated, oscillatory	0.971		0.914	0.980	0.974
Ackley ${[-10,30]}^d$	100	6	complicated, oscillatory	0.459		0.789	0.442	0.712
Corrupted Holder Table	100	2	complicated, oscillatory	0.844		0.889	0.822	0.918
Corrupted Exponential	100	8	complicated, oscillatory	0.580		0.847	0.581	0.806
HPO: NN Boston	100	9	unknown	0.720		0.761	0.810	0.770
HPO: NN Climate Model Crashes	100	9	unknown	0.629		0.717	0.683	0.678
Active learning: Robot Pushing	100	4	unknown	0.877		0.745	0.907	0.932

Table 2. Table 2: Function properties as defined in SigOpt  (McCourt, 2016 ) .

Property	Meaning
boring	A mostly boring function that only has a small region of action
oscillatory	A function with a general trend and an short range oscillatory component
complicated	These are functions that may fit a behavior, but not in the most obvious or satisfying way

Table 3. Table 3: Benchmark function domains used as specified in SigOpt  (McCourt, 2016 ) ’evalset’. The domains of Corrupted Holder Table and Corrupted Exponential correspond to the ones of Holder Table and Exponential, respectively.

Benchmark	Dim	Properties	Domain
Branin01	2	none	$[[-5,10],[0,15]]$
Branin02	2	none	$[[-5,15],[-5,15]]$
Powell Triple Log	12	none	${[-4,1]}^{12}$
Beale	2	boring	${[-4.5,4.5]}^2$
Hartmann	6	boring	${[0,1]}^6$
Griewank	2	oscillatory	${[-50,20]}^2$
Shubert01	2	oscillatory	${[-10,10]}^2$
Levy13	2	oscillatory	${[-10,10]}^2$
Cosine Mixture	10	oscillatory	${[-1,1]}^{10}$
Drop-Wave	10	oscillatory	${[-2,5.12]}^{10}$
Deflected Corrugated Spring	10	oscillatory	${[0,7.5]}^{10}$
Weierstrass	8	complicated	${[-0.5,0.2]}^8$
Cross In Tray	2	complicated, oscillatory	${[-10,10]}^2$
Holder Table	2	complicated, oscillatory	${[-10,10]}^2$
Ackley ${[-10,30]}^d$	2	complicated, oscillatory	${[-10,30]}^2$
Ackley ${[-10,30]}^d$	6	complicated, oscillatory	${[-10,30]}^6$
Corrupted Holder Table	2	complicated, oscillatory	${[-10,10]}^2$
Corrupted Exponential	8	complicated, oscillatory	${[-0.7,0.2]}^8$
NN Boston	9	unknown	Table 4
NN Climate Model Crashes	9	unknown	Table 4
Robot Pushing	4	unknown	Table 5

Table 4. Table 4: Neural Network hyperparameter domains as from  (Head et al., 2018 ) . Categorical options are set as specified below.

Parameter	Domain
hidden_layer_sizes	pow2int([1.0, 10.0])
alpha	pow10([-5.0, -1])
batch_size	pow2int([5.0, 10.0])
max_iter	pow2int([5.0, 8.0])
learning_rate_init	pow10([-5.0, -1])
power_t	[0.01, 0.99]
momentum	[0.1, 0.98]
beta_1	[0.1, 0.98]
beta_2	[0.1, 0.9999999]
learning_rate	constant
solver	adam
activation	relu
nesterovs_momentum	False

Table 5. Table 5: Active learning task of ’robot pushing’ as from  (Kaelbling & Lozano-Pérez, 2017 ) . Code is available at https://github.com/zi-w/Max-value-Entropy-Search . All instances of np.random.normal(0, 0.01) was replaced by np.random.normal(0, 1e-6) in push_world.py to make the function virtually noise-free. The goal location was set to a fixed location for reproducibility, as specified below.

Parameter	Domain
${\text{robot}}_x$	[-5, 5]
${\text{robot}}_y$	[-5, 5]
${\theta }_{\text{initial}}$	[0, 2$\pi$]
simulation steps	int([10., 300.])
${\text{goal}}_x$	3.0
${\text{goal}}_y$	2.0

Table 6. Table 6: Regret for various test functions; lower is better. The upper table shows the results after 50 objective function evaluations and the lower table after 100 evaluations. Due to computational cost, Warped GP results are only reported after 50 evaluations.

Benchmark	Evals	Dim	Func Max	Func Min	Initial regret	GP	Warped GP	Homosced GP	Heterosced GP	LGP
Hartmann	50	6	0.00	-3.32	2.764 $\pm$ 0.490	0.117 $\pm$ 0.109	1.360 $\pm$ 0.783	0.343 $\pm$ 0.575	0.074 $\pm$ 0.081	0.190 $\pm$ 0.436
Griewank	50	2	3.19	0.00	1.151 $\pm$ 0.441	0.104 $\pm$ 0.230	0.379 $\pm$ 0.179	0.220 $\pm$ 0.150	0.105 $\pm$ 0.101	0.102 $\pm$ 0.081
Shubert	50	2	210.45	-186.73	183.2 $\pm$ 9.820	114.6 $\pm$ 51.93	129.3 $\pm$ 31.52	113.1 $\pm$ 47.77	95.70 $\pm$ 60.38	75.34 $\pm$ 66.56
Ackley ${[-10,30]}^d$	50	2	22.27	0.00	16.45 $\pm$ 4.202	1.391 $\pm$ 2.101	12.79 $\pm$ 3.550	1.641 $\pm$ 0.723	1.282 $\pm$ 2.481	1.093 $\pm$ 0.741
Cross In Tray	50	2	-0.26	-2.06	0.393 $\pm$ 0.158	0.018 $\pm$ 0.052	0.195 $\pm$ 0.123	0.026 $\pm$ 0.032	0.009 $\pm$ 0.038	0.026 $\pm$ 0.062
Holder Table	50	2	0.00	-19.21	15.51 $\pm$ 2.571	0.983 $\pm$ 1.329	1.433 $\pm$ 0.687	1.451 $\pm$ 1.182	1.081 $\pm$ 1.376	0.112 $\pm$ 0.398
Corrupted Holder Table	50	2	3.46	-20.99	18.37 $\pm$ 2.444	4.816 $\pm$ 4.671	3.508 $\pm$ 0.983	3.138 $\pm$ 1.123	5.157 $\pm$ 5.433	1.836 $\pm$ 0.690
Branin01	100	2	308.13	0.40	6.975 $\pm$ 5.110	0.001 $\pm$ 0.000		0.001 $\pm$ 0.001	0.001 $\pm$ 0.001	0.000 $\pm$ 0.000
Branin02	100	2	506.98	5.56	27.666 $\pm$ 28.450	0.253 $\pm$ 0.500		0.736 $\pm$ 0.690	0.253 $\pm$ 0.501	0.380 $\pm$ 0.572
Beale	100	2	181853.61	0.00	1152.964 $\pm$ 3639.596	0.288 $\pm$ 0.302		0.259 $\pm$ 0.210	0.241 $\pm$ 0.253	0.250 $\pm$ 0.251
Hartmann	100	6	0.00	-3.32	2.764 $\pm$ 0.490	0.038 $\pm$ 0.082		0.162 $\pm$ 0.471	0.044 $\pm$ 0.066	0.058 $\pm$ 0.100
Griewank	100	2	3.19	0.00	1.151 $\pm$ 0.441	0.028 $\pm$ 0.020		0.099 $\pm$ 0.054	0.032 $\pm$ 0.023	0.055 $\pm$ 0.056
Levy	100	2	454.13	0.00	58.26 $\pm$ 37.30	0.289 $\pm$ 1.018		0.026 $\pm$ 0.028	0.205 $\pm$ 0.727	0.034 $\pm$ 0.025
Shubert ${[-10,10]}^d$	100	2	210.45	-186.73	183.2 $\pm$ 9.820	90.09 $\pm$ 56.44		88.72 $\pm$ 50.59	60.34 $\pm$ 58.95	22.70 $\pm$ 45.82
Deflected Corrugated Spring	100	10	25.87	-1.00	5.919 $\pm$ 1.364	3.921 $\pm$ 1.848		0.944 $\pm$ 0.574	3.411 $\pm$ 1.243	1.797 $\pm$ 0.880
Weierstrass	100	8	144.00	112.00	12.90 $\pm$ 1.430	5.070 $\pm$ 0.997		3.790 $\pm$ 0.767	5.401 $\pm$ 1.320	4.798 $\pm$ 1.436
Cross In Tray	100	2	-0.26	-2.06	0.393 $\pm$ 0.158	0.000 $\pm$ 0.000		0.001 $\pm$ 0.002	0.000 $\pm$ 0.000	0.000 $\pm$ 0.000
Holder Table	100	2	0.00	-19.21	15.51 $\pm$ 2.571	0.459 $\pm$ 1.048		0.495 $\pm$ 1.079	0.591 $\pm$ 1.176	0.005 $\pm$ 0.006
Ackley ${[-10,30]}^d$	100	2	22.27	0.00	16.45 $\pm$ 4.202	0.480 $\pm$ 0.720		1.318 $\pm$ 0.556	0.347 $\pm$ 0.597	0.395 $\pm$ 0.332
Ackley ${[-10,30]}^d$	100	6	22.27	0.00	19.89 $\pm$ 0.809	10.92 $\pm$ 6.656		4.191 $\pm$ 0.751	11.20 $\pm$ 6.438	5.819 $\pm$ 5.449
Corrupted Holder Table	100	2	3.46	-20.99	18.37 $\pm$ 2.444	2.921 $\pm$ 3.728		2.034 $\pm$ 0.781	3.396 $\pm$ 3.798	1.474 $\pm$ 0.342
Corrupted Exponential	100	8	-0.04	-0.99	0.391 $\pm$ 0.101	0.172 $\pm$ 0.109		0.057 $\pm$ 0.016	0.165 $\pm$ 0.087	0.074 $\pm$ 0.038
HPO: NN Boston	100	9	5.00	-0.85	0.886 $\pm$ 1.199	0.023 $\pm$ 0.012		0.023 $\pm$ 0.009	0.019 $\pm$ 0.013	0.021 $\pm$ 0.010
HPO: NN Climate Model Crashes	100	9	5.00	0.11	0.161 $\pm$ 0.099	0.047 $\pm$ 0.017		0.035 $\pm$ 0.017	0.040 $\pm$ 0.021	0.039 $\pm$ 0.017
Active learning: Robot Pushing	100	4	unknown	0.00	3.783 $\pm$ 1.842	0.468 $\pm$ 0.729		0.829 $\pm$ 0.673	0.369 $\pm$ 0.537	0.193 $\pm$ 0.128

Equations

$${\mathbf{x}}_{\text{min}}=\underset{\mathbf{x}\in \mathcal{X}}{\mathrm{argmin}}f(\mathbf{x}).$$

$$f(\mathbf{x}):=g(\mathbf{x},\mathbf{h}),\mathbf{h}:=h(\mathbf{x}),$$

$$p(\mathbf{H},\mathbf{\theta }\text{nonscript}|\text{nonscript}\mathbf{X},\mathbf{F})\propto p(\mathbf{F}\text{nonscript}|\text{nonscript}\mathbf{X},\mathbf{H},\mathbf{\theta })p(\mathbf{H})p(\mathbf{\theta })$$

$$\begin{array}{rl}& p({\mathbf{f}}_{\ast }\text{nonscript}|\text{nonscript}{\mathbf{x}}_{\ast },\mathbf{X},\mathbf{F})=\int p({\mathbf{f}}_{\ast }\text{nonscript}|\text{nonscript}{\mathbf{x}}_{\ast },{\mathbf{h}}_{\ast },\mathbf{X},\mathbf{F},\mathbf{H},\mathbf{\theta })\\ & p(\mathbf{H},\mathbf{\theta }\text{nonscript}|\text{nonscript}\mathbf{X},\mathbf{F})p({\mathbf{h}}_{\ast })\mathrm{d}\mathbf{H}\mathrm{d}\mathbf{\theta }\mathrm{d}{\mathbf{h}}_{\ast },\end{array}$$

$$\begin{array}{rl}& \alpha ({\mathbf{x}}_{\ast })=\int \mathbb{U}({\mathbf{f}}_{\ast },{\mathbf{x}}_{\ast },\mathbf{X},\mathbf{F})p({\mathbf{f}}_{\ast }|{\mathbf{x}}_{\ast },\mathbf{X},\mathbf{F})\mathrm{d}{\mathbf{f}}_{\ast },\end{array}$$

$$f(\mathbf{x})=g(\mathbf{x})+h,h\sim \mathcal{N}(0,\sigma (\mathbf{x}{)}^2),$$

$$f(\mathbf{x})=g(\mathbf{x},\mathbf{h}),g\sim \mathcal{G}\mathcal{P},\mathbf{h}\sim \mathcal{N}(0,\mathbb{I}).$$

$$\begin{array}{r}p(\mathbf{F}\text{nonscript}|\text{nonscript}\mathbf{X},\mathbf{H},\mathbf{\theta })=\mathcal{N}(\mathbf{F}|0,\mathbf{K}),\\ p(\mathbf{H})=\prod _{n=1}^N\mathcal{N}({\mathbf{h}}_n|0,{\sigma }_h^2\mathbb{I}),\end{array}$$

$$\alpha ({\mathbf{x}}_{\ast })$$

$$f(\mathbf{x})=g(\mathbf{x},\mathbf{h}),\mathbf{h}\sim \mathcal{N}(0,\mathbb{I}),$$

$$\text{pow10}(x)=10^x$$

$$\text{pow2int}(x)=\text{int}(2^x)$$

Full text

Machine Learning, Bayesian Optimization, Gaussian Processes, Robust, Challenging structures, ICML

1 Introduction

Bayesian optimization (BO) (Snoek et al., 2012a) is a method for finding the optimum of functions that are unknown and expensive to evaluate. By fitting a surrogate model to the samples of an unknown objective, the BO procedure iteratively picks the new samples of the objective believed to be the most informative about where the optimum is located.

Model misspecification has significant negative implications for any machine learning tasks. This is especially true for sequential decision making tasks such as BO, where the model is used not only to locate the optimum based on the collected data but also to decide where to collect data for future decisions. If the surrogate model is misspecified, it is likely to acquire samples that are less informative about the optimum, which will lead to a less efficient optimization. Therefore the quality of the surrogate model is essential to achieve both efficient and reliable results.

Many works have been done towards avoiding model misspecification in the surrogate model for BO, such as handling non-stationary objective functions with warpings (Snoek et al., 2014), tree-structured dependencies in the search space (Jenatton et al., 2017), and searching the optimum from piecewise comparisons (González et al., 2017). Comparing with the Gaussian process (GP) regression model in the standard BO setting, these methods avoid model misspecification in real-world problems by using more sophisticated surrogate models that are suitable for the corresponding problems. Bayesian inference with more sophisticated surrogate models will often require additional data to reduce uncertainty and confirm beliefs, because it considers more possibilities. Importantly the ultimate goal of BO is to find the optimum, not to model the unknown objective as precisely as possible. In practice, this means that using a surrogate with high complexity might perform worse compared to a simpler class even if the former contains the true objective function.

Instead of building a complex surrogate model with minimal model misspecification, we propose an alternative approach which allows trading off accuracy in modeling the objective with efficiency of capturing informative structures from small amounts of data. For example, we observe that structures such as local oscillations and discontinuities are less important to capture for the purposes of BO. Such details often require a lot of data to be closely captured in a surrogate model but do not help the search for the optimum, unless the search reaches the last stage of pinpointing the exact location of the optimum. To ignore these details, we associate an independent random input variable with every evaluation of the unknown function. As the random variables associated with new evaluations are conditionally independent of the posterior random variables associated with observed data given the function, this is referred to as irreducible uncertainty. Such variables are similar to the noise variables in regression models, which are used to capture measurement noise and the data variance that cannot be attributed to the input variables. In contrast to noise variables for noisy outcomes, where there is irreducible uncertainty about the data, there is now irreducible uncertainty in the model of the function.

We propose to use the surrogate models that are specified over well-behaved approximations of the objectives, which can be more useful for the search of the optimum (see Figure 1), augmented with flexible “noise” distributions to treat the nuisance parameters. We will demonstrate that, using the same function approximation, a surrogate model with a more flexible nuisance parameter distribution is more robust against challenging structures. In this paper we focus on noise-free objectives with complicated, oscillatory or discontinuous structures. In particular, we propose to use a Latent Gaussian process (LGP) (Pfingsten et al., 2006; Wang & Neal, 2012; Yousefi et al., 2016; Bodin et al., 2017) as the surrogate model due to its flexible nuisance parameter distribution and show that it outperforms the surrogate models with less flexible distributions such as GPs with additive likelihoods. LGP allows us to disentangle the complicated structures a GP surrogate struggles to model while highlighting important structures.

Our main contributions are:

• •

We propose to address challenging objective functions for BO by using a distribution in the surrogate model to explain structure that is challenging to model with few observations.

• •

We propose to use latent Gaussian processes (LGP) as surrogate models, which support non-stationary and non-Gaussian residuals.

• •

With experiments on multiple BO benchmarks, we show that our method significantly outperforms existing approaches.

2 Modulating Surrogates

Let $f:\mathcal{X}\to\mathbb{R}$ be an unknown, noise-free objective function defined on a bounded subset $\mathcal{X}\subset\mathbb{R}^{Q}$. The goal of BO is to solve the global optimization problem of finding

[TABLE]

In real world problems, the objective function is often not a well-behaved function and a suitable model is difficult to specify. Instead of applying an automated model selection method (Malkomes & Garnett, 2018), we propose to model only the essential structure of the objective function that is well-behaved and leave the rest of the function details to be absorbed in a noise distribution.

We consider the family of objective functions $f$ that can be represented as a composition of a well-behaved function and another arbitrary, latent function capturing the challenging details, i.e.

[TABLE]

where $g$ is a well-behaved function that can be nicely modeled by a surrogate model of choice, which is a Gaussian process (GP) in this paper, and the vector-valued function $h(\bm{x})$ encodes the structures which the surrogate model struggles to capture. In general this composition allows for complicated interactions between $\bm{x}$ and $\bm{h}$, producing complicated realizations of the function which is observed through data. A simple, special case of a function composition is additive structure 111Note that in the additive case, $h(\bm{x})$ must match the output in shape, i.e. be one-dimensional., i.e. $f(\bm{x})=g(\bm{x})+h(\bm{x})$.

Instead of modeling $h(\bm{x})$ as part of the surrogate model, we propose to ignore the structure of the objective function in $h(\bm{x})$ by replacing $h(\bm{x})$ with a random variable $\bm{h}$ per data point. The random variables $\bm{h}$ for different data points are independent among each other. The objective function becomes a function of two variables $g(\bm{x},\bm{h})$, in which $\bm{h}$ is a random variable which explain the data variance that cannot be explained by $\bm{x}$. In this paper, we use a normal distribution for the prior of $\bm{h}$, $\bm{h}\sim\mathcal{N}(0,\mathbb{I})$. Note that, although the distribution of $h(\bm{x})$ induced by the data distribution for $\bm{x}$ may not be zero-mean and unit-variance, it is easy to reformulate it as a linear transformation of a normal distribution with zero-mean and unit-variance and the resulting linear transformation can be absorbed into the function $g$. For further details on the definition, see the supplement.

With the above formulation, a BO method can be developed by constructing a model of the well-behaved function $g$ and a model of $h$. At each step of the BO optimization, a set of input and output pairs of the objective function has been collected, denoted as $\bm{X}=(\bm{x}_{1},\ldots,\bm{x}_{N})^{\top}$ and $\bm{F}=(\bm{f}_{1},\ldots,\bm{f}_{N})^{\top}$. The output $\bm{F}$ denotes the noise-free observations of the objective function. The Bayesian inference of the model aims at inferring the posterior distribution

[TABLE]

where $\bm{\theta}$ are the hyperparameters of the surrogate model and $\bm{H}=(\bm{h}_{1},\ldots,\bm{h}_{N})^{\top}$ is the concatenation of the nuisance parameters associated with the individual data points. The location of the next evaluation is determined according to an acquisition function, which uses the predictive distribution $\operatorname{p}(\bm{f}_{*}\nonscript\mkern 2.0mu|\nonscript\mkern 2.0mu\mathopen{}\allowbreak\bm{x}_{*},\bm{X},\bm{F})$ of the surrogate model,

[TABLE]

where $\bm{x}_{*}$ is the input of the prediction and $\bm{f}_{*}$ is the noise-free observation at the location $\bm{x}_{*}$. The predictive distribution of the latent variable $\operatorname{p}(\bm{h}_{*})$ associated with new evaluations is as of the i.i.d. assumption equal to the prior. As such $\operatorname{p}(\bm{h}_{*})$ contains model uncertainty irreducible by the active sampling loop, which we suggest to ignore via augmentation, see the supplement for details.

With the predictive distribution Eq. 4, the expectation of the acquisition function is derived as

[TABLE]

where the acquisition function of choice is denoted $\mathbb{U}$. Note that the predictive distribution due to the marginalization over $\bm{H}$ and $\bm{\theta}$ generally has a complicated form and that the above integral often requires approximate methods.

3 Latent GP surrogates and other choices

In the previous section we presented the BO formulation. We will now proceed to implement the formulation, and address the choice of surrogate model for the function (Eq. 2).

Additive noise model. As briefly mentioned in the previous section, a simple case of the composition (Eq. 2) is an additive structure, $f(\bm{x})=g(\bm{x})+h(\bm{x})$. Following the process of replacing $h(\bm{x})$ with the random variable $\bm{h}$, the resulting surrogate model of the objective function is

[TABLE]

where the variance of $h$ is assumed to be $\sigma^{2}$ in order to adapt to the value range of $f$. With a GP surrogate model for $g$, the above model recovers the GP regression model with the Gaussian likelihood.

A typical choice in the above model is to assume $\sigma^{2}$ to be constant, leading to a homoscedastic model. A limitation of noise variances being the same across all the datapoints is that it limits the capability of the model in terms of absorbing irregular variance. A straight-forward extension of the above model is the GP with heteroscedastic noise, in which the noise variance $\sigma^{2}$ is allowed to be different among data points (Goldberg et al., 1998; Lázaro-Gredilla & Titsias, 2011). Another choice could be specifying a GP prior for $h(\bm{x})$ and thus recover an additive GP model for $f$ (Bernardo et al., 1998; Duvenaud et al., 2011). Other available choices for an additive noise model include Student’s t-distribution (Jylänki et al., 2011), Laplace (Kuss, 2006) or mixture of Gaussian likelihoods as (Kuss, 2006; Stegle et al., 2008; Naish-Guzman & Holden, 2008) where (Naish-Guzman & Holden, 2008) considers the heteroscedastic case.

Latent Gaussian process. A major limitation of the additive noise models in general is the inability to capture the interaction between the input $\bm{x}$ and the noise $\bm{h}$. Another choice that produces a more flexible noise distribution is to introduce additive noise in the input of a GP (McHutchon & Rasmussen, 2011; Girard et al., 2003; Girard, 2004). This would correspond to the case of $f(\bm{x})=g(\bm{x}+\bm{h})$. A further more general case of the proposed methodology is to allow non-linear interactions between the random variable $\bm{h}$ and $\bm{x}$. This can be formulated as

[TABLE]

This formulation aligns with the general assumptions proposed in the previous section. In particular, the well-behaved function $g$ is assumed to follow a GP prior distribution, and the random variable derived from the challenging details of the objective function $\bm{h}$ feeds directly into the GP surrogate model. This allows for an arbitrary interaction between $\bm{h}$ and $\bm{x}$, as specified by the covariance function. The introduction of the random variable $\bm{h}$ in the input results in a flexible noise distribution, as the GP model can warp the normal distribution of $\bm{h}$ into a sophisticated distribution and allow non-linear interactions between $\bm{h}$ and $\bm{x}$. This GP model in (7) is also known as a latent Gaussian process (LGP) (Pfingsten et al., 2006; Wang & Neal, 2012; Yousefi et al., 2016; Bodin et al., 2017), which is developed for regression with heteroscedastic noise and non-Gaussian residuals. The non-Gaussian marginals arise as a consequence of the latent covariates and their nonlinear transformation through the covariance function.

Function modulation via $\mathcal{H}$

If we assume a stationary kernel over the product space $\mathcal{X}\times\mathcal{H}$, a constant $\bm{h}_{n}$ for all observations can be interpreted as the $\mathcal{H}$ subspace having no influence. This is due to the stationary property of the kernel, where covariances are determined only by the distances between points.

With everything else held constant, if an observation is moved away from other observations in the $\mathcal{H}$ space, the covariances between that observation and the others are reduced. Similarly, if the length scale in $\mathcal{X}$-direction is shortened, the covariances between that observation and the others can be equally reduced, but that also reduces the covariances between all other observations due to the global influence of the hyperparameter.

Structures in the data could be explained solely by reducing the $\mathcal{X}$-direction length scale adequately. In that case, evaluating the posterior at $\bm{h}_{\ast}=\bm{0}$ would yield exactly the posterior of a standard GP. Conversely, structures could be explained solely as observations being adequately far from each other in the $\mathcal{H}$ space while maintaining a longer $\mathcal{X}$ length scale. Evaluating the posterior at $\bm{h}_{\ast}=\bm{0}$ then yields a posterior that is both influenced by the longer length scale and which has lower covariances with the data, effectively producing a posterior over smoother functions. If the posterior inputs $\bm{H}$ are sufficiently far away with respect to the $\mathcal{H}$-direction length scale, all data variation will be captured in $\mathcal{H}$ and the posterior of the function at $\bm{h}_{\ast}=\bm{0}$ will in effect ignore the data.

The posterior weighting over this range of solutions is determined by the trade-off between the GP function prior and the prior of the latent inputs. As such, by controlling this trade-off, we can control properties of structures to be ignored and the ones to be used for search (see Figure 2). Important to note is that the mentioned data ignorance effect affects individual data points via the posterior of the corresponding latent input $\bm{h}_{n}$, which is influenced by the local and global fit of the function prior.

Reparameterization of LGP for ease of specifying the modulation prior

In a BO setting, some prior knowledge about what constitutes a significant change in the input space is often available. We would like to specify a joint prior of the GP and the latent inputs to ignore structures at the appropriate scale. In order to do this, we (re)-parametererise it in the following way. We set the lengthscale in the $\mathcal{H}$-direction to be the same as in the $\mathcal{X}$-direction and parameterize the latent input prior as $\mathcal{N}(\bm{h}_{n}|\bm{0},\sigma_{h}^{2}\mathbb{I})$ instead of a unit Gaussian. There is an equivalence between parameterizing $\sigma_{h}$ or setting this trade-off via a separate lengthscale for $\mathcal{H}$ as in (Wang & Neal, 2012). However, by using the above parameterization there is a direct correspondence in covariance reduction from moving an observation in $\mathcal{H}$ as in $\mathcal{X}$, and the prior for the latent inputs can be interpreted as a prior over the coarseness of the function we wish to exploit for search. As such, intuitions about the scale in $\mathcal{X}$ directly translate into the parameterization of the prior. See Figure 2 for a visualization of how changing $\sigma_{h}$ in the prior affects the modulated function posterior. To make the $\sigma_{h}$ parameterization relevant across input sizes and dimensionalities, we rescale the input domain $\mathcal{X}$ to be a unit hyper-cube and set $\sigma_{h}$ proportionally to the length of the diagonal of the domain $\sqrt{Q}$ (where $Q$ is the number of dimensions).

Posterior inference and acquisition calculation. In BO we assume that $N$ pairs of inputs and outputs $\bm{X}=(\bm{x}_{1},\ldots,\bm{x}_{N})^{\top}$ and $\bm{F}=(\bm{f}_{1},\ldots,\bm{f}_{N})^{\top}$ have been collected. To suggest the location for the next evaluation, we first need to infer the posterior distribution of the latent variables, which are $\bm{H}$ and $\bm{\theta}$ in LGP, and then search for the maximum of the acquisition function $\alpha(\bm{x})$.

Given the observed data, the probabilistic model of LGP is formulated as

[TABLE]

where $\bm{K}$ is the covariance matrix computed using a chosen kernel function $k(\cdot,\cdot)$ over the set of data points $\{\bar{\bm{x}}_{n}\}_{n=1}^{N}$, and $\bar{\bm{x}}_{n}$ is the concatenation of two vectors $(\bm{x}_{n}^{\top},\bm{h}_{n}^{\top})^{\top}$. Because $\bm{h}_{n}^{\top}$ enters the kernel function non-linearly, it is clear that the posterior distribution $\operatorname{p}(\bm{H},\bm{\theta}\nonscript\mkern 2.0mu|\nonscript\mkern 2.0mu\mathopen{}\allowbreak\bm{F},\bm{X})$ is intractable. To ensure the quality of the acquisition function, usually, BO methods draw posterior samples of latent variables via Markov Chain Monte Carlo (MCMC) methods such as slice sampling (Snoek et al., 2012a) or Hamiltonian Monte Carlo (Duane et al., 1987). We follow this practice and provide details in the supplement.

With the approximate posterior samples $\{\bm{H}_{i},\bm{\theta}_{i}\}_{i=1}^{M}$, we approximate the acquisition function with LGP in (LABEL:eq:utility) with Monte Carlo samples,

[TABLE]

where $\hat{\alpha}(\bm{x}_{*},\bm{H}_{i},\bm{\theta}_{i})$ is the acquisition function given the latent variables of LGP, which is closed-form for common acquisition functions such as expected improvement (EI) and upper confidence bound (UCB).

4 Related Work

Performing BO on an objective function that is not well-behaved is very challenging. Our method takes a Bayesian approach by incorporating a flexible noise distribution and utilising Bayesian inference to assign the challenging details of the objective function to the noise distribution. An alternative approach to this problem is to perform a model selection for the surrogate model, such that the choice of the surrogate model becomes a trade-off between the complexity of the model and the ability to locate the optimum under limited data, which has been explored in (Malkomes & Garnett, 2018). The approach uses a compositional kernel grammar from (Duvenaud et al., 2013) to induce a space of GP models to choose from. Although this and other model selection procedures (Malkomes et al., 2016; Duvenaud et al., 2013; Grosse et al., 2012; Gardner et al., 2017) themselves have shown promise, in addition to the computational overhead, the procedures are still reliant on the existence of suitable models in this space. It remains challenging to handle cases where the objective function contains structure that is both hard to specify a-priori, and that is unhelpful in guiding the search to the optimum.

The idea of making use of noise models for dealing with model mismatch to noise-free data is not in itself new. In (Gramacy & Lee, 2012) it was shown that introducing noise in the modelling of noise-free computer experiments can lead to models with better statistical properties such as predictive accuracy and coverage. In that work, homoscedatic noise was addressed and used in a regression context.

In this paper we consider noise-free functions and address model misspecification of the function surrogate, but many works have been done to make BO resilient to noisy experiments. For example, robust noise distributions such as Student’s t-distribution have been used to make BO more resilient to noise outliers (Martinez-Cantin et al., 2017a, b). Approaches to address noisy experiments, via the addition of likelihood functions, can be combined with our approach.

Hierarchical surrogate models with input warpings have been proposed to tackle BO for non-stationary objective functions (Snoek et al., 2014; Oh et al., 2018; Calandra et al., 2016). A particularly successful application is hyperparameter optimization for machine learning methods, in which the parameters are often presented in logarithmic scales. In this case, the Beta cumulative distribution function, which only has two parameters, serves as a good warping function (Snoek et al., 2014). Such augmentation in surrogate models requires strong domain knowledge of the objective function, and one often still has to control the increased complexity of the surrogate model, which is orthogonal to our approach.

5 Experiments

In this section we will demonstrate the benefit of our approach empirically. As the approach is motivated by robustness to the presence of challenging structures in the objective function, we will test its ability to improve search efficiency on a range of functions exhibiting such structure. Visual examples of functions with typical properties are shown in Figure 3. As we will show, our approach increases reliability in the search when faced with detrimental structure (see Figure 4) that has a large negative impact on traditional surrogates.

Baselines and metric

We compare with and without function modulation (Section 2) - implemented as in Section 3 - on a popular GP model setup for BO. In addition, we compare the LGP against other methods of handling challenging structure in the objective function, namely (i) a noiseless GP, (ii) a GP with homoscedastic noise, (iii) a GP with heteroscedastic noise and (iv) a non-stationary, Warped GP (Snoek et al., 2014).

We follow the standard practice to compare across benchmarks and provide the mean gap estimated over 20 runs as in (Malkomes & Garnett, 2018). The gap measure is defined as $\frac{f(x_{\text{first}})-f(x_{\text{best}})}{f(x_{\text{first}})-f(x_{\text{optimum}})}$, where $f(x_{\text{first}})$ is the minimum function value among the first initial random points, $f(x_{\text{best}})$ is the best function value found within the evaluation budget and $f(x_{\text{optimum}})$ is the function’s true optimum. Methods are judged to have very similar or equivalent performance to the best performing if not significantly different, determined by a two-sided paired Wilcoxon signed-rank test at 5% significance (Malkomes & Garnett, 2018). We also report regret (with mean and standard deviation) in the supplement.

We use the Matérn 5/2 kernel for all surrogates, the expected improvement acquisition function (where not otherwise stated) and Bayesian hyperparameter marginalisation as in (Snoek et al., 2012a). For the maximization of the expected utility with respect to input location, we use $\delta$-cover sampling, as in (De Freitas et al., 2012). The Warped GP implementation and inference is from the Spearmint package (Snoek et al., 2014). For further details, we refer to the supplement.

Benchmark datasets

We perform the comparisons on benchmarks from (McCourt, 2016; Head et al., 2018) using the default domains provided by respective benchmark, detailed in the supplement. In addition, problems are marked with the descriptive properties given in (McCourt, 2016) and in the supplement that can reflect the relative difficulty of the task.

Priors on the latent input variables

The prior $\operatorname{p}(\bm{h}_{n})=\operatorname{\mathcal{N}}(\bm{0},\sigma_{h}^{2}\mathbb{I})$ can be parameterized in relation to the relative scale of the characteristics to be ignored. We specify the function prior over the product space $\mathcal{X}\times\mathcal{H}$ using a kernel with common parameters for $\bm{x}_{n}$ and $\bm{h}_{n}$. Thus, the standard deviation of the prior $\sigma_{h}$ relates directly to distances in the $\mathcal{X}$-direction. When domain-specific knowledge is available, $\operatorname{p}(\bm{H})$ may be specified at an appropriate scale. However, we often do not have access to such knowledge. In all our experiments, we adopt a hierarchical prior approach whereby $\sigma_{h}$ is sampled uniformly from a small candidate set at each evaluation. Specifically, $\sigma_{h}\sim\operatorname{\mathbb{U}}(\{0.1d,0.01d,0\})$ where $d=\sqrt{Q}$, the length of the diagonal of the unit Q-dimensional hypercube. We found that this approach performed well empirically and is applied consistently across all our experiments where not otherwise specified. A choice of $\sigma_{h}\to 0$ corresponds to a noiseless GP without latent covariates.

Evaluation on benchmark suite

Table 1 presents results across a wide range of benchmark functions consisting of the SigOpt benchmark suite (McCourt, 2016). Three additional real-world benchmarks (Head et al., 2018; Malkomes & Garnett, 2018; Kaelbling & Lozano-Pérez, 2017) are included in the bottom section of the table. The benchmarks from (McCourt, 2016) are popular functions used in both black-box optimization as well as classic optimization literature. As of the focus of the paper, benchmarks from the literature exhibiting challenging properties such as oscillatory local structures were included, in addition to simpler functions for reference.

In general, the noise-free, homoscedastic and heteroscedastic noise GPs tend to either share best place with the LGP or be outperformed by it. The Warped GP, which warps the input space to obtain a tight fit to the data, consistently struggle with the complicated and oscillatory benchmarks. On some benchmarks there are large differences in favour of using noisy surrogates on the noise-free benchmarks. Such an example is Ackley 6D, which in the dataset is described as “technically oscillatory, but with such a short wavelength that its behavior probably seems more like noise” (McCourt, 2016). Another example is the Shubert function, which has multiple sharp local optima surrounded by large oscillations. On 2 of 18 benchmarks, the LGP was not best (or within the two-sided Wilcoxon test), but instead the homoscedastic GP. These functions were Weierstrass, which has a homoscedastic characteristic (see Figure 3), and Deflected Corrugated Spring, on both of which the LGP obtained the second highest mean gap. In contrast to the LGP and the heteroscedastic GP, the noise model of the homoscedastic GP sometimes hurt performance in relation to the noise-free GP. Given the black-box nature of functions in BO, it is important that the surrogate noise model ’turns off’ adequately when not needed. The heteroscedastic GP provided significant benefit on two benchmarks over the GP, whereas the LGP provided such benefit on eight benchmarks.

Real-world

Apart from widely used synthetic functions, we also compare our method on three real-world problems. The results are shown in the last three rows of Table 1. One of the benchmarks is an active learning task of a robot pushing a box within a simulation. As we show in Figure 5, the benchmark’s response surface is both nonsmooth and oscillatory. The LGP reliably found good solutions on the benchmark, while the other surrogates sometimes failed, resulting in high variances. The homoscedastic GP performed the worst, which we suspect is due to the nonsmooth and heavily oscillatory structures forces a high global noise level, which may lead to failure in utilising informative structure in other regions.

Other aquisition functions

We suggest that the problem with structures challenging to model is relevant to address irrespective of the acquisition function. To confirm that the method is applicable also using other acquisition functions we ran the Corrupted Exponential benchmark using both Expected Improvement (EI) and Lower Confidence Bound (LCB) with the default exploration weight ($=2.0$) from GPyOpt (GPyOpt, 2016). As can be seen in Table 1, in the case of EI, the GP and the heteroscedastic GP performed worse than the homoscedastic GP and the LGP. The homoscedastic GP achieved the highest mean gap, but the difference was not significant under the Wilcoxon test to the LGP which obtained a similar mean gap. Using the LCB acquisition function the performance for the homescedastic GP decreased to $0.818$ and the LGP increased to $0.858$, and their difference in rank in favour of the LGP was significant under the test. The heteroscedastic GP increased using LCB to $0.797$ and the GP to $0.751$, remaining as worst performing.

As the experimental evaluation demonstrates, our suggested approach for handling challenging structures in the objective function consistently improved reliability and performance over the traditional surrogate on a wide range of benchmarks. Importantly, on benchmarks where the extended methodology were not needed the performance aligned with that of the traditional surrogate. When it was needed, it was shown to often have a large positive impact on overall efficiency of the search.

6 Conclusion

We have presented an approach to Bayesian Optimization where the surrogate model is alleviated from needing to explain the observed objective function values perfectly, which is challenging for complicated or nonsmooth functions. Instead, we model the essential structure of the objective function that is well-behaved and leave the rest of the function details to be absorbed in a noise distribution. We show experimentally how our approach is able to solve synthetic and real-world benchmarks with challenging local structures reliably. Importantly our methodology can be applied to any surrogate model used for BO, and the specific case addressed in the paper can be included in any Gaussian process-based surrogate.

7 Acknowledgements

This project was supported by Engineering and Physical Sciences Research Council (EPSRC) Doctoral Training Partnership, the Hans Werthén Fund at The Royal Swedish Academy of Engineering Sciences, EPSRC CDE (EP/L016540/1), EPSRC CAMERA project (EP/M023281/1), the German Federal Ministry of Education and Research (project 01 IS 18049 A), and the Royal Society.

Appendix A Appendix

A.1 $\operatorname{p}(\bm{h}_{*})$ augmentation

The predictive distribution of the latent variable for new evaluations $\operatorname{p}(\bm{h}_{*})$, even at observed inputs locations, will always by the i.i.d. assumption be equal to the prior. As such, this source of uncertainty cannot be reduced by active sampling, nor does it reflect observational stochasticity which is made clear by the noiseless experiments. Including it would simply include unhelpful artifacts in the decision upon where to collect data and in the worst case the search would be stuck, see Figure 6. The importance of special treatment of similar uncertainty within active sampling has been noted in (González et al., 2017).

Noting the risk of ill-effects of simply marginalizing $\operatorname{p}(\bm{h}_{*})$, one might be tempted to introduce statistical dependencies in the model such that the belief about $\bm{h}_{*}$ associated with a new evaluation is updated from e.g. neighbouring observations. However, such dependencies does not come for free, as they would by necessity limit the model’s capability to explain away via $\bm{H}$. For example, by constraining nuisance parameters to be similar within neighbourhoods of $\mathcal{X}$ the ability to be robust to discontinuities such as step functions would be reduced.

We suggest to remove the effect $\bm{h}_{*}$ has on the predictive distribution of $\bm{f}_{*}$. In the case of additive models this is easily achieved. For example, in the GP regression context with noisy observations the predictive distribution of the noise-free latent function is easily derived by removal of the noise variable after the noise has been considered in the data. The reason why it is easy in the additive case is because the expectation is trivially separable, as $\mathbb{E}[g+h]=\mathbb{E}[g]+\mathbb{E}[h]$. In this paper we consider the introduction of nuisance parameters in surrogates generally, including those with non-linear interactions. In the general non-linear case how to remove the interaction is an open question. For these cases, when there is no available model-specific treatment, we propose the heuristic of collapsing the prior distribution of $\bm{h}_{*}$ into a Dirac delta distribution centered at its mode $\bm{h}_{*}=\bm{0}$, which is consistent with the additive case.

A.2 Further details on the objective function definition

[TABLE]

where $g$ is a well-behaved function that can be nicely modeled by a surrogate model of choice, which is a Gaussian process (GP) in this paper, and the vector-valued function $h(\bm{x})$ encodes the structures which the surrogate model struggles to capture. Let’s assume an uninformative prior distribution of $\bm{x}$ over $\mathcal{X}$, e.g. a uniform distribution $\bm{x}\sim\mathcal{U}(\mathcal{X})$. We denote the mean and variance of $h(\bm{x})$ under the prior distribution as $\mu_{\bm{h}}=\mathbb{E}_{\operatorname{p}(\bm{x})}[h(\bm{x})]$ and $\Sigma_{\bm{h}}=\mathbb{E}_{\operatorname{p}(\bm{x})}[(h(\bm{x})-\mu_{\bm{h}})(h(\bm{x})-\mu_{\bm{h}})^{\top}]$, respectively. An important step to convert $\bm{h}$ into being part of a noise distribution is to treat it as random and independent of $\bm{x}$, i.e. $\bm{h}$ becomes an independent random variable with respect to each data point, just like a standard noise term. In this paper, we use a normal distribution for $\bm{h}$, $\bm{h}\sim\mathcal{N}(\mu_{\bm{h}},\Sigma_{\bm{h}})$. In this way, the objective function becomes a function of two variables $g(\bm{x},\bm{h})$, in which $\bm{h}$ is a random variable independent of $\bm{x}$ and which explain the data variance that cannot be explained by $\bm{x}$. Note that the random variable $\bm{h}$ can be equivalently rewritten as $\bm{h}=\mu_{\bm{h}}+\bm{L}\bar{\bm{h}},\,\Sigma_{\bm{h}}=\bm{L}\bm{L}^{\top},\,\bar{\bm{h}}\sim\mathcal{N}(\bm{0},\mathbb{I})$ and then the objective function becomes $g(\bm{x},\mu_{\bm{h}}+\bm{L}\bar{\bm{h}})$. As $g$ is a non-linear function inferred during BO, the linear transform can be absorbed into the formulation.

A.3 Step function example

A.4 Model setup, inference and auxiliary optimization

We use the Warped GP (Snoek et al., 2014) surrogate model with default settings as provided by the Spearmint library (spe, ). We use the Matérn 5/2 kernel for all surrogates and the model-marginalised expected improvement acquisition function as in (Snoek et al., 2012b). For all baselines we make us of hyper-parameter marginalisation via Markov Chain Monte Carlo (MCMC) (Shahriari et al., 2016).

For the noise-free, homoscedastic GP and Warped GP (all with few parameters) we use slice sampling, as recommended for BO in (Snoek et al., 2012a) due to its automatic adjustment of the step size to match the local shape of the density function. The heteroscedastic GP has a latent noise variance per observation. Similarly, LGP has the latent inputs $\bm{H}$ associated with the observations, making inference impractical for both of these models using slice sampling as of comparably large dimensionalites. For these models, in all BO experiments, we use Hamiltonian Monte Carlo with step size adaptation. A burn-in of $30,000$ steps and a thinning rate of $x50$ (select every 50th value) were used, and $100$ posterior samples collected. Step-size adaptation were made during $80\%$ of the burn-in phase, with a target acceptance probabilty of $0.75$ which is the center of asymptotically optimal rate for HMC (Betancourt et al., 2014).

Observed outputs are normalized to have zero mean and unit variance (standard normalization) at each iteration. For the GP surrogate with and without the noise models as well using the latent input extension, we use a LogNormal(0, 1) prior for the lengthscale and noise variance parameters and use unit signal variance. We use one dimension for the latent variables $\bm{h}_{n}$ in the comparisons with the baselines. For the maximization of the expected utility with respect to input location, we use $\delta$-cover sampling, as in (De Freitas et al., 2012), for all models (where in our case the expected utility is the model-marginalised expected improvement). The sampling scheme works by iteratively sampling the utility more densely in $\mathcal{X}$ around the location of current highest obtained utility. Specifically, we double the concentration of an uniform sample density at each auxiliary iteration by multiplication of each side of the current sampling hypercube by a factor of $2^{-\nicefrac{{1}}{{Q}}}$ for $30$ iterations. Without loss of generality, all function domains are re-scaled to unit hypercubes for a consistent parameterization of priors across functions. In all experiments we start from $2$ uniformly drawn initial observations and stop after $50$ and $100$ observations, respectively, as indicated in the experiments.

A.5 Benchmarks

A.5.1 Property labels

For benchmark property labels, see Table 2.

A.5.2 Domains

[TABLE]

[TABLE]

For benchmark domains, see Table 3.

A.5.3 Regret version of results

For regret version of results table, see Table 6.

A.5.4 Corrupted Holder Table and Corrupted Exponential

The following corruption functions in Figure 8 was used for the benchmarks Corrupted Holder Table and Corrupted Exponential. small_corruption_func and large_corruption_func was applied to the SigOpt benchmarks (McCourt, 2016) Holder Table 2D and Exponential 8D, respectively. The input dimensions are re-scaled to be between 0 and 1 before the corruption is applied. The new function minimum of each function (due to the corruption) was estimated via $1e6$ uniformly drawn samples in the domains.