Bayesian Search for Robust Optima

TL;DR

This paper introduces a Bayesian optimisation method for robust black-box function minimisation, focusing on finding regions of insensitivity with respect to inputs, and demonstrates superior performance over existing methods in tests and real-world robot tasks.

Contribution

It proposes a novel robust optimisation approach using Gaussian processes and evolutionary algorithms to identify insensitive regions, improving stability and performance in black-box problems.

Findings

Outperforms state-of-the-art methods on benchmark functions.

Efficiently locates insensitive regions with high expected improvement.

Successfully applied to real-world robot arm control tasks.

Abstract

Many expensive black-box optimisation problems are sensitive to their inputs. In these problems it makes more sense to locate a region of good designs, than a single-possibly fragile-optimal design. Expensive black-box functions can be optimised effectively with Bayesian optimisation, where a Gaussian process is a popular choice as a prior over the expensive function. We propose a method for robust optimisation using Bayesian optimisation to find a region of design space in which the expensive function's performance is relatively insensitive to the inputs whilst retaining a good quality. This is achieved by sampling realisations from a Gaussian process that is modelling the expensive function, and evaluating the improvement for each realisation. The expectation of these improvements can be optimised cheaply with an evolutionary algorithm to determine the next location at which to evaluate the expensive function. We describe an efficient process to locate the optimum expected improvement. We show empirically that evaluating the expensive function at the location in the candidate uncertainty region about which the model is most uncertain, or at random, yield the best convergence in contrast to exploitative schemes. We illustrate our method on six test functions in two, five, and ten dimensions, and demonstrate that it is able to outperform two state-of-the-art approaches from the literature. We also demonstrate our method one two real-world problems in 4 and 8 dimensions, which involve training robot arms to push objects onto targets.

Figures

Tables

Table 1. Table 1. Summary of the benchmark functions used for experimentation, the domains 𝒳 𝒳 \mathcal{X} over which they were evaluated. A † † \dagger signifies that the stated equation has been modified from the referenced source (see Section 4 for full explanation).

Name	Equation	$𝒳$
Bumped Bowl† (Mirjalili and Lewis, 2016)	$f_1\left(𝐱\right)=\log {\Vert 𝐱\Vert }^2+e^{-10{\Vert 𝐱\Vert }^2}$	${[-4,4]}^D$
Levy03 (Laguna and Martí, 2005)		${[-4,4]}^D$
Styblinski-Tang (Styblinski and Tang, 1990)	$f_3\left(𝐱\right)=\frac{1}{2}{\sum }_{d=1}^D\left(x_d^4-16x_d^2+5x_d\right)$	${[-5,5]}^D$
Robust Problem 4 (Mirjalili and Lewis, 2016)	$f_4\left(𝐱\right)=1.3-\frac{1}{d}{\sum }_{d=1}^{D}H\left(x_d\right);$	${[-2,2]}^D$
Stepped Sphere† (Mirjalili and Lewis, 2016)	$f_5\left(𝐱\right)=D-D{\prod }_{d=1}^{D}G\left(x_d\right)+\frac{1}{100}{\sum }_{d=1}^D{x}_d^2;$	${[-10,10]}^D$
Quintic (Al-Roomi, 2015)	$f_6\left(𝐱\right)=\sum _{d=1}^D\left(x_d^5-3x_d^4+4x_d^3+2x_d^2-10x_d-4\right)$	${[-10,10]}^D$

Equations

$$I(\mathbf{x})=\max \left(f^{\star }-\hat{f}(\mathbf{x}),0\right),$$

$$f^{\star }=\underset{{\mathbf{x}}_n\in \mathcal{D}}{\min }f({\mathbf{x}}_n)=\underset{n}{\min }{f}_n$$

$$EI(\mathbf{x};\mathcal{D})={\int }_{-\infty }^{\infty }I(\mathbf{x})p(\hat{f}|\mathbf{x},\mathcal{D})\text{d}\hat{f}(\mathbf{x}).$$

$$EI(\mathbf{x};\mathcal{D})=(f^{\star }-\mu (\mathbf{x}))\Phi (Z)+\sigma (\mathbf{x})\varphi (Z),$$

$$p(\hat{f}|\mathbf{x},\mathcal{D})=\mathcal{N}(\hat{f}|\mu (\mathbf{x}),{\sigma }^2(\mathbf{x})),$$

$$\mu (\mathbf{x})$$

$${\sigma }^2(\mathbf{x})$$

$$\log p(\mathcal{D}|\eta )=-\frac{1}{2}\log |\mathbf{K}|-\frac{1}{2}{\mathbf{f}}^T{\mathbf{K}}^{-1}\mathbf{f}-\frac{N}{2}\log (2\pi ).$$

$${\Delta }_{ϵ}\triangleq \{\mathbf{\delta }\in \mathcal{X}|d(\mathbf{\delta })\le ϵ\}$$

$$Q_{avg}\left(\mathbf{x},{\Delta }_{ϵ}\right)={\int }_{{\Delta }_{ϵ}}f(\mathbf{x}+\mathbf{\delta })\pi (\mathbf{x}+\mathbf{\delta })\text{d}\mathbf{\delta },$$

$$Q_{max}\left(\mathbf{x},{\Delta }_{ϵ}\right)=\underset{\mathbf{\delta }\in {\Delta }_{ϵ}}{\sup }f(\mathbf{x}+\mathbf{\delta }).$$

$$\underset{\mathbf{x}\subset \mathcal{X}}{\min }Q\left(\mathbf{x},{\Delta }_{ϵ}\right).$$

$${\mathbf{x}}^{\prime }=\underset{\mathbf{x}\in \mathcal{X}}{\mathrm{argmin}}\underset{\mathbf{\delta }\in {\Delta }_{ϵ}}{\max }\text{lcb}(\mathbf{x}+\mathbf{\delta }),$$

$$\text{lcb}(\mathbf{x})=\mu (\mathbf{x})-\sigma (\mathbf{x})$$

$${\mathbf{x}}_{new}=\underset{\mathbf{\delta }\in {\Delta }_{ϵ}}{\mathrm{argmax}}\text{ucb}(\mathbf{x}+\mathbf{\delta }).$$

