Randomised Gaussian Process Upper Confidence Bound for Bayesian Optimisation

1 Introduction

Global optimisation is a cornerstone of modern scientific innovation. Optimisation of alloys and other materials has allowed us to create massive vehicles that are both strong and light enough to fly. Optimisation in medical science has seen us live longer and healthier lives than previously thought possible. This optimisation usually involves a trial-and-error approach of repeated experiments with different inputs to determine which input produces the most desirable output. Unfortunately, many system are expensive to sample, and the heuristic methods commonly used to select inputs are not sample-efficient. This can lead to these optimisation experiments being prohibitively costly. As such, methods that can select inputs in a sample-efficient manner can lead to faster and cheaper innovation in a wide range of fields.

Bayesian optimisation is one of the most sample efficient methods for optimising expensive, noisy systems [Brochu et al., 2010; Shahriari et al., 2015]. It has shown excellent performance on a range of practical problems, including problems in biomedical science [Turgeon et al., 2016; Gonzalez et al., 2015], materials science [Li et al., 2017; Kikuchi et al., 2018; Ju et al., 2017], and machine learning [Snoek et al., 2012; Klein et al., 2017; Xia et al., 2017]. It does so by using the data from previous samples to generate a statistical model of the system. This model is then used to suggest the next input through an acquisition function. The design of the acquisition function is non-trivial; and is critically important to the algorithms performance. It must balance selecting points with the goal of improving the statistical model (exploration), and selecting points with the goal of utilising the improved statistical model to find the global optima (exploitation). In addition to this, the costly nature of these problems means that it is desirable to have a theoretical guarantee of the algorithm's performance.

∗Contact Author

There are a large range of acquisition functions, all with different balances of exploration and exploitation, but we focus on Gaussian process upper confidence bound (GPUCB) [Srinivas et al., 2010] in this work. This controls its exploration-exploitation trade-off with a single hyperparameter, βt. It has strong theoretical guarantees on the overall convergence rate, but the bound they give is fairly loose. This causes the value of βt to be too large, causing significant over-exploration and hence poor practical performance. In practice, the theoretical guarantees need to be weakened by selecting a far smaller βt.

We introduce a novel modification to the GP-UCB acquisition function that significantly improves its explorationexploitation balance while still having a strong convergence guarantee. This is done by sampling its trade-off parameter with a distribution that allows for a range of exploration factors. However, the distribution is chosen such that convergence is guaranteed to be sub-linear while the sampled βt is generally smaller than the traditional GP-UCB. This reduction leads to a direct improvement on the convergence efficiency in practice. We demonstrate this improved performance over the standard GP-UCB implementation in a range of both synthetic benchmark functions and real-world applications. In summary, our contributions are:

• The development of a modified acquisition function: RGP-UCB.
• A convergence analysis of Bayesian optimisation using RGP-UCB.
• The demonstration of the performance of our method on a range of synthetic and real-world problems.

2 Background

In this section, we provide a brief overview of Bayesian optimisation, with an emphasis on acquisition functions and regret bounds. For a more in-depth overview of Bayesian optimisation, we refer readers to [Rasmussen and Williams, 2006], [Shahriari et al., 2015], and [Brochu et al., 2010].

2.1 Bayesian Optimisation

Bayesian optimisation is a method for optimising expensive, noisy black-box functions. It represents the system being optimised as an unknown function, f. As this is black-box, it is impossible to directly observe. However, we can sample it with an input variable, x, to obtain a noisy output, y = f(x) + ϵnoise, where ϵnoise ∼ N(0, σnoise) is the random noise corrupting the measurement. Bayesian optimisation seeks to efficiently find the optimal input for such systems over a bounded search space, X :

To do this, it creates a statistical model of f using all previously sampled input-output pairs, Dt-1 = {xi, yi}t-1 i=1 . This statistical model is usually a Gaussian process, but other models can be used. The statistical model is then used to create an acquisition function, α(x). This is essentially a map of our belief of how useful a given input will be for optimising the system. As such, it can be used to suggest xt; the next input with which to sample the system. This input and its corresponding output can then be added to the previous data, Dt = Dt-1 ∪ {xt, yt}. This process can be iterated, with the Gaussian process improving with each iteration, until a stopping condition has been met. As Bayesian optimisation is generally used on expensive systems, this stopping condition is often a maximum number of experiments, T .

2.2 Gaussian Process

