Preferential Bayesian optimisation with Skew Gaussian Processes

Bayesian optimization (BO) is a powerful tool for global optimisation of expensive-to-evaluate black-box objective functions (Brochu et al.,, 2010; Mockus,, 2012). However, in many realistic scenarios, the objective function to be optimized cannot be easily quantified. This happens for instance in optimizing chemical and manufacturing processes, in cases where judging the quality of the final product can be a difficult and costly task, or simply in situations where only human preferences are available, like in A/B tests (Siroker and Koomen,, 2013). In such situations, Preferential Bayesian optimization (PBO) (González et al.,, 2017) or more general algorithms for active preference learning should be adopted (Brochu et al.,, 2008; Zoghi et al.,, 2014; Sadigh et al.,, 2017; Bemporad and Piga,, 2019). These approaches require the users to simply compare the final outcomes of two different experiments and indicate which they prefer. Indeed, it is well known that humans are better at comparing two options rather than assessing the value of “goodness” of an option (Thurstone,, 1927; Chernev et al.,, 2015).

This contribution is focused on PBO. As in state-of-the-art approach for PBO (González et al.,, 2017), we use a Gaussian process (GP) as a prior distribution of the latent preference function and a probit likelihood to model the observed pairwise comparisons. However, our contribution differs from and improves González et al., (2017) in several directions.

First, the state-of-the-art PBO methods usually approximate the posterior distribution of the preference function via Laplace’s approach. On the other hand, we compute the exact posterior distribution, which we prove to be a Skew Gaussian Process (SkewGP) (recently introduced in Benavoli et al., (2020) for binary classification).¹¹1In particular, the present works extends the results in Benavoli et al., (2020) by showing that SkewGPs are conjugate to probit affine likelihoods. This allows us to apply SkewGPs for preference learning. Through several examples, we show that the posterior has a strong skewness, and thus any approximation of the posterior that relies on a symmetric distribution (such as Laplace’s approximation) results in sub-optimal predictive performances and, thus, slower convergence in PBO.

Second, we propose computationally efficient methods to draw samples from the posterior that are then used to calculate the acquisition function.

Third, we extend standard acquisition functions used in BO to deal with preference observations and propose a new acquisition function for PBO obtained by combining the dueling information gain with the expected probability of improvement.

Fourth, we define an affine probit likelihood to model the observations. Such a likelihood allows us to handle, in a unified framework, mixed categorical and preference observations. These two different types of information are usually available in manufacturing, where some parameters may cause the process to fail and, therefore, producing no output. In standard BO, where the function evaluation is a scalar, a common way to address this problem is by penalizing the objective function when no output is produced, but this approach is not suitable in PBO as the output is not a scalar.

The rest of the paper is organized as follows. A review on skew-normal distributions and SkewGP is provided in Section 2. The main results of the paper are reported in Section 3, where we show that the posterior distribution of the latent preference function is a SkewGP under the proposed affine probit likelihood. The marginal likelihood is derived and maximized to chose the model’s hyper-parameters. An illustrative example is presented in 4 to show the drawbacks of Laplace’s approximation, thus highlighting the benefits of computing the exact SkewGP posterior distribution. PBO through SkewGP is discussed in Section 5, where extensive tests with different acquisition functions are reported and clearly show that PBO based on SkewGP consistently outperforms Laplace’s approximation both in terms of convergence speed and computational time.

In this section we provide details on the skew-normal distribution. The skew-normal (O’Hagan and Leonard,, 1976; Azzalini,, 2013) is a large class of probability distributions that generalize a normal by allowing for non-zero skewness. A univariate skew-normal distribution is defined by three parameters location $\xi\in\mathbb{R}$, scale $\sigma>0$ and skew parameter $\alpha\in\mathbb{R}$ and has the following (O’Hagan and Leonard,, 1976) Probability Density Function (PDF)

where $\phi$ and $\Phi$ are the PDF and Cumulative Distribution Function (CDF), respectively, of the standard univariate Normal distribution. Over the years many generalisations of this distribution were proposed, in particular Arellano and Azzalini, (2006) provided a unification of those generalizations in a single and tractable multivariate Unified Skew-Normal distribution. This distribution satisfies closure properties for marginals and conditionals and allows more flexibility due the introduction of additional parameters.

Consider $n$ input points $X=\{{\bf x}_{i}:i=1,\ldots,n\}$, with ${\bf x}_{i}\in\mathbb{R}^{d}$, and a data-dependent matrix $W\in\mathbb{R}^{m\times n}$. We define an affine probit likelihood as