$${\hat{Q}}_{max}(F)=\underset{F_i\in F}{\max }{F}_i$$

$${\hat{Q}}_{avg}(F)=\frac{1}{|F|}\sum _{F_i\in F}{F}_i.$$

$${\Delta }_{ϵ}(\mathbf{x})=\left\{\mathbf{x}+\mathbf{\delta }|\mathbf{\delta }\in {\Delta }_{ϵ}\right\}$$

$${\mathbf{x}}^{\star }=\underset{\mathbf{x}\in \mathcal{N}(X_N)}{\mathrm{argmin}}\hat{Q}(\{\mu (\mathbf{x}+\mathbf{\delta })|\mathbf{\delta }\in T\}),$$

$$f(x)=\sin (3\pi x^3)-\sin (8\pi x^3),$$

$$w({\mathbf{x}}^{\prime })={\mathbf{x}}^{\prime }.$$

$$w_{{\sigma }^2}({\mathbf{x}}^{\prime })=\underset{\mathbf{x}\in {\Delta }_{ϵ}({\mathbf{x}}^{\prime })}{\mathrm{argmax}}{\sigma }^2(\mathbf{x}),$$

$$w_{\mu }({\mathbf{x}}^{\prime })=\underset{\mathbf{x}\in {\Delta }_{ϵ}({\mathbf{x}}^{\prime })}{\mathrm{argmax}}\mu (\mathbf{x}),$$

$$w_{\mathcal{U}}$$

$$w_{UCB}$$

$$\bar{\mathcal{N}}(X_N)=\bigcup _{\mathbf{x}\in \mathcal{N}(X_N)}{\Delta }_{ϵ}(\mathbf{x}).$$

$${\mathbf{x}}^{\star }=\underset{\mathbf{x}\in \mathcal{N}(X_N)}{\mathrm{argmin}}\hat{Q}\left({\Delta }_{ϵ}(\mathbf{x})\right),$$

$$\hat{Q}({\Delta }_{ϵ}(\mathbf{x}))=\underset{{\mathbf{x}}^{\prime }\in {\Delta }_{ϵ}(\mathbf{x})}{\max }\mu ({\mathbf{x}}^{\prime }).$$

$$I_m({\Delta }_{ϵ}(\mathbf{x}))=\max \left(0,Q_m^{\star }-\underset{{\mathbf{x}}_k\in {\Delta }_{ϵ}(\mathbf{x})}{\max }{f}_m({\mathbf{x}}_k)\right)$$

$$Q_m^{\star }=\underset{\mathbf{x}\in {\Delta }_{ϵ}({\mathbf{x}}^{\star })}{\max }{f}_m(\mathbf{x}).$$

Full text

Expensive black-box functions can be optimised effectively with Bayesian optimisation, where a Gaussian process is a popular choice as a prior over the expensive function. We propose a method for robust optimisation using Bayesian optimisation to find a region of design space in which the expensive function’s performance is relatively insensitive to the inputs whilst retaining a good quality. This is achieved by sampling realisations from a Gaussian process that is modelling the expensive function, and evaluating the improvement for each realisation. The expectation of these improvements can be optimised cheaply with an evolutionary algorithm to determine the next location at which to evaluate the expensive function. We describe an efficient process to locate the optimum expected improvement. We show empirically that evaluating the expensive function at the location in the candidate uncertainty region about which the model is most uncertain, or at random, yield the best convergence in contrast to exploitative schemes.

We illustrate our method on six test functions in two, five, and ten dimensions, and demonstrate that it is able to outperform two state-of-the-art approaches from the literature. We also demonstrate our method one two real-world problems in 4 and 8 dimensions, which involve training robot arms to push objects onto targets.

Robust optimisation, Bayesian optimisation, Gaussian processes

††copyright: acmcopyright††journalyear: 2021††doi: 00.0000/0000000.0000000††journal: JACM††journalvolume: 0††journalnumber: 0††article: 000††publicationmonth: 0††ccs: Computing methodologies Gaussian processes

1. Introduction

Optimisation is the search for the best-performing design with respect to a predefined objective function. In an ideal scenario the objective function would be well behaved and insensitive to small changes in the design parameters. Any loss of performance due to the mis-specification of the parameters or the inaccuracy of the model itself would be negligible and largely go unnoticed. But, for most real-world problems this is not the case. Often the landscape of real objective functions varies rapidly as the design parameters change, so even a small perturbation to these parameters could lead to an unacceptably diminished performance. Such perturbations are frequently the manifestation of uncertainties in the design process, and arise for a number of reasons. For example: tolerances in the manufacturing process may mean that the realised design is slightly different from the modelled one; once the design has been produced it may be operated away from its design conditions; the model used to optimise the design may not accurately reflect reality; or there may be uncertainties in the environmental conditions (e.g. air temperature or pressure).

The goal of robust optimisation is to locate the best designs that have stable performance irrespective of small perturbations to the design parameters. Classic global optimisers are often ineffective at solving this problem, because robust optima do not necessarily coincide with global optima; in fact optimal robust designs potentially exist in a different region of the domain altogether!

A straightforward way to quantify the robustness of a proposed optimum is to evaluate the objective function for a large number of design parameters in the vicinity of the optimum. However, this strategy is ineffective when the objective function is expensive to evaluate. These expensive-to-evaluate functions are common to many disciplines, including tuning the parameters of machine learning algorithms (Bergstra et al., 2011; Snoek et al., 2012), robotics (Tesch et al., 2011; Lizotte et al., 2007), and other engineering design problems (Daniels et al., 2018; Anthony and Keane, 2003; Wiesmann et al., 1998; Volz et al., 2019; Shourangiz-Haghighi et al., 2020). Bayesian optimisation is a principled and efficient technique for the global optimisation of these sorts of functions. The idea underlying Bayesian optimisation is to place a prior distribution over the objective function and then update that prior with observations of the objective function (obtained by expensively evaluating it) in order to yield a posterior predictive distribution. This posterior distribution thus encodes the optimiser’s knowledge of the objective function landscape. It is then used to inform where to make the next observation of the objective function through the use of an acquisition function, which balances the exploitation of regions known to have good performance with the exploration of regions where there is little information about the function’s response. A Gaussian process is a popular choice of prior, because they are intuitive to understand, capable of modelling the objective function accurately with few data, and cheap to evaluate.

