Towards Practical Preferential Bayesian Optimization
with Skew Gaussian Processes

Preferential Bayesian optimization (BO) has been an attractive approach for solving problems where reliable feedback is limited to pairwise comparison, the so-called duels. This preference setting often appears in human-in-the-loop optimization problems such as visual design optimization (Koyama et al., 2020) and generative melody composition (Zhou et al., 2021) because it is easier for humans to judge which one is better than to give an absolute rating (Kahneman and Tversky, 2013). The system (i.e., the optimization method) in the human-in-the-loop optimization presents choices and receives preferential feedback interactively. To reduce the waiting time for users, the system is required to quickly present the new options to users by learning from the observed feedback information.

An important challenge in preferential BO, which uses the preferential Gaussian process (GP) model to represent preference structure, is that the posterior distribution is computationally intractable skew GP. The existing approaches to this difficulty are twofold: The first approach is the Gaussian approximation, such as Laplace approximation (LA) and expectation propagation (EP) (Chu and Ghahramani, 2005b, a). Actually, most of the preferential BO algorithms proposed so far are based on this Gaussian approximation (Brochu et al., 2010; Nielsen et al., 2015; González et al., 2017; Siivola et al., 2021; Fauvel and Chalk, 2021), leading to good computational efficiency. However, the predictive accuracy can be poor, as discussed in Section 3.1, because the Gaussian approximation ignores the skewness of the exact posterior represented as the skew GP (Benavoli et al., 2020, 2021a, 2021b). The second approach is to directly employ the skew GP model using the Markov chain Monte Carlo (MCMC) method (Benavoli et al., 2020, 2021a, 2021b). Benavoli et al. (2020, 2021a, 2021b) argued that the Gaussian approximation has poor estimation accuracy, and their MCMC-based estimation is effective in preferential BO. However, the experimental evidence was insufficient to show why the Gaussian approximation degrade the preferential BO.

The aim of this study is to build truly practical preferential BO with experimental evidence to support it. To this end, we first evaluate the accuracy of the Gaussian approximation-based methods, which have been most used in the preferential BO literature. Our experiments reveal that those methods fail to correctly predict the probability of duel outcomes, which is critical in the preferential BO setting. This result implies that it is important to consider the skewness of the true posterior, and then we also improve MCMC for skew GP estimation. Specifically, we present MC estimators that reduce the estimation variance for skew GP. Furthermore, we empirically clarify the high computational efficiency of the Gibbs sampling compared with an MCMC method employed in (Benavoli et al., 2020, 2021b). These improvements make the preferential BO algorithm with skew GP more attractive. On the other hand, the preferential BO algorithms using MCMC can still be unrealistic in preferential BO, where the computational time may become a bottleneck in many cases.

To build practical preferential BO, we further develop a simple and computationally efficient preferential BO method that incorporates the skewness of the true posterior in a randomized manner. The basic idea is to use a posterior additionally conditioned by a random sample from the original posterior itself, called hallucination, by which we show that computing acquisition functions (AFs) in the proposed method is computationally efficient while reflecting the skewness of the original posterior. Noteworthy, any powerful AFs developed so far in the standard BO literature (e.g., upper confidence bound (UCB) (Srinivas et al., 2010) and expected improvement (EI) (Mockus et al., 1978)) can be integrated, which improves the flexibility of the proposed method.

Our contributions are summarized as follows:

1. 1.

   Rethinking existing preferential BO algorithms: We reveal the poor accuracy of Gaussian approximation-based methods in the true posterior approximation, which implies the importance of considering skewness in the preferential BO setting. Motivated by this observation, we improve MCMC for skew GP estimation.

2. 2.

   Building practical preferential BO algorithm: We develop a simple and computationally efficient method, called hallucination believer, that can reflect the skewness of the original posterior and integrate any powerful AFs developed so far in the standard BO literature.

3. 3.

   Extensive empirical validations: Numerical experiments on 12 benchmark functions show that the proposed method achieves better or at least competitive performance in both terms of computational efficiency and sample complexity compared with Gaussian approximation-based preferential BO (González et al., 2017; Siivola et al., 2021; Fauvel and Chalk, 2021) and MCMC-based preferential BO (Benavoli et al., 2021a, b), respectively.

Our experimental codes are publicly available at https://github.com/CyberAgentAILab/preferentialBO.

First, we investigate the limitation of the Gaussian approximation, which is the inaccurate predictive probability. Second, we tackle the improvement of MCMC-based estimation for skew GP. We show the practical effectiveness of Gibbs sampling for skew GP. Furthermore, we newly derive MC estimators, showing that the estimation variance is lower than the estimators used in (Benavoli et al., 2021a, b).

