{'fraction_non_alphanumeric': 0.04662565191593017, 'fraction_numerical': 0.015083851222635145, 'mean_word_length': 4.3343328335832085, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 13, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Preferential Bayesian optimisation with Skew Gaussian Processes

Bayesian optimisation (BO) is a very effective approach for sequential black-box optimization where direct queries of the objective function are expensive. However, there are cases where the objective function can only be accessed via preference judgments, such as"this is better than that"between two candidate solutions (like in A/B tests or recommender systems). The state-of-the-art approach to Preferential Bayesian Optimization (PBO) uses a Gaussian process to model the preference function and a Bernoulli likelihood to model the observed pairwise comparisons. Laplace's method is then employed to compute posterior inferences and, in particular, to build an appropriate acquisition function. In this paper, we prove that the true posterior distribution of the preference function is a Skew Gaussian Process (SkewGP), with highly skewed pairwise marginals and, thus, show that Laplace's method usually provides a very poor approximation. We then derive an efficient method to compute the exact SkewGP posterior and use it as surrogate model for PBO employing standard acquisition functions (Upper Credible Bound, etc.). We illustrate the benefits of our exact PBO-SkewGP in a variety of experiments, by showing that it consistently outperforms PBO based on Laplace's approximation both in terms of convergence speed and computational time. We also show that our framework can be extended to deal with mixed preferential-categorical BO, typical for instance in smart manufacturing, where binary judgments (valid or non-valid) together with preference judgments are available.


Introduction
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). 1 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.