The chief contribution of this paper is the introduction of a novel acquisition function for the Bayesian robust optimisation of expensive black-box functions. Although we have phrased robust optimisation as seeking an optimum robust to “small perturbations” of design parameters, the technique we present is applicable to arbitrarily large perturbations. We evaluate the method on 6 benchmark functions and a real-world problem showing that it provides state-of-the-art performance compared with two competing algorithms.

We begin by outlining background material and reviewing similar techniques in Section 2. Section 3 builds upon the previous section by introducing the Bayesian optimisation of the robust domain, and giving a demonstration on a toy function in one dimension. Results of five- and ten-dimensional test problems are presented alongside analysis in Section 4, where we also evaluate the algorithms on two active learning robot pushing problems. Finally, the conclusion and suggestions for future work can be found in Section 5.

2. Background

This section comprises background material in Bayesian optimisation (Section  2.1), Gaussian processes (Section  2.2), and robust optimisation for expensive-to-evaluate functions (Section 2.3).

2.1. Bayesian optimisation

Although stochastic search algorithms, such as evolutionary algorithms, have been popular for the optimisation of black-box functions, Bayesian optimisation is often more attractive, particularly for expensive-to-evaluate functions. Through explicitly modelling the expensive function and accounting for the uncertainty in the model, the search can be guided efficiently to promising areas of the decision space: either those with high certainty of being better than the current best solution, or those with high uncertainty that may be better than the current best. See (Brochu et al., 2010) for an introduction to Bayesian optimisation, and (Shahriari et al., 2016) for a recent comprehensive review.

To be definite and without loss of generality, we assume that the goal of optimisation is to minimise a function $f(\mathbf{x})$, where $\mathbf{x}$ are the design parameters in the feasible space $\mathcal{X}\subset\mathbb{R}^{D}$.

Bayesian optimisation relies on constructing a probabilistic model of $f(\mathbf{x})$. Assume that $f(\mathbf{x})$ has been (expensively) evaluated at $N$ locations $\{\mathbf{x}_{n}\}_{n=1}^{N}$ so that data $\mathcal{D}=\{(\mathbf{x}_{n},f_{n}\triangleq f(\mathbf{x}_{n}))\}_{n=1}^{N}$ are available from which to learn a model. Then Bayesian modelling is used to construct a posterior predictive distribution $p(f\,|\,\mathbf{x},\mathcal{D})$ at any desired location $\mathbf{x}$. Crucially, Bayesian modelling gives not only a prediction of the function value at $\mathbf{x}$, but the posterior distribution quantifies the uncertainty in the prediction as well. Where next to expensively evaluate in $\mathcal{X}$ is determined by an acquisition function, which balances the exploitation of predicted good values of $f(\mathbf{x})$ with exploring uncertain and potentially good regions. Here we use the popular expected improvement (Jones et al., 1998), which has been shown to be effective in practice and for which some theoretical guarantees exist (Bull, 2011). Alternatives such as the probability of improvement (Kushner, 1964) or upper-confidence bound (Srinivas et al., 2010; Brochu et al., 2010) could also be used.

If $f(\mathbf{x})$ is modelled to take the value $\hat{f}(\mathbf{x})$, then the improvement at $\mathbf{x}$ is defined as

[TABLE]

where

[TABLE]

is the best function value from the evaluations thus far. The expected improvement is then

[TABLE]

Gaussian processes are commonly used for modelling $f(\mathbf{x})$ in which case the posterior predictive distribution is itself a Gaussian density (see Section 2.2) with mean $\mu(\mathbf{x})$ and variance $\sigma^{2}(\mathbf{x})$. In this case the expected improvement has the closed analytical form (Jones et al., 1998):

[TABLE]

where $Z=(f^{\star}-\mu(\mathbf{x}))/\sigma(\mathbf{x})$ and $\phi(\cdot)$ and $\Phi(\cdot)$ are the standard Normal density and cumulative distribution functions respectively.

The next (expensive) evaluation is then chosen as that with the greatest expected improvement: $\mathbf{x}^{\prime}=\operatornamewithlimits{argmax}_{\mathbf{x}\in\mathcal{X}}EI(\mathbf{x})$. This location is often discovered by using an evolutionary algorithm to maximise $EI(\mathbf{x})$, which is rapid since $EI(\mathbf{x})$ is computationally cheap to evaluate. The evaluated location and its function value are added to $\mathcal{D}$ and the optimisation proceeds iteratively until some stopping criterion is met or, more commonly, the available computational resources are exhausted.

2.2. Gaussian processes

Gaussian processes (GPs) (Rasmussen and Williams, 2006) are commonly used for Bayesian optimisation due to their flexibility and the simple Gaussian predictive posterior distributions. Briefly, a GP is a collection of random variables, and any finite number of these have a joint Gaussian distribution. Given data $\mathcal{D}$ and a feature vector $\mathbf{x}$, the GP posterior predictive density of the target $\hat{f}$ is Gaussian:

[TABLE]

where the mean and variance of the prediction are given by

[TABLE]

Here $\mathbf{f}=(f_{1},f_{2},\ldots,f_{N})^{T}$ is the vector of evaluated function values at $\mathbf{x}_{1},\mathbf{x}_{2},\ldots,\mathbf{x}_{N}$. Non-linearity in the GP enters through a kernel function $k(\mathbf{x},\mathbf{x}^{\prime}),$ which models the covariance between two feature vectors. The $N\times N$ covariance matrix $\mathbf{K}$ collects these covariances together, $K_{ij}=k(\mathbf{x}_{i},\mathbf{x}_{j})$, and $\mathbf{k}=\mathbf{k}(\mathbf{x})$ is the $N$-dimensional vector of covariances between the training data and $\mathbf{x}$: $k_{n}=k(\mathbf{x},\mathbf{x}_{n})$. There are a number of kernels that could be used, for example the squared exponential function or the Matérn family of covariance functions (Rasmussen and Williams, 2006); here we used the Matérn covariance function with smoothing parameter $\nu=5/2$ which has been recommended for modelling realistic functions (Snoek et al., 2012).