Although we tackled the speed-up of MCMC-based estimation, its computational time can still be unrealistic in preferential BO. Furthermore, Gaussian approximation often suffers poor prediction performance ignoring the skewness of the posterior. In this work, we propose a computationally efficient preferential BO method while considering the skewness in a randomized manner.

Although dueling bandit (Yue et al., 2012) considers online learning from the duels, it does not consider the correlation between arms in general, as reviewed in (Sui et al., 2018). In contrast, preferential BO (Brochu et al., 2010) aims for more efficient optimization using the preferential GP, which can capture a flexible correlational structure in an input domain. Since the true posterior is computationally intractable, most prior works employed Gaussian approximation-based preferential GP models (Chu and Ghahramani, 2005b, a). Typically, preferential BO sequentially duels $\boldsymbol{x}_{t}^{(1)}$, which is the most exploitative point (e.g., maxima of the posterior mean), and $\boldsymbol{x}_{t}^{(2)}$, which is the maxima of some AFs.

The design of AFs for $\boldsymbol{x}_{t}^{(2)}$ is roughly twofold. First approaches (Brochu et al., 2010; Siivola et al., 2021) use the AFs for standard BO, which assumes that the direct observation $y$ can be obtained (e.g., EI (Mockus et al., 1978)). However, this approach selects similar duels repeatedly since the preferential GP model’s variance is hardly decreased due to the weakness of the information obtained from the duel. Second approach considers the relationship between $\boldsymbol{x}_{t}^{(1)}$ and $\boldsymbol{x}_{t}^{(2)}$. For example, bivariate EI (Nielsen et al., 2015) uses the improvement from $f(\boldsymbol{x}_{t}^{(1)})$, and MUC (González et al., 2017; Fauvel and Chalk, 2021) employed the uncertainty of the probability that $\boldsymbol{x}_{t}^{(2)}$ wins $\boldsymbol{x}_{t}^{(1)}$. Considering the correlation to $\boldsymbol{x}_{t}^{(1)}$, these AFs are expected to explore appropriately. However, as discussed in Section 3.1, Gaussian approximation (LA and EP) can be inaccurate for evaluation of the joint distribution of $f(\boldsymbol{x}_{t}^{(1)})$ and $f(\boldsymbol{x}_{t}^{(2)})$.

Instead of using Gaussian approximation, Benavoli et al. (2021a) directly estimates the true posterior using MCMC. However, MCMC-based AFs require heavy computational time, which can be a bottleneck in practical applications.

Another well-used AF for preferential BO is TS (Russo and Van Roy, 2014). González et al. (2017) proposed to select $\boldsymbol{x}_{t}^{(1)}$ by TS, and Siivola et al. (2021); Sui et al. (2017) proposed to select both $\boldsymbol{x}_{t}^{(1)}$ and $\boldsymbol{x}_{t}^{(2)}$ by TS, both of which are based on Gaussian approximation. Benavoli et al. (2021a) also proposed to select $\boldsymbol{x}_{t}^{(2)}$ by TS, which is generated from skew GP. Although the merit of Gaussian approximation is low computational time, the posterior sampling from skew GP only once is sufficiently fast (comparable to EP). Therefore, the advantage of Gaussian approximation for TS-based AF is limited. In Section 6, we will show that TS can cause over-exploration when the input dimension is large in preferential BO, as with the standard BO reviewed in (Shahriari et al., 2016).

We investigate the effectiveness of the proposed preferential BO method through comprehensive numerical experiments. We employed the 12 benchmark functions. In this section, we show the results for the 8 functions, called the Branin, Holder table, Bukin, Eggholder, Ackley, Harmann3, Hartmann4, and Hartmann6; others are shown in Appendix G.

