Batch Bayesian optimisation via density-ratio estimation with guarantees

Bayesian optimisation (BO) algorithms approach the problem in Equation 1 via sequential decision making [1]. At each iteration, BO selects a query point by maximising an acquisition function $a$:

The acquisition function encodes information provided by the observations collected so far ${\mathcal{{D}}}_{t}:=\{{\boldsymbol{\mathbf{x}}}_{i},y_{i}\}_{i=1}^{t-1}$ using a probabilistic model over $f$, typically a Gaussian process (GP) [2], conditioned on the data. After collecting an observation $y_{t}$, the dataset is updated with the new query-observation pair ${\mathcal{{D}}}_{t}:={\mathcal{{D}}}_{t-1}\cup{{\boldsymbol{\mathbf{x}}}_{t},y_{t}}$. This process then repeats for a given number of iterations $T$.

The expected improvement (EI) [4, 14] is a popular acquisition function in the BO literature and the basis for many BO algorithms. At each iteration $t\geq 1$, one can define $\tau:=\min_{i<t}y_{t}$ as an incumbent target. EI is then defined as:

In the case of a GP prior on $f|{\mathcal{{D}}}_{t-1}\sim\mathcal{GP}(\mu_{t-1},k_{t-1})$, the EI is available in closed form as a function of the GP posterior. However, the EI may be reformulated without the need for a prior.

Under mild assumptions, Bergstra et al. [15] showed that the EI can be formulated as a density ratio between two probability distributions. Let $\ell({\boldsymbol{\mathbf{x}}}):=p({\boldsymbol{\mathbf{x}}}|y\leq\tau)$ represent the probability density over ${\boldsymbol{\mathbf{x}}}\in{\mathcal{{X}}}$ conditioned on the observation $y$ being below a threshold $\tau\in\mathbb{{R}}$. Conversely, let $g({\boldsymbol{\mathbf{x}}}):=p({\boldsymbol{\mathbf{x}}}|y>\tau)$. For $\gamma\in[0,1]$, the $\gamma$-relative density ratio between these two densities is:

noting that $\gamma=0$ leads to the ordinary probability density ratio definition, ${}_{0}({\boldsymbol{\mathbf{x}}})=\nicefrac{{\ell({\boldsymbol{\mathbf{x}}})}}{{g({\boldsymbol{\mathbf{x}}})}}$. Now if we choose $\tau:=\Phi^{-1}(\gamma)$, where $\Phi(s):=p(y\leq s)$ represents the cumulative distribution function of the marginal distribution of observations,¹¹1Note that $p(y\leq s)=\int_{{\mathcal{{X}}}}p(y\leq s|{\boldsymbol{\mathbf{x}}})p({\boldsymbol{\mathbf{x}}}){\mathop{}\operatorname{d}}{\boldsymbol{\mathbf{x}}}$, where we may assume $p({\boldsymbol{\mathbf{x}}})$ uniform. for $s\in\mathbb{{R}}$, and then replace $\tau$ in Equation 3, Bergstra et al. [15] have shown that²²2Bergstra et al. [15] and Tiao et al. [9] also rely on the mild assumption that $p({\boldsymbol{\mathbf{x}}}|y)\approx\ell({\boldsymbol{\mathbf{x}}})$ for all $y\leq\tau$. $a_{\operatorname{EI}}({\boldsymbol{\mathbf{x}}})\propto{}_{\gamma}({\boldsymbol{\mathbf{x}}})$, for ${\boldsymbol{\mathbf{x}}}\in{\mathcal{{X}}}$. Based on this fact, Tiao et al. [9] showed:

where $\pi({\boldsymbol{\mathbf{x}}}):=p(y\leq\tau|{\boldsymbol{\mathbf{x}}})$ can be approximated by a probabilistic classifier trained with a proper scoring rule, such as the binary cross-entropy loss:

Other examples of proper scoring rules include the least-squares loss, which leads to probabilistic least-squares classifiers [16], and the zero-one loss. We refer the reader to Gneiting and Raftery [17] for a review and theoretical analysis on this topic.