In addition the kernel function depends upon a number of hyper-parameters, $\eta$. Training the GP comprises inferring these hyper-parameters by maximising the marginal likelihood of the data $p(\mathcal{D}\,|\,\eta)$ given by

[TABLE]

Although the log marginal likelihood function landscape may be non-convex and multi-modal, we adopt the standard practice of using a gradient-based optimiser (L-BFGS) with several random starts to estimate good hyper-parameter values (GPy, 2012).

2.3. Robust Optimisation of Expensive Functions

The focus of robust optimisation is to determine function optima in the face of uncertainty. These uncertainties are typically considered to arise in one of three categories: (a) design uncertainties, (b) model uncertainties, or (c) environmental uncertainties (Beyer and Sendhoff, 2007). A variety of methods exist to alleviate the added complexity of searching for robust optima; comprehensive surveys are given by Beyer and Sendhoff (2007) and Gabrel et al. (2014).

In this work we focus on one category of robustness: uncertainties or mis-specification in the design parameters. In this case robust optimisation seeks to find an $\mathbf{x}\in\mathcal{X}$ where some small perturbation $\boldsymbol{\delta}$ does not cause the objective function’s response at $\mathbf{x}+\boldsymbol{\delta}$ to become unacceptably poor. We define the set of all possible perturbations to be

[TABLE]

where $\epsilon$ is the robustness parameter, and $d(\cdot)$ is a function controlling the shape of the robust region, which is often a distance function. The choice of $\epsilon$ and $d(\cdot)$ is determined by the problem owner.

A quality function $Q(\mathbf{x},\Delta_{\epsilon})$ can be used to quantify the robustness of $\mathbf{x}$. There are a number of ways one could formulate $Q$, for example the average performance over the perturbations:

[TABLE]

where $\pi(\mathbf{x}+\boldsymbol{\delta})$ denotes the probability that in practice $f(\cdot)$ will be evaluated at $\mathbf{x}+\boldsymbol{\delta}$ rather than $\mathbf{x}$. Although for analytic convenience $\Delta_{\epsilon}$ is sometimes taken to be unbounded and $\pi(\mathbf{x})$ is taken to be a Gaussian distribution centred on $\mathbf{x}$, in most real applications $\Delta_{\epsilon}$ is taken to be a (small) compact set and $\pi$ is uniform on it.

An alternative quality function, which may be more useful in practice, is to guarantee the worst-case performance (Beyer and Sendhoff, 2007),

[TABLE]

Having defined a quality function, robustness parameter, and shape function the robust optimisation problem may be written as:

[TABLE]

Here we further assume that the objective function is a computationally-expensive-to-evaluate black box. For a general review of robust optimisation with expensive functions see (Chatterjee et al., 2019). When the objective function is expensive to evaluate, optimisers such as evolutionary algorithms (Branke, 1998; Paenke et al., 2006) or particle swarm optimisation (Dippel, 2010) will not be viable due to the large number of function evaluations they demand. Therefore it is essential to apply methods that necessitate only small numbers of observations. In spite of this requirement, relatively few methods exist in the literature to address this problem (Chatterjee et al., 2019).

A related field is that of level set estimation (Gotovos et al., 2013), where the aim is to determine a set of points of equal objective function value, which encompass a region where the performance is guaranteed to be better than some given threshold (either as a specific value or percentage of the unknown optimum). Bogunovic et al. (2016) extended the idea of level set estimation to work in a unified fashion with Bayesian optimisation.

There are a few methods in the literature that use GPs to develop a surrogate model of the expensive function (Jin, 2011; Ong et al., 2006; Lee and Park, 2006), which reduces the computational cost of the optimisation in two ways. Firstly, it enables the surrogate model (rather than the expensive function) to be searched using, for example, an evolutionary algorithm or simulated annealing. Secondly, the use of a surrogate has the clear benefit of curtailing the computational burden of evaluating the robustness of solutions, because the surrogate model can be interrogated instead of the true objective function. Although these methods lighten the computational load, they do not take advantage of the uncertainty in the surrogate model, which is available when using a GP and could be used to help guide the search in subsequent iterations. Picheny et al. (2013) present a review of robust acquisition functions for use with a GP, but these only account for noise in the function’s response.

A state-of-the-art Bayesian approach was presented by ur Rehman et al. (2014). This method exploits a GP with a modified formulation of the expected improvement, which aims to account for the robust performance over a region of the design space. Whilst this technique is shown to be useful for expensive robust optimisation, there are two drawbacks: (a) the uncertainty of the GP is largely disregarded when calculating the modified expected improvement, as only the uncertainty at the estimated worst performing location is considered; and (b) this method is demonstrated with a somewhat substantial number of initial observations (100 in 10 dimensions), which makes it rather unsuitable for very expensive functions. Since this method is considered to be state-of-the-art we have elected to include our own implementation of it for comparison during experimentation.

More recently the StableOpt algorithm has been presented (Bogunovic et al., 2018). This is a confidence-bounds approach that exploits a Gaussian process model to locate a region of good inputs. In essence, at each time-step $n$, a candidate robust location $\mathbf{x}^{\prime}$ is found as

[TABLE]

where $\text{lcb}(\cdot)$ denotes the lower confidence bound acquisition function:

[TABLE]

where $\mu(\mathbf{x})$ and $\sigma(\mathbf{x})$ are the posterior predictive mean and variance at $\mathbf{x}$ (equations (6) and (7)). Adding $\sigma(\mathbf{x})$ to $\mu(\mathbf{x})$ instead of subtracting it obtains the upper confidence bound function, $\text{ucb}(\mathbf{x})$. After identifying a candidate robust location, an expensive evaluation is made at the location $\mathbf{x}_{new}$ such that

[TABLE]

Performing the search in this way ensures that the search operates pessimistically when identifying robust optima, yet optimistically when determining the next location at which to evaluate the expensive function.

3. Bayesian Search for a Robust Optimum