The Gaussian process is the most common statistical model used in Bayesian optimisation. It models each point in x ∈ X as a Gaussian random variable. As such, it is completely characterised by a mean function, µ(x), and a variance function, σ2 (x). However, in order to model sensibly behaved functions there must be correlation between neighbouring points. This is done by conditioning the distribution on the data via a kernel function. There are many viable kernel functions, but one of the simplest and most popular is the squared exponential kernel as it only depends on a single hyperparameter, the lengthscale l. This is given by kSE(xi, xj) = exp (- ‖xi-xj‖2 2l2 ). Using a kernel such as this, a predictive distribution, f(x) ∼ N(µ(x)σ2 (x)), can be obtained by conditioning the Gaussian process on the data, Dt:

where Kt,(i,j) = k(xi, xj) is the kernel matrix, k∗∗ = kt(x, x), and k∗ = [k(x1, x), k(x2, x), . . . , k(xt, x)]. Here I is the identity matrix with the same dimensions as Kt, and σnoise is the output noise standard deviation.

2.3 Acquisition Functions

Once the Gaussian process has been generated, it must then be used to create an acquisition function. This is chosen such that its global maxima in X will be the best next point to sample:

However, the design of such a function is non-trivial. It must first suggest points spread over the search space to improve the Gaussian process. This is called exploration. Once the Gaussian process has been improved enough, it must then transition to suggesting points in regions that have a high probability of containing the global optima. This is called exploitation. If an acquisition function does not explore enough, it may get stuck exploiting sub-optimal regions and never find the global optima. However, if it explores too much, it may waste costly evaluations improving an already adequate Gaussian process. This makes balancing exploration and exploitation vital. There is a wide range of common acquisition functions, all of which have different balances of exploration and exploitation. These include entropy search (ES) by [Hennig and Schuler, 2012], predictive entropy search (PES) by [Hern´andez-Lobato et al., 2014], knowledge gradient (KG) by [Scott et al., 2011], and others. However, we will only consider the Gaussian process upper confidence bound (GP-UCB) by [Srinivas et al., 2010], Thompson sampling by [Russo and Van Roy, 2014], max-value entropy search (MES) by [Wang and Jegelka, 2017], and expected improvement (EI) by [Jones et al., 1998], with the latter three only being used as baselines.

GP-UCB

GP-UCB [Srinivas et al., 2010] is one of the most intuitive acquisition functions. It balances exploration and exploitation through a single hyperparameter, βt:

Increasing βt makes the acquisition function favour points with high variance, causing more exploration. Decreasing βt will make the acquisition function favour points with high mean, causing more exploitation. However, the selection of βt is not done to optimally balance exploitation and exploration, but is done such that the cumulative regret is bounded. It has been proved that, assuming the chosen kernel satisfies

for some constants a, b > 0, then with probability 1 - δ, the algorithmwill have sub-linear regret if

While the regret bound provided by this choice of βt is desirable, it unfortunately is far larger than needed. This leads to sub-optimal real world performance due to over-exploration. In their own paper, the authors divided the suggested βt by a factor of 5 to achieve better performance [Srinivas et al., 2010].

3 ProposedMethod

We describe our improved GP-UCB acquisition function, randomised Gaussian process upper confidence bound (RGPUCB), and prove that it has a sub-linear regret bound.

Algorithm 1 Bayesian Optimisation with RGP-UCB

Input:Dt0 = {xi, yi}t0 i=1 , Γ scale parameter, θ, # Iterations T, Kernel lengthscale l

• 1: for t = t0 to T do
• 2: Build a Gaussian Process (GP) with Dt.
• 3: Set κt = log( 1 √2π (t2+1)) log(1+θ/2) and draw βt from Γ(κt, θ)
• 4: Optimise the acquisition function, αt(x) = µt(x) + √βtσt(x) to obtain xt+1 and use it to sample yt+1 = f(xt+1) + ϵn from the system.
• 5: Augment the data, Dt+1 =Dy ∪ {xt+1, yt+1}.
• 6: end for
• 7: return x+ = xi such that yi = max t∈[1,T ]yt .

3.1 RGP-UCB

While the standard GP-UCB method has a desirable regret bound, it has relatively poor performance. This is due to the βt used to satisfy this bound being far too large, forcing significant over-exploration. As such, a method for selecting a smaller βt while maintaining a regret bound is desirable. We show that this can be done by sampling βt from a distribution, and as such, we call our method randomised Gaussian process upper confidence bound (RGP-UCB). Doing so means that we bound the Bayesian regret instead of the regret, but it allows for far greater freedom in selecting βt, letting it be set far smaller while still maintaining convergence guarantees. However, we require βt > 0 since a negative βt will punish exploration. We also do not want our distribution to suggest a very large βt as that will cause over-exploration. As such, we draw βt from a Γ(κt, θ) distribution. Examples of this distribution can be seen in Figure 1 and the complete algorithm is shown in Algorithm 1. Much like with standard GP-UCB, the parameters of this distribution are chosen to satisfy a regret bound, as per Theorem 3. However, we show that we only need to set one of the two distribution parameters for the bound to hold. Unlike the standard GP-UCB, this allows us to tune βt to substantially improve our methods performance without compromising its theoretical guarantees.

