Randomised Gaussian Process Upper Confidence Bound for 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)+\epsilon_{noise}$, where $\epsilon_{noise}\sim N(0,\sigma_{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, $\mathcal{X}$:

To do this, it creates a statistical model of $f$ using all previously sampled input-output pairs, $D_{t-1}=\{x_{i},y_{i}\}_{i=1}^{t-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, $\alpha(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 $x_{t}$; the next input with which to sample the system. This input and its corresponding output can then be added to the previous data, $D_{t}=D_{t-1}\cup\{x_{t},y_{t}\}$. 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$.

The Gaussian process is the most common statistical model used in Bayesian optimisation. It models each point in $x\in\mathcal{X}$ as a Gaussian random variable. As such, it is completely characterised by a mean function, $\mu(x)$, and a variance function, $\sigma^{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

Using a kernel such as this, a predictive distribution, $f(x)\sim\mathcal{N}(\mu(x)\sigma^{2}(x))$, can be obtained by conditioning the Gaussian process on the data, $D_{t}$:

where $\mathbf{K}_{t,(i,j)}=k(x_{i},x_{j})$ is the kernel matrix, $k_{*\!*}=k_{t}(x,x)$, and $\mathbf{k}_{*}=[k(x_{1},x),k(x_{2},x),\ldots,k(x_{t},x)]$. Here $\mathbf{I}$ is the identity matrix with the same dimensions as $\mathbf{K}_{t}$, and $\sigma_{noise}$ is the output noise standard deviation.

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 $\mathcal{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ández-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), and expected improvement (EI) by Jones et al. (1998), with the latter two only being used as baselines.

While the standard GP-UCB method has a desirable regret bound, it has relatively poor performance. This is due to the $\beta_{t}$ used to satisfy this bound being far too large, forcing significant over-exploration. As such, a method for selecting a smaller $\beta_{t}$ while maintaining a regret bound is desirable. We show that this can be done by sampling $\beta_{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 $\beta_{t}$, letting it be set far smaller while still maintaining convergence guarantees. However, we require $\beta_{t}>0$ since a negative $\beta_{t}$ will punish exploration. We also do not want our distribution to suggest a very large $\beta_{t}$ as that will cause over-exploration. As such, we draw $\beta_{t}$ from a $\Gamma(\kappa_{t},\theta)$ 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 $\beta_{t}$ to substantially improve our methods performance without compromising its theoretical guarantees.

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. $|\mathcal{X}|<\infty$. 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 $\tau^{d}$ grid of equally spaced points, $\mathcal{X}_{\mathsf{dis}}$. We denote $[x]_{\tau}$ as the closest point to $x$ in $\mathcal{X}_{\mathsf{dis}}$.

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

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.

Next, we attempt to bound the first component.

Finally, we need to bound components $R_{3}$ and $R_{4}$. 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.

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. 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 $\theta=8$ (RGP-UCB $\theta=8$).

• •

  Our randomised Gaussian process upper confidence bound with $\theta=1$ (RGP-UCB $\theta=1$).

• •

  Our randomised Gaussian process upper confidence bound with $\theta=0.5$ (RGP-UCB $\theta=0.5$).

• •

  Standard Gaussian process upper confidence bound (GP-UCB) Srinivas et al. (2010).

• •

  Expected improvement (EI) Jones et al. (1998).

• •

  Thompson sampling (Thompson) Russo and Van Roy (2014).

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.

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 $\theta$ parameter in the $\beta_{t}\sim\Gamma(\kappa_{t},\theta)$ distribution. Increasing $\theta$ will increase the expected $\beta_{t}$, increasing exploration.

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

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

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²²2All benchmark functions use the recommended parameters from https://www.sfu.ca/~ssurjano/optimization.html. 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.

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 $\theta$ 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).

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

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, $\epsilon$, the kernel parameter, $\gamma$, and a soft margin parameter, $C$. All experiments are done with the public Space GA scale dataset ³³3 Dataset can be found at https://www.csie.ntu.edu.tw~cjlin/libsvmtools/datasets/regression.html. The results are shown in Figure 3.

We can see that 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 $\theta$.

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.

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.

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.