The search for a robust optimum can be distilled to the minimisation of the selected quality function $Q(\mathbf{x},\Delta_{\epsilon})$ (Section 2.3) for a given robust set $\Delta_{\epsilon}$. The obstacle to straightforward optimisation is that evaluating $Q$ for any candidate location is unachievable for continuous domains, because the evaluation of $f(\mathbf{x}+\boldsymbol{\delta})$ is required for all $\boldsymbol{\delta}\in\Delta_{\epsilon}$. In spite of this, one can envisage that a good approximation to $Q$ could be made by aggregating $f(\mathbf{x}+\boldsymbol{\delta})$ over many $\boldsymbol{\delta}\in\Delta_{\epsilon}$. Although such an approach is practicable for cheap-to-evaluate objective functions, the need to optimise expensive functions renders this approximation infeasible as well.

Here we propose to search for the robust optimum by constructing a Gaussian process model of the expensive function $f$, and then using it to estimate the chosen quality function $Q$ to determine the robust quality at a particular location. This allows approximation of the expected improvement for candidate locations by drawing realisations from the Gaussian process and then evaluating the improvement for each realisation based on these estimated quality values. Algorithm 1 shows the main steps in our robust optimisation procedure.

The process is initialised with a small number $N$ of evaluations of $f(\mathbf{x})$, usually chosen using a low discrepancy sampling scheme (Matoušek, 1998), such as a Sobol’ sequence (Sobol’, 1967) or Latin hypercube sampling (Morris and Mitchell, 1995; McKay et al., 1979). These allow a Gaussian process to be constructed (line 2).

In Algorithm 1 we make use of a template $T$ which is a discretised representation of the robust set $\Delta_{\epsilon}$. The exact structure of $T$ is left to the practitioner, but there are several considerations for constructing a template that would be useful in practice. Firstly the set of points constituting the template could be distributed uniformly within the bounds of the robust set, alternatively it would be possible to arrange the set of points to be more densely distributed near the centre of the robust set, thus giving greater weight to perturbations closer to the centre of the robust set; cf equation (10). Finally, for worst case performance guarantees (equation (11)) one might consider placing a larger proportion of points on the boundary of the robust set, because this is where one might assume the extremes of the function response within the robust set would be.

Next, an appropriate quality function (Section 2.3) must be selected to estimate the robust quality of a location. Again, the choice of this function is left to the practitioner, but here we provide two such examples of modified quality functions which operate on sets of response values for each location within the robust set template. Denoting the set of responses estimated from the GP surrogate by $F$, worst-case quality function is estimated as

[TABLE]

and the average quality function by

[TABLE]

Over lines 12 to 18 of Algorithm 1, a realisation $f_{m}$ is drawn from the fitted Gaussian process model for the set of $\mathbf{x}$ that constitute the candidate and best-so-far templates for the robust region ($X_{\delta}$ and $X^{\star}_{\delta}$ respectively). Note that a realisation evaluated at a set of locations $\{\mathbf{x}_{k}\}$ is a draw from a multivariate Gaussian $\mathcal{N}(\mathbf{0},\mathbf{K})$ where $K_{kl}=k(\mathbf{x}_{k},\mathbf{x}_{l})$. Then the improvement for each realisation can be calculated as per line 16 and the expected improvement is the mean of the $M$ realisation-specific improvements (line 19).

Ordinarily, in non-robust Bayesian optimisation, the best returned location is one of the observations $\mathbf{x}_{n}\in X_{N}$, which provides some guarantee about the quality of the returned location. In order to provide a similar guarantee we suggest that it is reasonable to enforce that the returned best robust location should be within the vicinity of at least one expensively evaluated location, that is for $\mathbf{x}^{\star}$ to be considered a valid solution there must exist $\mathbf{x}_{n}\in X_{N}$ such that $d(\mathbf{x}^{\star}-\mathbf{x}_{n})\leq\epsilon$, where $X_{N}$ is the set of $N$ evaluated locations, and $d(\cdot)$ and $\epsilon$ are the distance function and robustness parameter defining the robust set (9). For convenience we write

[TABLE]

for the robust set located at $\mathbf{x}$. Then the requirement that the robust optimum should be near an evaluated location means that the set of possible locations for the best robust optimum found so far is $\mathcal{N}(X_{N})=\cup_{\mathbf{x}_{n}\in X_{N}}\Delta_{\epsilon}(\mathbf{x}_{n})$, which is illustrated in Figure 4. With expensively evaluated locations $X_{N}$, the best robust location is identified as:

[TABLE]

and can be found using, for example, an evolutionary optimiser to search over the $X_{N}$ (Algorithm 1 line 3). We discuss the computation of $\mathbf{x}^{\star}$ in more detail in Section 3.3.

As we show below, although we demand that a robust location is within the vicinity of an evaluated location, it can be advantageous to evaluate $f$ at a location other than that with maximum expected improvement. In Section 3.2 we explore a number of criteria for choosing the location to evaluate. In Algorithm 1 $f$ is expensively evaluated at the location provided by the function $w(\cdot)$ (lines 5 and 7). This sequence is then repeated until convergence is achieved or computational resources are exhausted.

We now illustrate the full procedure with a toy example.

3.1. Toy Example

We illustrate the procedure using the toy one-dimensional function

[TABLE]

for $x\in\mathcal{X}=[0,1]$. For simplicity we restrict the robust set to be an interval, where the robustness parameter $\epsilon=0.1$ defines the span of the interval, and the shape function is given by $d(\delta)=\|\delta\|$, where $x+\delta\in\mathcal{X}$.

Figure 1a shows the toy target function $f$ and the induced robust landscape for the worst-case quality $Q_{max}$ as given by (11). This toy function illustrates how the optimal single point location, namely the minimum of $f(\mathbf{x})$, can exist in a distinct location from the optimal robust region.

The first step is to fit a Gaussian process to an initial set of observations of the expensive function $f$. In this instance we have used an initial set of $N=8$ observations; example draws from the resulting Gaussian process can be seen in Figure 1b.