3.2 Theoretical Analysis

Bayesian optimisation is commonly used in high value problems. As such, theoretical measures of its convergence are desirable. The cumulative regret is one such measure, and is the cornerstone of GP-UCB [Srinivas et al., 2010]. Regret is simply the difference between the current sampled value and the global optima. The cumulative regret is the sum of the regret over all iterations:

As RGP-UCB is probabilistic, we instead need to use the Bayesain regret by Russo et al [Russo and Van Roy, 2014]:

However, their proof was for Gaussian processes with a finite search space, i.e. |X| < ∞. As such, we follow a method similar to [Srinivas et al., 2010] and [Kandasamy et al., 2017] and introduce a discretisation of our search space into an τ d grid of equally spaced points, Xdis. We denote [x]τ as the closest point to x in Xdis.

With this, we can begin bounding the Bayesian regret of our algorithm by decomposing it into components that are easier to bound.

Lemma 1. The Bayesian regret of a probabilistic RGP-UCB algorithm, αt-1, over T iterations can be decomposed as

Proof. This simply follows from the fact that, as xt = argmaxx∈X αt-1(x), we have that α(xt) ≥ α(x∗ ).

With this decomposition, we simply need to find a bound for each term to bound the Bayesian regret. We will start with the second term.

Theorem 1. Assuming that βt is drawn from a Γ(κt, θ) distribution with κt = log( 1 √2π (t2+1)) log(1+θ/2) , the following bound holds

Proof. As the posterior distribution of f(x) at iteration t - 1 is N(µt-1, σ2 t-1 ), the distribution of f(x) - αt-1(x)|βt, is simply N (-√βtσt-1(x),σ2 t-1 (x)). Hence

We note that the exponential term is simply the moment generating function of √βt, which has a closed form for a gamma distribution. This lets us express our inequality as

As we want this to decay at a sub-linear rate, we need

Setting κt-1 = log( 1 √2π (t2+1)) log(1+θt-1/2) to satisfy this, equation R2 is bounded by the following:

Next, we attempt to bound the first component.

Theorem 2. Assuming that βt is drawn from a Γ(k, θ) distribution, the following bound holds

Proof. Using Jensen's Inequality we have

We can then use the Cauchy-Schwartz inequality to get

As βt is a gamma distribution with shape parameter κt and scale parameter θ, its maximum is given by

where γ is the Euler-Mascheroni constant and F -1 (x) is the inverse CDF of β. This finally gives us the following bound:

Finally, we need to bound components R3 and R4. For these, we can use lemma 10 from [Kandasamy et al., 2017].

This means that we have that

With this we can finally find our Bayesian regret bound.

Theorem3. If βt is sampled froma Γ(κt, θ) distributionwith

then the RGP-UCB acquisition function has its Bayesian regret bounded by

where γ is the Euler-Mascheroni constant and F -1 (x) is the inverse CDF of βt

Proof. The result follows simply by combining the bounds for the various components.

4 Results

In this section we present results that demonstrate the performance of RGP-UCB in comparison to other common acquisition functions. We also demonstrate the impact of varying the θ parameter of the gamma distribution used to sample βt. The Python code used for this paper can be found at https://github.com/jmaberk/RGPUCB.

4.1 Experimental Setup

We test our method against a selection of common acquisition functions on a range of Bayesian optimisations problems. These include a range of synthetic benchmark functions and real-world optimisation problems. These are all transformed into continuous maximisation problems for consistency. In each case, the experiment was run for 40d iterations and repeated 10 times with 3d+1 different initial points. The initial points are chosen randomly with a Latin hypercube sample scheme [Jones, 2001]. The methods being tested are:

• Our randomised Gaussian process upper confidence bound with θ = 8 (RGP-UCB θ = 8).
• Our randomised Gaussian process upper confidence bound with θ = 1 (RGP-UCB θ = 1).
• Our randomised Gaussian process upper confidence bound with θ = 0.5 (RGP-UCB θ = 0.5).
• Standard Gaussian process upper confidence bound (GP-UCB) [Srinivas et al., 2010].
• Expected improvement (EI) [Jones et al., 1998].
• Thompson sampling [Russo and Van Roy, 2014].
• Max-value entropy search (MES) [Wang and Jegelka, 2017]

Note that we turn all functions that are traditionally minimised into maximisation problems by taking their negative for consistency. As such, higher results are always better.

4.2 Selection of the Trade-off Parameter