BORE is summarised in Algorithm 1. As seen, the marginal observations distribution CDF $\Phi(s):=p(y\leq s)$ is replaced by the empirical approximation $\hat{\Phi}_{t}(s):=\frac{1}{t}\sum_{i=1}^{t}\mathbb{I}[y_{i}\leq s]$ and its corresponding quantile function $\hat{\Phi}_{t}^{-1}$. At each iteration, observations are labelled according to the estimated $\gamma$th quantile $\tau$, and a classifier $\hat{\pi}_{t}$ is trained by minimising the loss $\mathcal{L}[\pi|\tilde{\mathcal{{D}}}_{t}]$ over the data points $\tilde{\mathcal{{D}}}_{t}$. A query point ${\boldsymbol{\mathbf{x}}}_{t}$ is chosen by maximising the classifier’s probabilities, which in our case corresponds to maximising the expected improvement. A new observation is collected, and the algorithm continues running up to a given number of iterations $T$. As demonstrated, no explicit probabilistic model for $f$ is needed, only a classifier, which can be efficiently trained via, e.g., stochastic gradient descent.

We consider the case of probabilistic least-squares (PLS) classifiers [20, 21]. In particular, we model a probabilistic classifier $\pi:{\mathcal{{X}}}\to[0,1]$ as an element of a reproducing kernel Hilbert space (RKHS) ${\mathcal{{H}}}$ associated with a positive-definite kernel $k:{\mathcal{{X}}}\times{\mathcal{{X}}}\to\mathbb{{R}}$. A RKHS is a space of functions equipped with inner product $\langle\cdot,\cdot\rangle_{k}$ and norm $\lVert\cdot\rVert_{k}:=\sqrt{\langle\cdot,\cdot\rangle_{k}}$ [22]. For the purposes of this analysis, we will also assume that $k({\boldsymbol{\mathbf{x}}},{\boldsymbol{\mathbf{x}}})\leq 1$, for all ${\mathcal{{X}}}$.⁴⁴4This assumption can always be satisfied by proper scaling. This setting allows for both linear and non-parametric models. Gaussian assumptions on the function space would lead us to GP-based PLS classifiers [2], but we are not restricted by Gaussianity in our analysis. If the kernel $k$ is universal, as $\Phi_{\epsilon}$ is injective, we can also see that the RKHS assumption allows for modelling any continuous function.

For a given $\tau\in\mathbb{{R}}$, a PLS classifier is obtained by minimising the regularised squared-error loss:

where $\lambda>0$ is a given regularisation factor and $z_{i}:=\mathbb{I}[y_{i}\leq\tau]\in\{0,1\}$. In the RKHS case, the solution to the problem above is available in closed form [23, 16] as:

where ${\boldsymbol{\mathbf{k}}}_{t}({\boldsymbol{\mathbf{x}}}):=[k({\boldsymbol{\mathbf{x}}},{\boldsymbol{\mathbf{x}}}_{1}),\dots,k({\boldsymbol{\mathbf{x}}},{\boldsymbol{\mathbf{x}}}_{t})]^{\mathsf{T}}\in\mathbb{{R}}^{t}$, ${\boldsymbol{\mathbf{{K}}}}_{t}:=[k({\boldsymbol{\mathbf{x}}}_{i},{\boldsymbol{\mathbf{x}}}_{j})]_{i,j=1}^{t}\in\mathbb{{R}}^{t\times t}$ and ${\boldsymbol{\mathbf{z}}}_{t}:=[z_{1},\dots,z_{t}]^{\mathsf{T}}\in\mathbb{{R}}^{t}$. This PLS approximation may not yield a valid classifier, since it is possible that $\hat{\pi}_{t}({\boldsymbol{\mathbf{x}}})\notin[0,1]$ for some ${\boldsymbol{\mathbf{x}}}\in{\mathcal{{X}}}$. However, it allows us to place a confidence interval on the optimal classifier’s prediction, as presented in the following theorem, which is based on theoretical results from the online learning literature [24, 25]. Our proofs can be found in the supplement.