Each realisation of the Gaussian process can be thought of as a potential $f$ whose quality can be evaluated using the chosen quality function. Figure 1c shows the result of applying the quality function $Q_{max}$ to a drawn realisation $f_{m}$ together with the resulting improvement $I_{m}(x,\Delta_{\epsilon})$ from that realisation.

For simplicity in this example, we have constrained the best uncertainty set to be centred on an observation $x_{n}\in X_{8}$; the best uncertainty set so far (for the depicted realisation) is centred at $x_{m}^{\star}=0.842$, including the observation at $x=0.842$. The robust quality for each realisation is calculated as the difference between the robust quality of the best-so-far location and the candidate location; as shown in Algorithm 1 on line 16 and the robust expected improvement is approximated as an average over all of the realisations, using the procedure given between lines 10 and 19 of Algorithm 1. This is the acquisition function used for determining where to sample $f$ next. Figure 1 (bottom) compares the robust expected improvement with the (usual) single-point expected improvement (4), which clearly demonstrates that the robust expected improvement gives greater weight to searching the more robust regions of design space.

Figure 2 shows the result of continuing the optimisation procedure for three additional iterations; the objective function is evaluated at the location of maximum expected improvement, $\mathbf{x}^{\prime}$. The robust optimiser quickly locates the region of the robust optimum, whereas the single-point optimiser searches the region of the (fragile) global minimum.

3.2. Sampling Location

Non-robust acquisition functions, such as the expected improvement described in Section 2.1, determine a point of maximum acquisition. In contrast, their robust counterparts yield a region, which presents an additional decision in the optimisation process: where within this region should the next observation of the target function be made?

We constrain the location of the new observation $\mathbf{x}_{new}$ to be within the robust set of maximum acquisition $\Delta_{\epsilon}(\mathbf{x}^{\prime})$. Further, we propose the use of a sampling function $w(\mathbf{x}^{\prime})$, which determines where within $\Delta_{\epsilon}(\mathbf{x}^{\prime})$ to locate the next expensive evaluation: $\mathbf{x}_{new}=w(\mathbf{x}^{\prime})$. As we show empirically in section 4 on a range of test functions, the choice of $w$ has a significant impact on the algorithm’s ability to converge to the best uncertainty region. Here we present five suggestions for the sampling function $w$; Figure 3 illustrates each of them.

**Centred observation.: **

Usually uncertainty sets ${\Delta}(\mathbf{x})$ are symmetrical about $\mathbf{x}$ and an obvious choice is to observe the objective function at the location of maximum expected improvement—the centre of the uncertainty set:

[TABLE]

**Most uncertain observation.: **

A maximally explorative approach to improving the estimate of the quality of the predicted best uncertainty set is to observe the expensive function at the location of maximum uncertainty within $\Delta_{\epsilon}(\mathbf{x}^{\prime})$:

[TABLE]

where $\sigma^{2}(\mathbf{x})$ is the predicted variance of the Gaussian process at $\mathbf{x}$, see (7).

**Worst-case prediction.: **

An alternative to improve the estimate of the uncertainty set’s quality is to query at the location of the worst-case predicted value:

[TABLE]

where $\mu(\mathbf{x})$ is the predicted mean of the Gaussian process at $\mathbf{x}$. This strategy has the benefit of confirming or revising the predicted worst case performance with an actual evaluation.

**Uniformly at random.: **

Finally, draw $\mathbf{x}_{new}$ uniformly at random within $\Delta_{\epsilon}(\mathbf{x}^{\prime})$:

[TABLE]

This approach may also be expected to promote exploration, but not in such a directed way as the “most uncertain observation” scheme.

**UCB.: **

Taking inspiration from StableOpt, we locate $\mathbf{x}_{new}$ such that it maximises the upper confidence bound within the proposed robust region:

[TABLE]

where, as with StableOpt, we take $\beta=2$.

3.3. Best-So-Far Robust Location

In standard Bayesian optimisation the best-so-far location $\mathbf{x}^{\star}$ and its function value $f(\mathbf{x}^{\star})$ are simply available because $f$ has been evaluated at all $\mathbf{x}\in X_{N}$, so deciding on the current $\mathbf{x}^{\star}$ is merely a matter of comparing the evaluated locations and then selecting the best one. However, in this robust scheme the improvement for a particular realisation $f_{m}$ requires a procedure to search for the quality of the best robust region: $\max_{\mathbf{x}^{\star}\in\mathcal{N}(X_{N})}Q(\Delta_{\epsilon}(\mathbf{x}^{\star}))$ so that candidate robust centres can be compared with it. Since the evaluation of $Q(\Delta_{\epsilon}(\mathbf{x}))$ requires evaluating the modelled $f$ at a set of locations covering $\Delta_{\epsilon}(\mathbf{x})$, this optimisation in turn requires evaluating the modelled $f$ over all locations that might be covered by the robust region, that is over the extended neighbourhood of $X_{N}$:

[TABLE]

The extended neighbourhood is illustrated in Figure 4.

To avoid this potentially expensive optimisation for every draw of a realisation of $f$, we instead identify the best robust region as:

[TABLE]

where $\hat{Q}$ is evaluated from the mean of the modelled $f$:

[TABLE]

As shown in Algorithm 1, $\mathbf{x}^{\star}$ is determined once each new observation is acquired (line 4).

By evaluating $f_{m}$ at a number of locations $\mathbf{x}_{k}$ in a candidate robust set, the improvement for a particular realisation is then evaluated as

[TABLE]

with the best quality found so far estimated as:

[TABLE]

3.4. Convergence

Essentially we aim to perform Bayesian optimisation on the induced robust landscape $Q(\mathbf{x})$ (11). As a result all of the usual theoretical guarantees for the convergence of Bayesian optimisation already presented in the literature for Bayesian optimisation apply. In particular Bull (2011) shows that Bayesian optimisers using the expected improvement acquisation function as here converge under certain conditions; see also (Vazquez and Bect, 2010).

However, it should be recognised that there are two points where this algorithm makes approximations that affect guarantees of convergence, which we discuss briefly here.

Firstly, the expected improvement is approximated by averaging over sample realisations drawn from the Gaussian process modelling $f(\mathbf{x})$. In the large sample limit this approximation converges to the desired value, like $O(1/M)$ in the number of realisations $M$. However, in practice the acquisition function is approximated with relatively few samples (up to $M=1000$ in the experiments reported here.)