An advantage of our method is that it can change its exploration-exploitation balance without compromising its convergence guarantee. This is done by changing the θ parameter in the βt ∼ Γ(κt, θ) distribution. Increasing θ will increase the expected βt, increasing exploration.

As different problems favour different explorationexploration balances, we tested a range of θ values on a range of different problems. In Figure 2, we show the performance of a range of θ values on an exploitation favouring problem, the Alpine 2 (5D) function, and an exploration favouring problem, the Dropwave (2D) function.

It was found that θ = 8 gives good performance on both the above exploration-favouring problem and other similar problems tested. Likewise, θ = 0.5 is a good choice for exploitation favouring problems. We also note that θ = 1 has decent performance on both problems, making it a good choice for problems where the required exploration-exploitation balance is completely unknown.

4.3 Synthetic Benchmark Functions

The first demonstration of our methods performance is the optimisation of several common synthetic benchmark functions. These are the Dropwave (2D), Sphere (4D), Alpine 2 (5D), and Ackley (5D) functions 1 . Results for these are shown in Figure 3.

Here we can see that RGP-UCB has competitive performance in all of the above cases. In general, it does significantly better than the standard GP-UCB and the Thompson sampling acquisition functions. EI has better early performance in many cases, as it starts exploiting earlier. However, RGP-UCB tends to have better exploration and therefore often able to beat it in the long-term. MES, on the other hand, has poor early performance in general but catches up to other methods later. However, explorative RGP-UCB (θ = 8) beat it in the exploration favouring dropwave function. In the other synthetic functions, both the balanced (θ = 1) and exploitative (θ = 0.5) RGP-UCB achieved superior results.

The leftmost functions were chosen to be pathological cases which disproportionately favours exploration (Dropwave 2D) and exploitation (Sphere 4D). These can be seen as best-case examples for GP-UCB and EI respectively. However, RGP-UCB is able to out-perform them even on these if its θ parameter is chosen properly, and does so while maintaining its convergence guarantee. This formulation of EI does not have a known regret bound as it follows the standard implementation and hence doesn't satisfy the assumptions required by the current bounds [Bull, 2011; Ryzhov, 2016; Wang and de Freitas, 2014; Nguyen et al., 2017].

1 All benchmark functions use the recommended parameters from https://www.sfu.ca/∼ssurjano/optimization.html

In the middle two plots, RGP-UCB is superior even with the conservative parametrisation of θ = 1.

4.4 Machine Learning Hyperparameter Tuning

Our first demonstration of the real-world performance of RGP-UCB is the hyperparameter tuning of a support vector regression (SVR) [Drucker et al., 1997] algorithm. This is the support vector machine classification algorithm extended to work on regression problems, with performance measured in root mean squared error (RMSE). It has three hyperparameters, the threshold, ϵ, the kernel parameter, γ, and a soft margin parameter, C. All experiments are done with the public Space GA scale dataset 2 . The results are shown in Figure 3.

We can see that, while MES is superior here, the final performance of all three variants of our method exceeds that of standard GP-UCB. The high exploitation and balanced variants are competitive with EI, with the former achieving higher final performance. As with many real-world problems, SVR is known to favour higher exploitation, and is therefore an example of when the user would know to try a smaller θ.

4.5 Materials Science Application: Alloy Heat Treatment

The second demonstration of RGP-UCB's performance is the optimisation of a Aluminium-Scandium alloy heat treatment simulation [Robson et al., 2003]. The goal of the simulation is to optimise the resulting alloys hardness, measured in MPa. The hardening process is controlled through multiple cooking stages, each with two hyperparametrs, the duration and a temperature. As we use a two-stage cooking simulation, there is a total of four hyperparameters to optimise through Bayesian optimisation. The results are shown in Figure 3.

2 Dataset can be found at https://www.csie.ntu.edu.tw∼cjlin/ libsvmtools/datasets/regression.html

Outside of MES, the results are very similar to the previous SVR example, with the high-exploitation method having the best performance and the balance method being competitive with EI.

4.6 Conclusion

We have developed a modified UCB based acquisition function that has substantially improved performance while maintaining a sub-linear regret bound. We have proved that this bound holds in terms of Bayesian regret while allowing for some flexibility in the selection of its parameters. We have also demonstrated the impact of said parameters on the performance. Moreover, we have shown that its performance is competitive or greater than existing methods in a range of synthetic and real-world applications. In particular, we show that it out-performs the standard GP-UCB implementation on most problems, even when its hyperparamter is chosen conservatively.

Acknowledgements

This research was supported by an Australian Government Research Training Program (RTP) Scholarship awarded to JMA Berk, and was partially funded by the Australian Government through the Australian Research Council (ARC). Prof Venkatesh is the recipient of an ARC Australian Laureate Fellowship (FL170100006).