We performed HB combining EI Mockus et al. (1978) and UCB Srinivas et al. (2010), denoted as HB-EI and HB-UCB, respectively. We employed the baseline methods, LA-EI(Brochu et al., 2010), EP-EI Siivola et al. (2021), EP-MUC Fauvel and Chalk (2021), EP-TS-MUC (González et al., 2017)¹¹1This method was originally proposed as “dueling TS,” whose name is same as (Benavoli et al., 2021a). To distinguish them, we denote TS-MUC since this method selects the first input by TS and the second input by MUC. Furthermore, although González et al. (2017) used the GPC model, we use the EP-based preferential GP model to concentrate on the difference from the AF., DuelTS Benavoli et al. (2021a), DuelUCB Benavoli et al. (2021a), and EIIG²²2 Benavoli et al. (2021a, b) proposed EIIG that is EI minus information gain (IG), which is probably a typo since both EI and IG should be large. Thus, we used a modified one that EI plus IG. Benavoli et al. (2021a), where DuelTS, DuelUCB, and EIIG were based on skew GP and other methods with Suffix “EP” and “LA” employed EP and LA, respectively. The first input $\boldsymbol{x}^{(1)}_{t}$ for EP-EI, EP-MUC, and LA-EI is the maxima of the posterior mean over the training duels, and other methods, including HB, use the winner so far as $\boldsymbol{x}^{(1)}_{t}$, following (Benavoli et al., 2020). For preferential GP models, we use RBF kernel with automatic relevance determination (Rasmussen and Williams, 2005), whose lengthscales are selected by marginal likelihood maximization per 10 iterations, and set fixed noise variance $\sigma^{2}_{\rm noise}=10^{-4}$. We use marginal likelihood estimation by LA, which is sufficiently precise if $\sigma^{2}_{\rm noise}$ is small, as discussed in (Kuss and Rasmussen, 2005, Section 7.2). We show the accuracy of the marginal likelihood estimations of LA and EP in Appendix G. Other detailed settings are shown in Appendix F.

Figure 4 show the results for (a) computational time and (b) iterations. As a performance measure, we used the regret defined as $f(\boldsymbol{x}_{*})-f(\tilde{\boldsymbol{x}}_{t})$, where $\tilde{\boldsymbol{x}}_{t}$ is a recommendation point at $t$-th iteration, for which we used $\boldsymbol{x}^{(1)}_{t}$. We report the mean and standard error of the regret over 10 random initialization, where the initial duel is $3d$ uniformly random input pair.

First, we focus on the comparison against Gaussian approximation-based methods. LA-EI is the worst method due to the poor approximation accuracy. Note that, since we set the plotted interval for the small regret to focus the differences between HB and baseline methods except for LA-EI, LA-EI is out of the range in Hartmann 3, 4, and 6 functions. The regret of EP-EI, which selects the duel between the maximums of the posterior mean and EI, is often stagnated due to over-exploitation. The regrets of EP-MUC and EP-TS-MUC also stagnate in the Bukin and Ackley functions since EP is often inaccurate in predicting the joint distribution of several outputs, as shown in Section 3.1. Furthermore, EP-TS-MUC deteriorates in the high-dimensional functions such as Hartmann 4 and 6 functions due to over-exploration from the nature of TS. Overall, incorporating the skewness of the true posterior, HB-EI and HB-UCB outperformed Gaussian approximation-based methods in Figure 4 (b). Furthermore, the same tendency can also be confirmed in Figure 4 (a) since the computational time of HB-EI and HB-UCB is comparable with Gaussian approximation.

Second, DuelUCB and EIIG are MCMC-based AFs based on skew GP. Therefore, although they often show rapid convergence in Figure 4(b) because of the use of skew GP, they show slow convergence in Figure 4(a) due to the heavy computational time of MCMC. Hence, HB-EI and HB-UCB are comparable to and outperformed MCMC-based AFs in Figure 4(b) and (a), respectively, which suggest the practical effectiveness of the proposed methods.

Last, since DuelTS requires only one sample path, its performance is not largely changed in Figure 4 (a) and (b). Although the regret of DuelTS converges rapidly in low-dimensional functions, it converges considerably slowly in relatively high-dimensional functions, such as Hartmann4 and Hartmann6, in which the proposed methods clearly outperformed DuelTS. In our experiments, we confirmed that TS results in over-exploration in relatively high-dimensional functions, which is the well-known problem of TS for the standard BO, as reviewed in (Shahriari et al., 2016).

Towards building a more practical preferential BO, this work has made the following important progresses. We first revealed the poor performance of the Gaussian approximation for skew GP, such as Laplace approximation and expectation propagation, despite being the method most widely used in preferential BO. This result motivated us to improve the MCMC-based estimation of the skew GP, for which we showed the practical effectiveness of Gibbs sampling and derived the low variance MC estimators. Finally, we developed a new preferential BO method that achieves both high computational efficiency and low sample complexity, and then verified its effectiveness through extensive numerical experiments. The regret analysis of the proposed method is one of the interesting future works.

This work was supported by MEXT KAKENHI 21H03498, 22H00300, MEXT Program: Data Creation and Utilization-Type Material Research and Development Project Grant Number JPMXP1122712807, and JSPS KAKENHI Grant Number JP21J14673.