where $\Phi_{m}({\bf x}):=\Phi_{m}({\bf x};I_{m})$ is the standard $m$-variate Gaussian CDF evaluated at ${\bf x}\in\mathbb{R}^{m}$ with identity covariance matrix. Note that this likelihood model includes the classic GP probit classification model (Rasmussen and Williams,, 2006) with binary observations $y_{1},\ldots,y_{n}\in\{0,1\}$ encoded in the matrix $W=\text{diag}(2y_{1}-1,\ldots,2y_{n}-1)$, where $m=n$. Moreover, as we will show in Corollary 1, the likelihood in (4) is equal to the preference likelihood for a particular choice of $W$. This model however also allows to seamlessly mix classification and preference information as we will show below. Here we prove how the skew-normal distribution is connected to the affine probit likelihood, extending a a result proved in (Durante,, 2019, Th.1 and Co.4) for the parametric setting for a standard probit likelihood.²²2(Durante,, 2019, Th.1 and Co.4) assumes that the matrix $W$ is diagonal, but the same results can straightforwardly extend to generic $W$.

All the proofs are in Appendix B. We now prove that, a-posteriori, for a new test point $\bf x$, the function $f(\bf x)$ is SkewGP distributed under the affine probit likelihood in (4).

This is the main result of the paper and allows us to show that, in the case of preference learning, we can compute exactly the posterior and, therefore, Laplace’s approximation is not necessary.

We provide a one-dimensional illustration of the difference between GPL and SkewGP. Consider the non-linear function $g(x)=cos(5x)+e^{-\frac{x^{2}}{2}}$ which has a global maximum at $x=0$. We assume we can only query this function through pairwise comparisons. We generate 7 pairwise random comparisons: the query point $x_{i}$ is preferred to $x_{j}$ (that is $x_{i}\succ x_{j}$) if $g(x_{i})>g(x_{j})$. Figure 3(top-left) shows $g(x)$ and the location of the queried points.⁶⁶6The preferences between the queried points are $1.25\succ-1.8$, $-1.23\succ 1.25$, $0.18\succ-1.23$, $0.18\succ-2.52$, $-2.52\succ 2.18$, $-1.8\succ-0.5$, $-1.8\succ 0.67$. Figure 3(bottom-left) shows the predicted posterior preference function $f(x)$ (and relative 95% credible region) computed according to GPL and SkewGP. Both the methods have the same prior: a GP with zero mean and RBF covariance function (the hyperparameters are the same for both methods and have been set equal to the values that maximise Laplace’s approximation to the marginal likelihood, $l=0.35$ and $\sigma^{2}=0.02$). Therefore, the only difference between the two posteriors is due to the Laplace’s approximation. The true posterior (SkewGP) of the preference function is skewed, this can be seen from the density plot for $f(-0.51)$ in Figure 3(top-right). Figure 3(bottom-right) shows an example, $f(0.19)$, where SkewGP and Laplace’s approximation differ significantly: Laplace’s approximation heavily underestimates the mean and the support of the true posterior (SkewGP) also evident from Figure 3(bottom-left). These differences determine the poor performance of PBO based on GPL as we will see in the next sections.

In sequential BO, our objective is to seek a new data point ${\bf x}$ which will allow us to get closer to the maximum of the target function $g$. Since $g$ can only be queried via preferences, this is obtained by optimizing w.r.t. ${\bf x}$ a dueling acquisition function $\alpha({\bf x},{\bf x}_{r})$, where ${\bf x}_{r}$ (reference point) is the best point found so far, that is the point that has the highest probability of winning most of the duels (given the observed data $\mathcal{D}$) and, therefore, it is the most likely point maximizing $g$.⁷⁷7 By optimizing the acquisition function $\alpha({\bf x},{\bf x}_{r})$, we aim to find a point that is better than ${\bf x}_{r}$ (also considering the trade-off between exploration and exploitation). After computing the optimum of the the acquisition function, denoted with ${\bf x}_{n}$, we query the black-box function for ${\bf x}_{n}\stackrel{{\scriptstyle?}}{{\succ}}{\bf x}_{r}$. If ${\bf x}_{n}\succ{\bf x}_{r}$ then ${\bf x}_{n}$ becomes the new reference point (${\bf x}_{r}$) for the next iteration. We consider three pairwise modifications of standard acquisition functions: (i) Upper Credible Bound (UCB); (ii) Thompson sampling (TH); (iii) Expected Improvement Info Gain (EIIG).

UCB: The dueling UCB acquisition function is defined as the upper bound of the minimum width $\gamma$% (in the experiments we use $\gamma=95$) credible interval of $f({\bf x})-f({\bf x}_{r})$.

TH: The dueling Thompson acquisition function is $f_{j}({\bf x})-f_{j}({\bf x}_{r}),$ where $f_{j}$ is a sampled function from the posterior.

EIIG: We now propose the dueling EIIG that is the combination of the expected probability of improvement (in log-scale) and the dueling information gain:

where $IG({\bf x},{\bf x}_{r})=h\left(E_{f\sim p(f|\mathcal{D})}\left(\Phi\left(\tfrac{f({\bf x})-f({\bf x}_{r})}{\sqrt{2}\sigma}\right)\right)\right)-E_{f\sim p(f|\mathcal{D})}\left(h\left(\Phi\left(\tfrac{f({\bf x})-f({\bf x}_{r})}{\sqrt{2}\sigma}\right)\right)\right),$ with $h(p)=-p\log(p)-(1-p)\log(1-p)$ being the binary entropy function of $p$. This definition of dueling information gain is an extension to preferences of the information gain, formulated for GP classifiers in (Houlsby et al.,, 2011). This last acquisition function allows us to balance exploration-exploitation by means of the nonnegative scalar $k$ (in the experiments we use $k=0.1$ (more exploration) and $k=0.5$). To compute these acquisition functions, we make explicitly use of the generative capabilities of our SkewGP surrogated model as well as the efficiency of sampling the learned preference function.

Note that, González et al., (2017) use different acquisition functions based on the Copland’s score (to increase exploration). Moreover, they optimize $\alpha({\bf x}_{a},{\bf x}_{b})$ with respect to both ${\bf x}_{a},{\bf x}_{b}$, while ${\bf x}_{b}$ is fixed and equal to ${\bf x}_{r}$ in our setting. SkewGP can easily be employed as surrogated model in González et al., (2017) PBO setting (and improve their performance due to the limits of the Laplace’s approximation). We have focused on the above acquisition functions (UCB, TH, EIIG) because they can easily be computed – instead the Copland’s score requires to numerically compute an integral with respect to ${\bf x}$.

In this section we present numerical experiments to validate our PBO-SkewGP and compare it with PBO based on the Laplace’s approximation (PBO-GPL).

First, we consider again the maximization of $g(x)=cos(5x)+e^{-\frac{x^{2}}{2}}$ and the same $7$ initial preferences used in Section 4.

We run PBO for 4 iterations and, at each step, we query the point that maximises the UCB of GPL. We also compute the true posterior and the true maximum of UCB using SkewGP for comparison. Both the methods have the same prior: a GP with zero mean and RBF covariance function (the hyperparameters are fixed to the same values for both methods, that is the values that maximise Laplace’s approximation to the marginal likelihood, $l=0.35$ and $\sigma^{2}=0.02$). Therefore, the only difference between the two posteriors is due to the Laplace’s approximation.

Figure 4 shows, for each iteration, a row with three plots. The left plot reports $g(x)$ and the queried points ($x_{r}$ is in orange). The central plot shows the GPL and SkewGP posterior predictive means of $f(x)-f(x_{r})$ and the relative 95% credible intervals. The maximum of each UCB is showed with a star marker. The right plot shows the skewness statistics for the SkewGP predictive posterior distribution of $f(x)-f(x_{r})$ as a function of $x$, defined as:

with $\mu:=\operatorname{E}\left[f(x)\right]$, and computed via Monte Carlo sampling from the posterior.

Figure 4 shows that the Laplace approximation for $f(x)-f(x_{r})$ is much worse than SkewGP. Moreover, the posterior of $f(x)-f(x_{r})$ is heavily skewed. The maximum magnitude of $SS(f(x)-f(x_{r}))$ is $-1.3$ (see the relative marginal posteriors in Figure 5(top)).

Note from Figure 4(1st-row, central) that, while the maximum of UCB for GPL and SkewGP almost coincides in the initial iteration, they significantly differ in the second iteration, Figure 4(2nd-row, central): SkewGP’s UCB selects a point very close to the global maximum, while GPL explores the border of the search space. This again happens in the subsequent two iterations, see 4(3nd-row and 4th-row, central).

This behavior is neither special to this trial nor to the UCB acquisition function, we repeat this experiment 20 times starting with $10$ initial (randomly selected) duels and using all the three acquisition functions. We compare GPL versus SkewGP with fixed kernel hyperparameters (same as above) so that the only difference between the two algorithms is in the computation of the posterior.⁸⁸8In our implementation we compute the acquisition functions via Monte Carlo sampling (2000 samples). We report the average (over the 20 trials) performance, which is defined as the value of $g$ evaluated at the current optimal point ${\bf x}_{r}$, considered to be the best by each method at each one of the 100 iterations.

The results are showed in Figure 5(bottom): SkewGP always outperforms GPL. This is only due to Laplace’s approximation (hyper-parameters are the same). In 1D, the differences are smaller for the Thompson (TH) acquisition function (due to the “noise” from the random sampling step). However, SkewGP-TH converges faster.

In this work we have shown that is possible to perform exact preferential Bayesian optimization by using Skew Gaussian processes. We have demonstrated that in the setting of preferential BO: (i) the Laplace’s approximation is very poor; (ii) given the strong skewness of the posterior, any approximation of the posterior that relies on a symmetric distribution will result to sub-optimal predictive performances and, therefore, slower convergence in PBO. We have also shown that we can extend this model to deal with mixed preferential-categorical Bayesian optimisation, while still providing the exact posterior. We envisage an immediate extension of our current approach. Many optimisation applications are subject to safety constraints, so that inputs cannot be freely chosen from the entire input space. This leads to so-called safe Bayesian optimization (Gotovos et al.,, 2015), that has been extended to safe preferential BO in (Sui et al.,, 2018). We plan to solve this problem exactly using SkewGP.