We now consider BORE with a PLS classifier. For this analysis, we will assume an ideal setting where $\tau$ is fixed, possibly corresponding to the true $\gamma$th quantile of the observations distribution. However, our results hold for any choice of $\tau\in\mathbb{{R}}$ and can therefore be assumed to approximately hold for a varying $\tau$ which is converging to a fixed value. In this setting, the algorithm’s choices are: given by:

where $\hat{\pi}_{t}$ is the estimator in Equation 13. we can then apply Theorem 1 to the classifier-based regret in Equation 11 to obtain a regret bound. For this result, we will also need the following quantity:

where the maximisation is taken over the discrete set of locations $\{{\boldsymbol{\mathbf{x}}}_{i}\}_{i=1}^{N}\subset{\mathcal{{X}}}$ and ${\boldsymbol{\mathbf{{K}}}}_{N}:=[k({\boldsymbol{\mathbf{x}}}_{i},{\boldsymbol{\mathbf{x}}}_{j})]_{i,j=1}^{N}$. This quantity denotes the maximum information gain of a Gaussian process model after $N$ observations. We are now ready to state our theoretical result regarding BORE’s regret.

As Theorem 2 shows, the regret of the BORE algorithm in the PLS setting is comprised of two components. The first term is related to the regret of a GP-UCB algorithm [see 26, Thr. 3] and its known to grow sub-linearly for a few popular kernels, such as the squared exponential and the Matérn class [3, 27]. The second term, however, reflects the uncertainty of the algorithm around the optimum location ${\boldsymbol{\mathbf{x}}}^{*}$. If the algorithm never samples at that location, this second summation might have a mostly linear growth, which will not lead to a vanishing regret. In fact, if we consider Equation 13 and a RKHS with a translation-invariant kernel, we see that, as soon as an observation $z_{t}=1$ is collected at a location ${\boldsymbol{\mathbf{x}}}_{t}\neq{\boldsymbol{\mathbf{x}}}^{*}$, that location will constitute the maximum of the classifier output. Then the algorithm would keep returning to that same location, missing opportunities to sample at ${\boldsymbol{\mathbf{x}}}^{*}$.

It is worth noting that Theorem 2 reflects the regret of BORE in an idealistic setting where the algorithm uses the optimal PLS estimator in the function space ${\mathcal{{H}}}$. However, if we train a parametric classifier, such as a neural network, via gradient descent, the behaviour will not necessarily be the same, and the algorithm might still achieve a good performance. In the original BORE paper, for instance, a parametric classifier is trained by minimising the binary cross-entropy loss [9] and leads to a successful performance in experiments. Neural network models trained via stochastic gradient descent are known to provide approximate samples of a posterior distribution [28, 29], instead of an optimal best-fit predictor, which might make BORE behave like Thompson sampling [30] (see discussion in the appendix). Nevertheless, Theorem 2 still shows us that BORE may get stuck into local optima, which is not ideal for BO methods. In the next section, we present an extension of the BORE framework which addresses this shortcoming.

We propose replacing $\hat{\pi}_{t}$ in Algorithm 1 by an upper confidence bound which is such that:

which with probability greater than $1-\delta$, given $\delta\in(0,1)$. Therefore, $\pi_{t,\delta}({\boldsymbol{\mathbf{x}}})$ represents an upper quantile over the optimal class probability $\pi^{*}({\boldsymbol{\mathbf{x}}})$. BORE++ selects ${\boldsymbol{\mathbf{x}}}_{t}\in\operatornamewithlimits{argmax}_{{\boldsymbol{\mathbf{x}}}\in{\mathcal{{X}}}}\pi_{t-1,\delta}({\boldsymbol{\mathbf{x}}})$.