Secondly, we approximate the quality of a robust set $Q(\Delta_{\epsilon}(\mathbf{x}))$ by evaluating the Gaussian process at a set of locations over the template covering the robust set (here we use 60 samples in 2 dimensions, 250 samples in 5 dimensions, and 400 samples in 10 dimensions). Clearly, this approximation may under-estimate the worst case quality and it relies on the fidelity of the Gaussian process approximation. Sufficiently dense sampling of $\Delta_{\epsilon}(\mathbf{x})$ can achieve a good approximation and an alternative is to use a search procedure such as direct (Jones et al., 1993) which can provide an upper bound to the worst-case performance. In practice, however, we have found this to be exorbitantly expensive.

4. Evaluation

We present results of the performance of our method in comparison to the state-of-the-art methods StableOpt (Bogunovic et al., 2018) and the method described by ur Rehman et al. (2014) with $D\in\{5,10\}$ over six common benchmark functions (Mirjalili and Lewis, 2016; Laguna and Martí, 2005; Styblinski and Tang, 1990). We also include some visualisations of results in two dimensions, which serve to demonstrate the differences in search regimes between the compared methods. Figure 5 presents two-dimensional visualisations of each function for reference, which are defined in Table 1.

Each benchmark was selected to test a different aspect of robust optimisation. Function $f_{1}$ (Bumped Bowl) presents a situation where the robust optimum is situated at a local maximum, which tests the ability of an algorithm to overlook the better-performing non-robust region. Our implementation of this benchmark has been modified from (Mirjalili and Lewis, 2016) to ensure that the robust optimum exists exactly at the peak of the local maximum. Benchmarks Levy 03 and Mirjalili & Lewis 04 ($f_{2}$and $f_{4}$) are examples of functions with multiple local minima. In the case of benchmarks $f_{3}$ (Styblinkski-Tang) and $f_{6}$ (Qunitic), the robust optimum resides just outside of the global optimum, which tests robust procedures’ resilience to non-robust regions. Optimisers that exploit the parabolic sphere have difficulty in finding the “step” in $f_{5}$ (Stepped Sphere) containing the optimum, which occupies a vanishingly small proportion of the domain as the number of dimensions increases. This function has been modified from (Mirjalili and Lewis, 2016) to ensure that the size of the step remains significant as $D$ increases: even so, the proportion of $\mathcal{X}$ containing the lower step is only $2^{-D}$.

The robust quality measure used for evaluating candidate locations was the worst-case quality $Q_{max}$ (11), the elected shape function was $d(\boldsymbol{\delta})=\|\boldsymbol{\delta}\|_{2}$, and the robustness parameter $\epsilon=\frac{|u-l|}{8}$, where $u$ and $l$ are the upper and lower bounds of the domain respectively (i.e. $\mathcal{X}=[u,l]^{D}$). See supplementary material S1.1 and S1.2 for comparative results using a “square-shaped” robust region $d(\boldsymbol{\delta})=\|\boldsymbol{\delta}\|_{1}$ in 5 and 10 dimensions respectively.

As noted above the functions were modelled with a Gaussian process with a Matérn 5/2 kernel; kernel parameters were inferred by optimising the marginal likelihood (8) using L-BFGS from 10 random restarts. $M\in\{100,500,1000\}$ realisations were used in the evaluation of the acquisition function, where $M$ increases if no improvement is found at lower values. The location of the acquisition function maximum was found by evaluating it for 1000 Latin hypercube samples and then optimising the most promising 10 of these with L-BFGS-B (Byrd et al., 1995); the optimisation budget was 100 evaluations of the expensive function.

Python code to generate figures and reproduce all experiments is available online111https:github.com/url/completed/on/publication

We evaluated the five sampling schemes proposed in Section 3.2: centred observation (21), most uncertain observation (22), maximised UCB (25), and uniformly at random (24).

To enable paired comparisons, each method was initialised using the same $D+1$ Latin hypercube samples. The experiments were repeated 30 times for statistical comparison. Figures 6 and 7 compare the convergence of each of the methods tested. For the functions $f_{2}$, $f_{3}$, and $f_{6}$, where it is difficult to analytically determine the robust optimum, we have completed 20 repeated trials of CMA-ES (Hansen et al., 2003) in order to locate an approximation for the global robust optimum.

4.1. Analysis

The convergence plots in Figure 6 and Figure 7 show the regret, namely the difference between the state of the optimiser and the value of the true robust minimum. They demonstrate that our robust optimisation procedure is generally capable of locating and exploiting robust optima with a small number of observations of the underlying expensive function. However, we note that all methods perform significantly less well in $D=10$ dimensions and none of the methods are able to locate a good value for the robust optimum for the Levy03 or stepped sphere functions.

A summary of all of the tested algorithms over all of the functions is shown as a critical difference plot (Demšar, 2006) in Figure 8. Of the five competing sampling functions, $w(\cdot)$, sampling at the most uncertain location within the robust region consistently enables better convergence in $D=5$ dimensions. In $D=10$ dimensions the performance of our algorithm with several sampling methods and StableOpt are statistically indistinguishable using the Wilcoxon signed rank test at $p=0.05$.

The success of the most uncertain sampling strategy with $D=5$ is largely because of the increased exploration of the best-so-far robust region, which leads to more even coverage of observations in that key region. The benefits of increased exploration are evident for the Stepped Sphere $f_{5}$ in $D=5$ dimensions where the most uncertain method is the only method of ours that explores the bottom of the parabolic bowl containing the step sufficiently to locate the optimum. In higher dimensions we conjecture that the model of the function is not sufficiently good to allow identification of the most uncertain point.

Figure 9 shows the typical search pattern for each of the five sampling methods, StableOpt, and ur Rehman et al.’s approach after 30 iterations on the two-dimensional Bumped Bowl function (Mirjalili and Lewis, 2016). Each run was initialised from the same set of 3 Latin hypercube samples. It is clear from this example that both the “most uncertain” sampling scheme has made a better estimate of the robust optimum. In addition, the “most uncertain” sampling scheme has been exploratory over the remainder of the domain, and has lead to the most even coverage of observations at the robust optimum.