Background on skew-normal distributions and skew-Gaussian processes
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 ξ ∈ R, scale σ > 0 and skew parameter α ∈ R and has the following (O'Hagan and Leonard, 1976) Probability Density Function (PDF) where φ and Φ 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.

Unified Skew-Normal distribution
A vector z ∈ R p is distributed as a multivariate Unified Skew-Normal distribution with latent skewness dimension s, z ∼ SUN p,s (ξ ξ ξ , Ω, ∆, γ γ γ, Γ), if its probability density function (Azzalini, 2013, Ch.7) is: where φ p (z − ξ ξ ξ ; Ω) represents the PDF of a multivariate Normal distribution with mean ξ ξ ξ ∈ R p and covariance Ω = D ΩΩ D Ω ∈ R p×p , withΩ being a correlation matrix and D Ω a diagonal matrix containing the square root of the diagonal elements in Ω. The notation Φ s (a; M) denotes the CDF of N s (0, M) evaluated at a ∈ R s . The parameters γ γ γ ∈ R s , Γ ∈ R s×s , ∆ p×s of the SUN distribution are related to a latent variable that controls the skewness, in particular ∆ is called Skewness matrix. Figure 1: Density plots for SUN 1,s (0, 1, ∆, γ, Γ). For all plots Γ is a correlation matrix, γ = 0, dashed lines are the contour plots of y ∼ N 1 (0, 1).
The PDF (1) is well-defined provided that the matrix i.e., M is positive definite. Note that when ∆ = 0, (1) reduces to φ p (z − ξ ξ ξ ; Ω), i.e. a skew-normal with zero skewness matrix is a normal distribution. Moreover we assume that Φ 0 (·) = 1, so that, for s = 0, (1) becomes a multivariate Normal distribution.  Figure 1 shows the density of a univariate SUN distribution with latent dimensions s = 1 (a1) and s = 2 (a2). The effect of a higher latent dimension can be better observed in bivariate SUN densities as shown in Figure 2. The contours of the corresponding bivariate normal are dashed. We also plot the skewness directions given byΩ −1 ∆.
For what follows however it is important to know that (see, e.g., Azzalini, 2013, Ch.7) the distribution is closed under marginalization and conditioning. We have reviewed these results in Appendix A together with an additive representation of the SUN that is useful for sampling from the posterior.
A SkewGP Benavoli et al. (2020) is a generalization of a skew-normal distribution to a stochastic process. We consider here a location function ξ : R d → R, a scale (kernel) function Ω : R d × R d → R, a skewness vector function ∆ : R d → R s and the parameters γ γ γ ∈ R s , Γ ∈ R s×s . A real function f : R d → R is a SkewGP with latent dimension s, if for any sequence of n points x 1 , . . . , x n ∈ R d , the vector [ f (x 1 ), . . . , f (x n )] ∈ R n is skew-normal distributed with parameters γ γ γ, Γ and location, scale and skewness matrices, respectively, given by .
The skew-normal distribution is well defined if the matrix M = Γ ∆(X) ∆(X) T Ω(X,X) is positive definite for all X = {x 1 , . . . , x n } ⊂ R d and for any n. In that case we write f ∼ SkewGP s (ξ , Ω, ∆, γ, Γ). Benavoli et al. (2020) shows that this is a well defined stochastic process. In the next section we connect this stochastic process to preference learning.

SkewGP and affine probit likelihood
Consider n input points X = {x i : i = 1, . . . , n}, with x i ∈ R d , and a data-dependent matrix W ∈ R m×n . We define an affine probit likelihood as where Φ m (x) := Φ m (x; I m ) is the standard m-variate Gaussian CDF evaluated at x ∈ 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 , . . . , y n ∈ {0, 1} encoded in the matrix W = diag(2y 1 − 1, . . . , 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.
Theorem 1. Let us assume that f (x) is GP distributed with mean function ξ (x) and covariance function , Ω(x, x )), and consider the likelihood p(W | f (X)) = Φ m (W f (X)) where W ∈ R m×n . The posterior distribution of f (X) is a SUN: where, for simplicity of notation, we denoted ξ (X), Ω(X, X) as ξ ξ ξ , Ω and Ω = D ΩΩ D Ω .
All the proofs are in Appendix B. We now prove that, a-posteriori, for a new test point x, the function f (x) is SkewGP distributed under the affine probit likelihood in (4).

Theorem 2. Let us assume a GP prior f
with W ∈ R m×n , then a-posteriori f is SkewGP with mean function ξ (x), covariance function Ω(x, x ), skewness function ∆(x, X) = Ω(x, X)W T , andγ,Γ as in (5). 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.

Exact preference learning
We now apply results of Theorem 1 and Theorem 2 to the case of preference learning. For two different inputs v k , u k ∈ X, a pairwise preference v k u k is observed, where v k u k expresses the preference of the instance v k over u k . A set of m pairwise preferences is given and denoted as D = {v k u k : k = 1, . . . , m}.
Likelihood. We assume that there is an underlying hidden function f : R d → R which is able to describe the observed set of pairwise preferences D. Specifically, given a preference v k u k , then the function f is To allow tolerances to model noise, we assume that value of the hidden function f is corrupted by a Gaussian noise with zero mean and variance σ 2 .
We use the likelihood introduced in Chu and Ghahramani (2005), which is the joint probability distribution of observing the preferences D given the values of the function f at X, i.e., For identifiability reasons, without loss of generality, we set σ 2 = 1 2 . 2 Posterior. The posterior distribution of the values of the hidden function f at all x ∈ X given the observations D is then: In state-of-art PBO (González et al., 2017), a Laplace's approximation of the posterior p( f (X) | D) is used to construct the acquisition function. The following Corollary shows that the posterior p( f (X) | D) is distributed as a SkewGP.
Then the posterior of f (X) is given by (5).
In PBO, in order to compute the acquisition functions, we must be able to draw efficiently independent samples from the posterior in Theorem 2.
Proposition 1. Given a test point x, posterior samples of f (x) can be obtained as: 2 Equivalently, we instead estimate the kernel variance.
Note that sampling U 0 can be achieved efficiently with standard methods, however using standard rejection sampling for the variable U 1 would incur in exponentially growing sampling time as the dimension m increases. Here we use the recently introduced sampling technique linear elliptical slice sampling (liness,Gessner et al. (2019)) which improves Elliptical Slice Sampling (ess, Murray et al. (2010)) for multivariate Gaussian distributions truncated on a region defined by linear constraints. In particular this approach derives analytically the acceptable regions on the elliptical slices used in ess and guarantees rejection-free sampling. Since lin-ess is rejection-free, 3 we can compute exactly the computation complexity of (15): O(n 3 ) with storage demands of O(n 2 ). SkewGPs have similar bottleneck computational complexity of full GPs. Finally, observe that U 1 does not depend on x and, therefore, we do not need to re-sample U 1 to sample f at another test point x . This is fundamental because acquisition functions are functions of x and, we need to optimize them in PBO.
Marginal likelihood. Here, we follow the usual GP literature (Rasmussen and Williams (2006)) and we consider a zero mean function ξ (x) = 0 and a parametric covariance kernel Ω(x, x ) indexed by θ ∈ Θ.Typically θ contains lengthscale parameters and a variance parameter. For instance, for the RBF kernel The parameters θ are chosen by maximizing the marginal likelihood, that for SkewGP is provided hereafter.
then the marginal likelihood of the observations D is withγ γ γ,Γ defined in Theorem 2 (they depend on θ ).
If the size of W is too large the evaluation of Φ m could become infeasible, therefore here we use the approximation introduced in Benavoli et al. (2020), see inequality in (9) where B 1 , . . . , B b are a partition of the training dataset into b random disjoint subsets, |B i | denotes the number of observations in the i-th element of the partition,γ γ γ B i ,Γ B i are the parameters of the posterior computed using only the subset B i of the data (in the experiments |B i | = 30). Details about the routine we use to compute Φ |B i | (·) and the optimization method we employ to maximise the lower bound in (9) are in Appendix C.

Mixed classification and preference information
Consider now a problem where we have two types of information: whether a certain instance is preferable to another (preference-like observation) and whether a certain instance is attainable or not (classificationlike observation). Such situation often comes up in industrial applications. For example imagine a machine that produces a product whose final quality depends on certain input parameters. Assume now that certain values of the input parameters produce no product. In this case we might want to evaluate the quality of the product with binary comparisons (preferences) along with a binary class that indicates whether the input configuration is valid. By using only a preference likelihood or a classification likelihood we would not be using all information. In this case, observations are in fact pairs and the space of possibility is where v k and u k are respectively the current and reference input. We propose a new likelihood function to model the above setting, which is defined as follows: It is then immediate to write the above likelihood in the form (4) and, therefore, use both sources of information. We associate to each point x i a binary output y i ∈ {0, 1} where the class 0 denotes a nonvalid output. In case of valid output, we assume that we conducted m comparisons obtaining the couples D = {v k u k : k = 1, . . . , m} where v k u k expresses the preference of the instance v k over u k . The likelihood is then a product of two independent probit likelihood functions . . .
Since we assume that the two likelihood are independent we can compute the overall likelihood: with W = W class W pref ∈ R (n+m)×n , where W pref is the matrix of preferences defined as in Corollary 1. Therefore, the results in Section 3 still holds in this mixed setting.

Comparison SkewGP vs. Laplace's approximation
We provide a one-dimensional illustration of the difference between GPL and SkewGP. Consider the nonlinear function g(x) = cos(5x) + e − 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 x j ) if g(x i ) > g(x j ). Figure 3(top-left) shows g(x) and the location of the queried points. 5 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 σ 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.