To derive an upper confidence bound on a classifier’s predictions $\pi({\boldsymbol{\mathbf{x}}})$, we can take a few different approaches. For a parametric model $\pi_{\boldsymbol{\mathbf{\theta}}}$, a Bayesian model updating the posterior $p({\boldsymbol{\mathbf{\theta}}}|{\mathcal{{D}}}_{t})$ leads to a corresponding predictive distribution over $\pi_{\boldsymbol{\mathbf{\theta}}}({\boldsymbol{\mathbf{x}}})$. This is the case of ensemble models [31], for instance, where we approximate predictions $p(y\leq\tau_{t}|{\boldsymbol{\mathbf{x}}},{\mathcal{{D}}}_{t})\approx\frac{1}{M}\sum_{i=1}^{M}\pi_{{\boldsymbol{\mathbf{\theta}}}^{i}}({\boldsymbol{\mathbf{x}}})$ with ${\boldsymbol{\mathbf{\theta}}}^{i}\sim p({\boldsymbol{\mathbf{\theta}}}|{\mathcal{{D}}}_{t})$. Instead of using the expected class probability, however, BORE++ uses an (empirical) quantile approximation for $\pi_{t,\delta}$ to ensure Equation 19 holds. Bayesian neural networks [32], random forests [33], dropout methods, etc. [34], also constitute valid approaches for predictive uncertainty estimation. An alternative approach is to place a non-parametric prior over $\pi^{*}$, such as a Gaussian process model [2], which allows for the modelling of uncertainty directly in the function space where $\pi^{*}$ lies. In the next section, we present a concrete derivation of BORE++ for the PLS classifier setting which takes the non-parametric perspective and allows us to derive theoretical performance guarantees.

In the PLS setting, the result in Theorem 1 gives us a closed-form expression for a classifier upper confidence bound satisfying the condition in Equation 19. Given $\delta\in(0,1)$, we set:

where $\sigma_{t}$ and $\beta_{t}$ are set according to Theorem 1. We then obtain the following result for BORE++.

According to Theorem 3, the regret of BORE++ vanishes if the maximum information gain $\xi_{T}$ grows sub-linearly, since $\lim_{T\to\infty}\frac{{{R}}_{T}}{T}=0$ and $\min_{t\leq T}r_{t}\leq\frac{{{R}}_{T}}{T}$. Sub-linear growth is known to be achieved for popular kernels, such as the squared exponential, the Matérn family and linear kernels [3, 27]. This result also tells us that theoretically BORE++ performs no worse than GP-UCB since they share similar regret bounds [3, 26]. However, in practice, the BORE++ framework offers a series of practical advantages over GP-UCB, such as no need for an explicit surrogate model, and a classifier which does not need to be a GP and can therefore be more flexible and scalable to high-dimensional problems and large amounts of data. The connection with GP-UCB, instead, brings us new insights into how the density-ratio BO algorithm can still share some of the well known guarantees of traditional BO methods.

Applying an optimisation-as-inference perspective to BORE, we can formulate a batch BO algorithm which does not require an explicit regression model for $f$. The classifier $\hat{\pi}({\boldsymbol{\mathbf{x}}})\approx p(y\leq\tau|{\boldsymbol{\mathbf{x}}})$ naturally turns out as a likelihood function over query locations ${\boldsymbol{\mathbf{x}}}\in{\mathcal{{X}}}$. Since the search space ${\mathcal{{X}}}$ is compact, we can assume a uniform prior distribution $p({\boldsymbol{\mathbf{x}}})\propto 1$. Also note that the normalisation constant in this case is simply $\int_{\mathcal{{X}}}p(y\leq\tau|{\boldsymbol{\mathbf{x}}})p({\boldsymbol{\mathbf{x}}}){\mathop{}\operatorname{d}}{\boldsymbol{\mathbf{x}}}=p(y\leq\tau)=\gamma$. Our posterior distribution then becomes:

Therefore, we formulate a batch version of BORE as an inference problem aiming for:

where $D_{\mathrm{KL}}(q||\ell)$ denotes the Kullback-Leibler (KL) divergence between $q$ and $\ell$, and ${\mathcal{{P}}}$ represents the space of probability distributions over ${\mathcal{{X}}}$. Sampling from $\ell$ would allow us to obtain the points of interest in the search space, including the optimum ${\boldsymbol{\mathbf{x}}}^{*}$ and other locations where $y\leq\tau$. However, as the true $p(y\leq\tau|{\boldsymbol{\mathbf{x}}})$ is unknown, we instead formulate a proxy inference problem with respect to a surrogate target distribution $\hat{p}_{t}$ based on the classifier model. For BORE, we set $\hat{p}_{t}({\boldsymbol{\mathbf{x}}})\propto\hat{\pi}_{t-1}({\boldsymbol{\mathbf{x}}})$, while for BORE++ the setting is $\hat{p}_{t}({\boldsymbol{\mathbf{x}}})\propto\pi_{t-1,\delta}({\boldsymbol{\mathbf{x}}})$. In contrast to $\ell({\boldsymbol{\mathbf{x}}})\propto p(y\leq\tau|{\boldsymbol{\mathbf{x}}})$, the normalisation constant for the surrogate distributions is unknown, leading us to a proxy problem of minimising $q_{t}\in\operatornamewithlimits{argmin}_{q\in{\mathcal{{P}}}}D_{\mathrm{KL}}(q||\hat{p}_{t})$ at each iteration $t\geq 1$. This variational inference problem above can be efficiently solved via Stein variational gradient descent (SVGD) [35], described next.

In our implementation, we apply SVGD to approximately sample a batch ${\mathcal{{B}}}_{t}:=\{{\boldsymbol{\mathbf{x}}}_{t,i}\}_{i=1}^{M}$ of $M\geq 1$ points from $\hat{p}_{t}$. Other approximate inference algorithms could also be applied. One of the main advantages of SVGD, however, is that it encourages diversification in the batch, capturing the possible multimodality of $\hat{p}_{t}$. Given the batch locations, observations can be collected in parallel and then added to the dataset to update the classifier model.

SVGD is an approximate inference algorithm which represents a variational distribution $q$ as a set of particles $\{{\boldsymbol{\mathbf{x}}}^{i}\}_{i=1}^{M}$ [35]. The particles are initialised as i.i.d. samples from an arbitrary base distribution and then optimised via a sequence of smooth transformations towards the target distribution, which in our case corresponds to $\hat{p}_{t}\propto\pi_{t-1,\delta}$. The SVGD steps are given by:

where $k:{\mathcal{{X}}}\times{\mathcal{{X}}}\to\mathbb{{R}}$ is a positive-definite kernel, and $\alpha>0$ is a small step size. Intuitively, the first term in the definition of ${\boldsymbol{\mathbf{\zeta}}}_{t}$ guides the particles to the modes of $\hat{p}_{t}$, while the second term encourages diversification by repelling nearby particles. Theoretical convergence guarantees [36, 37] and practical extensions, such as second-order methods [38, 39] and derivative-free approaches [40], have been proposed in the literature. Further details on SVGD can be found in Liu and Wang [35].

We follow the derivations in Oliveira et al. [41] to derive a distributional regret bound for batch BORE++ with respect to its target sampling distribution $\ell$, which is presented in the following result.

As in the case of non-batch BORE++, the distributional regret bounds for the batch algorithm also grow sub-linearly for most popular kernels, leading to an asymptotically vanishing simple regret. Although different, to compare the distributional regret of batch BORE++ with the non-distributional regret bounds for BORE++, we may consider a case where $\tau$ is set to the function minimum $\tau:=f({\boldsymbol{\mathbf{x}}}^{*})=\min_{{\boldsymbol{\mathbf{x}}}\in{\mathcal{{X}}}}f({\boldsymbol{\mathbf{x}}})$ and the observation noise is small. In this case, the batch sampling distribution would converge to a Dirac at the optimum, so that $\mathbb{E}_{{\boldsymbol{\mathbf{x}}}\sim\ell}[f({\boldsymbol{\mathbf{x}}})]\approx f({\boldsymbol{\mathbf{x}}}^{*})$. Compared to the regret of non-batch BORE++ (Theorem 3) after collecting an equivalent number of observations $T^{\prime}:=MT$, the expected regret of the batch version of BORE++ after $T$ iterations is then lower by a factor of $\xi_{T}/\xi_{MT}$, noting that $\xi_{T}\leq\xi_{T^{\prime}}=\xi_{MT}$. Therefore, batch BORE++ should be able to achieve lower regret per runtime than sequential BORE++ with an equivalent number of observations.