In all of the benchmark functions, our method has out-performed that of ur Rehman et al.’s competing method, with respect to the value of regret at the final iteration. It should be noted that ur Rehman et al.’s method was able to get closer to the optimum location on numerous occasions, but struggled to accurately settle its robust region, which resulted in a much worse value of regret (due to the fragility of the function landscape).

Generally, and as one might expect, the performance of all of the compared schemes worsens as the number of dimensions increases. Benchmark problem $f_{5}$ is difficult for a GP to model due to the discontinuous downward step in one corner of the domain. This problem is exacerbated in higher dimensions: whilst in two dimensions the downward step covers $\frac{1}{4}$ of the domain, as noted above, the step scales such that it occupies $\frac{1}{2^{D}}$ of the domain. The result is that in ten dimensions the downward step exists in less than one-thousandth of the domain. With the extremely limited number of observations made during our experiments it is unsurprising that none of the sampling schemes, StableOpt nor ur Rehman et al.’s method discovered the step. In five dimensions only the “most uncertain” and StableOpt were able to find the step, which indicates that the improved exploration offered by these schemes presents significant advantages in locating complicated response features.

In general ur Rehman et al.’s approach and our approach using the “centre” sampling scheme do not do well, and in fact perform similarly poorly in similar circumstances. This appears to be a result of both of the approaches tending towards making more exploitative observations, which means that they repeatedly sample in approximately the same location and fail to construct an accurate model of $f$.

Each of the presented methods make significant ground towards improving their quality of robustness, as seen in Figures 6 and 7. And each one generally exhibits a similar convergences curve for the initial 20 iterations as the vicinity of the robust optimum is approached. However, during the remainder of the run it is clear that “most uncertain” and “random” generally outperform the other approaches.

4.2. Real-world application

Here we illustrate our robust optimisation method on a real-world problem. We aim to optimise the input parameters for two active learning robot pushing problems (Wang et al., 2018). In the first problem a single robot hand aims to push an object towards an unknown target location; see Figure 10a. Once the robot has pushed the object, per its input parameters, the push is then evaluated as the distance between the target and the final location of the pushed object. The input parameters to this problem are the starting position of the robot hand, the orientation of the robot hand, and the number of time steps for which the robot hand moves. The direction of travel of the robot hand is constrained to be in the direction of the object’s initial location. We refer to this four-dimensional problem as Push4.

In the second problem, as shown in Figure 10b, two robot hands are tasked to push their respective object towards two unknown targets. This problem poses an additional layer of complexity in that both of the robot hands and objects may block each other’s path if they intersect. The resulting feedback from this problem is the sum of the distance between both of the objects and their respective targets. Due to the possibility of the robots and objects blocking each other, some problem instances do not have perfect solutions. There are eight input parameters to optimise in this case: the starting position of both of the robot hands, the orientation of both of the robot hands, and the number of time steps for which each of the robot hands move. The resulting eight-dimensional problem is referred to as Push8.

As with (Wang and Jegelka, 2017), the object’s initial location in Push4 is always at the centre of the domain, but the target location is changed for each of the optimisation runs. Each of the compared algorithms used the same target locations, so that the runs are directly comparable. Doing this means that we average over instances of the problem class, rather than repeatedly optimising the same function from a different set of starting observations. Likewise, in Push8 the starting position of the objects were fixed, and their respective targets were randomised for each optimisation run. In order to estimate the optimum for each of the different initialisations of the problems we initially sampled the feasible space with 1000 different robot parameters, and then locally optimised the landscape using CMA-ES from the best 20 sampled points.

The results of the experiments, shown in Figure 11, show that on the four-dimensional Push4 that the “most uncertain”, “random”, and “centre” selection methods achieve the best median value of regret. However, in eight dimensions all of the selection methods are capable of improving the function value, but there is little to discern between them. In fact all of the selection methods are considered statistically indistinguishable. This is perhaps unsurprising, because the landscape for the Push8 problem is very complicated: there are many plateaux and discontinuities throughout the landscape, which means that many of the problem instances are very difficult for a Gaussian process to model. It should be noted that in spite of this difficulty, all of our proposed selection methods make some headway to improving the fitness value, and are actually able to find better results than our initial estimates of the optima using random sampling and CMA-ES.

5. Conclusion

This paper has introduced a novel algorithm for the robust optimisation of expensive-to-evaluate functions in the context of Bayesian optimisation. Experiments on a range of commonly-used benchmark functions show that our method is effective at locating robust optima, and able to outperform two state-of-the-art methods from the literature. The algorithm depends upon building a model of the expensive function and then evaluating the improvement with respect to a chosen quality function over realisations drawn from the model. The expectation of these improvements is then used as an acquisition function, which is maximised to inform the next location of the expensive function to evaluate.

Standard (non-robust) Bayesian optimisation, StableOpt and ur Rehman et al.’s method only require the GP to be evaluated once per iteration; we require $M$ evaluations of the GP to form an estimate of the robust improvement. Consequently our methods take longer roughly $M$ times as long to decide on the next location for evaluation of the expensive function. While this additional burden is significant for benchmark functions like these which are trivial to evaluate, the additional time required is insignificant in comparison with the time required to evaluate real-world expensive functions. Note also that interrogating the GP can be made in parallel.

Subject to the demand that the expensive-to-evaluate function be evaluated in the putative robust region, we have demonstrated that the choice of sampling location can markedly affect the convergence rate and quality of the final solution. Methods that promote exploration are generally more effective than exploitative methods; sampling from the location about which the model is most uncertain is effective, although we suspect that in higher dimensions that quality of the surrogate model may not be good enough to properly identify the most uncertain location.

Future work entails the simultaneous optimisation of the location and the shape of the robust region, and improvements in uncertainty estimation in high dimensions.

Acknowledgements.

This work was supported by the Engineering and Physical Sciences Research Council [grant number EP/M017915/1].

S1. Square robust region

S1.1. 5-dimensional results

S1.2. 10-dimensional results