Dueling acquisition functions
In sequential BO, our objective is to seek a new data point 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 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: , 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 α(x a , x b ) with respect to both x a , x b , while x b is fixed and equal to 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 x.

Numerical experiments
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 − 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 σ 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: , 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. 7 We report the average (over the 20 trials) performance, which is defined as the value of g evaluated at the current optimal point 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.

Optimization benchmarks
We have considered for g(x) the same four benchmark functions used in (González et al., 2017): 'Forrester' (1D), 'Six-Hump Camel' (2D), 'Gold-Stein' (2D) and 'Levy' (2D), and additionally 'Rosenbrock5' (5D) and 'Hartman6' (6D). These are minimization problems. 8 Each experiment starts with 10 initial (randomly selected) duels and a total budget of 100 duels are run. Further, each experiment is repeated 20 times with different initialization (the same for all methods) as in (González et al., 2017). We compare PBO based on GPL versus SkewGP using the three different acquisition functions described in Section 5: UCB, EIIG (k = 0.1 and 0.5) and Thompson (denoted as TH). As before, we show plots of #iterations versus g(x r ). In these experiments we optimize the kernel hyperparameters by maximising the marginal likelihood for 7 In our implementation we compute the acquisition functions via Monte Carlo sampling (2000 samples). 8 We converted them into maximizations so that the acquisition functions in Section 5 are well-defined.   both GPL and SkewGP. Figure 6 reports the performance of the different methods. Consistently across all benchmarks PBO-SkewGP outperforms PBO-GPL no matter the acquisition function. PBO-SkewGP has also a lower computational burden as showed in the Table at the bottom of Figure 6 that compares the median (over 80 trials, that is 20 trials times 4 acquisition functions) computational time per 100 iterations in seconds (on a standard laptop).

Mixed preferential-categorical BO
We examine now situations where certain unknown values of the inputs produce no-output preference and address it as a mixed preferential-categorical BO as described in Section 3.2. We consider again g(x) = cos(5x) + e − x 2 2 and assume that any input x ≤ −0.2 produces a non-valid output. Figure 7(left) shows g(x), the location of the queried points 9 and the no-valid zone (in grey). Figure 7(right) shows the predicted posterior preference function f (x) (and relative 95% credible region) computed according to SkewGP. We can see how SkewGP learns the non-valid region: the posterior mean is negative for x −0.2 and positive otherwise (the oscillations of the mean for x −0.2 capture the preferences). This is consistent with the likelihood in (10).
We consider now the following benchmark optimization problem proposed in (Sasena et al., 2002): min 2 + 1 100 (x 2 − x 2 1 ) 2 + (1 − x 1 ) 2 + 2(2 − x 2 ) 2 + 7 sin( 1 2 x 1 ) sin( 7 10 x 1 x 2 ) We also assume that both the objective function and h(x) are unknown. The goal is to find the minimum: x = [2.7450 2.3523] with optimal cost −1.1743. Figure 8(left) shows the level curves of the objective function, the non-valid zone (grey bands) and the location of the minimizer (red star). We compare two approaches. The first approach uses PBO-SkewGP that minimizes the objective plus a penalty term, 10 8 max(0, h(x)) 2 , that is non-zero in the non-valid region. Adding a penalty for non-valid inputs is the most common approach to deal with this type of problems. The second approach uses a SkewGP based mixed preferentialcategorical BO (SkewGP-mixed) that accounts for the valid/non-valid points as in Section 3.2. Figure  8(right) shows the performance of the two compared methods, which is consistent across the three different acquisition functions: SkewGP-mixed converges more quickly to the optimum. This confirms that modelling directly this type of problems via the mixed preferential-categorical likelihood in (10) enables the model to fully exploit the available information. Also in this case, SkewGP allows us to compute the corresponding posterior exactly.

Conclusions
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 con-straints, 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.

A Background on the Skew-Normal distribution
In this section we provide more details on the skew-normal distribution. See, e.g. Azzalini (2013) for a complete reference on skew-normal distributions and their parametrizations.

A.1 Additive representations
The role of the latent dimension s can be briefly explained as follows. Consider a random vector (2) and define y as the vector with distribution (x 1 | x 0 + γ γ γ > 0). The density of y can be written as where the first equality comes from a basic property of conditional distributions, see, e.g. (Azzalini, 2013, Ch. 1.3), and the second equality is a consequence of the multivariate normal conditioning properties. Then we have that z = ξ ξ ξ + D Ω y ∼ SUN p,s (ξ ξ ξ , Ω, ∆, γ γ γ, Γ).
Proof. We can show that the representation x 0 , x 1 is equivalent to u 0 , u 1 . Define u 1 = x 0 and u 0 = x 1 − Then we have that u 0 and u 1 are independent. This can be verified by the fact that u 0 and u 1 are normally distributed with covariance Cov(u 0 , u 1 ) = 0 which can be verified with algebraic computations. Finally note that (u 0 + ∆Γ −1 u 1,−γ γ γ ) is distributed as (x 1 | x 0 γ γ γ > 0).
The additive representation introduced above is used in Section 2.4 to draw samples from the distribution.

A.3 Sampling from the posterior predictive distribution
Consider a test point x * and assume we have a sample from the posterior distribution f (X) | D.
Then by using the marginalization property introduced above we obtain the formula in (5).
B Proofs of the results in the paper Theorem 1 We aim to derive the posterior of f (X). The joint distribution of f (X), D is where we have omitted the dependence on X for easier notation. We note that From (12)-(13) and the definition of the PDF of the SUN distribution, we can easily show that we can rewrite (12) as a SUN distribution with updated parameters: Theorem 2 Consider the test point x ∈ R d and the vectorf = , and the predictive distribution is by definition We can then apply Lemma 1 withf and the likelihood p( Corollary 1 We can write the likelihood function as where V ∈ R m×n with V i, j = 1 if v i = x j and 0 otherwise and U ∈ R m×n with U i, j = 1 if u i = x j and 0 otherwise. Then we can apply Lemma 1 for the posterior distribution of f (X) and Theorem 2 for the posterior distribution of f at an unobserved point.

C.1 Laplace's approximation
The Laplace's approximation for preference learning was implemented as described in Chu and Ghahramani (2005). We use standard Bayesian optimisation to optimise the hyper-parameters of the kernel by maximising the Laplace's approximation of the marginal likelihood.

C.2 Skew Gaussian Process
To compute Φ |B i | (·) in (9), we use the routine proposed in Trinh and Genz (2015), that computes multivariate normal probabilities using bivariate conditioning. This is very fast. We optimise the hyper-parameters of the kernel by maximising the lower bound in (9) and we use simulated annealing.

C.3 Acquisition function optimisation
In sequential BO, our objective is to seek a new data point 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. x a dueling acquisition function α(x, x r ), where x r 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 D) and, therefore, it is the most likely point maximizing g. For both the models (Laplace's approximation (GPL) and SkewGP) we compute the acquisition functions α(x, x r ) via Monte Carlo sampling from the posterior. In fact, although for GPL some analytical formulas are available (for instance for UCB), by using random sampling (2000 samples with fixed seed) for both GPL and SkewGP removes any advantage of SkewGP over GPL due to the additional exploration effect of Monte Carlo sampling. Computing α(x, x r ) in this way is very fast for both the models (for SkewGP this is due to lin-ess). We then optimize α(x, x r ): (i) by computing α(x, x r ) for 5000 random generated value of x (every time we use the same random points for both SkewGP and GPL); (ii) the best random instance is then used as initial guess for L-BFGS-B.