1 Introduction
In many practical situations a user is able to select the best options among a finite set of choices, however they are unable to state explicitly the motivations for their choices. A notable example is in industrial applications where the manufactured product has to satisfy several qualitative requirements that are known to trained staff, but such requirements were never expressed explicitly. In such cases, the definition of quantitative objectives would allow for an explicit multiobjective optimization which would lead to better options. However, measuring the objectives in a quantitative way is often technically difficult and costly. In this context, we would like to improve the quality of the manufactured product using directly the feedback provided by the user’s choices (“these products are better than those”), i.e. we would like to learn the “choice function” of the user and find the inputs that optimize this function. In this paper we propose a Bayesian framework to learn choice functions from a dataset of observed choices. Our framework learns a latent mapping of objectives that are consistent with the given choices, therefore we are also able to optimize them with a multiobjective Bayesian optimization algorithm.
2 Background
The main contributions of this paper leverage four topics: (1) Bayesian Optimisation (BO); (2) preferential BO; (3) multiobjective BO; (4) choice functions learning. In this section we briefly review the state of the art of each topic.
2.1 Bayesian Optimisation (BO)
BO [jones1998efficient] aims to find the global maximum of an unknown function which is expensive to evaluate. For a scalar realvalued function on a domain , the goal is to find a global maximiser . BO formulates this as a sequential decision problem – a tradeoff between learning about the underlying function (exploration) and capitalizing on this information in order to find the optimum (exploitation). BO relies on a probabilistic surrogate model, usually a Gaussian Process (GP) [rasmussen2006gaussian], to provide a posterior distribution over given a dataset of previous evaluations of . It then employs an acquisition function (e.g. Expected Improvement [jones1998efficient, mockus1978application], Upper Credible Bound [srinivas2009gaussian]) to select the next candidate option (solution) . While the true function is expensivetoevaluate, the surrogatebased acquisition function is not, and it can thus be efficiently optimized to compute an optimal candidate to be evaluated on . This process is repeated sequentially until some stopping criterion is achieved.
2.2 Preferential Bayesian Optimisation (PBO)
In many applications, evaluating can be either too costly or not always possible. In these cases, the objective function may only be accessed via preference judgments, such as “this is better than that” between two candidate options like in A/B tests or recommender systems (pairwise comparisons are usually called duels in the BO and bandits literature). In such situations, PBO [shahriari2015taking] can be used. This approach requires the agent to simply compare the final outcomes of two different candidate options and indicate which they prefer, that is the evaluation is binary either “better than” or “worse than” .
In the PBO context, the stateoftheart surrogate model is based on a method for preference learning developed in [ChuGhahramani_preference2005]. This method assumes that there is an unobservable latent function value associated with each training sample , and that the function values preserve the preference relations observed in the dataset, that is whenever “better than” . As in the BO setting, by starting with a GP prior on and, by using the likelihood defined in [ChuGhahramani_preference2005], we obtain a posterior distribution over given a dataset of preferences. This posterior distribution is not a GP and several approximations [ChuGhahramani_preference2005, houlsby2011bayesian] were proposed. In [houlsby2011bayesian], the authors showed that GP preference learning is equivalent to GP classification with a transformed kernel function. By using this reformulation, the authors easily derive two approximations for the posterior based on (i) the Laplace Approximation (LP) [mackay1996bayesian, williams1998bayesian]; (ii) Expectation Propagation (EP) [minka2001family]. The LP approximation was then used to develop a framework for PBO [shahriari2015taking] and a new acquisition function, inspired by Thomson sampling, was proposed in [gonzalez2017preferential]. More recently, [benavoli2020preferential]
showed that the posterior of GP preference learning is a Skew GP
[Benavoli_etal2020, benavoli2021]. Based on this exact model, the authors derived a PBO framework which outperformed both LP and EP based PBO. Although in this work we focus on GPs as surrogate model, it is worth to mention alternative approaches for PBO developed by [BDG08, pmlrv32zoghi14, sadigh2017active, bemporad2019active].PBO was recently extended [siivola2020preferential] to the batch case by allowing agents to express preferences for a batch of options. However, as depicted in Figure 1 where an agent expresses preferences among 5 options, in batch PBO there can be only one batch winner. In fact, PBO assumes that two options are always comparable.^{1}^{1}1More precisely, the underlying GPbased model implies a total order and so two options may also be equivalent. When PBO is applied to the multiobjective case such as for instance , it is therefore assumed that the agent’s preferences are determined by a weighted combination of the objectives .
2.3 Multiobjective (MO) optimization
The goal of MO optimization is to identify the set of Pareto optimal options (solutions) such that any improvement in one objective means deteriorating another. Without loss of generality, we assume the goal is to maximize all objectives. Let be a vectorvalue objective function with , where is the number of objectives. We recall the notions of Pareto dominated options and nondominated set.
Definition 1 (Pareto dominate option).
Consider a set of options . An option is said to Pareto dominate another option , denoted as , if both the following conditions are true:

for all , ;

, such that .
Definition 2 (Nondominated set).
Among a set of options , the nondominated set of options are those that are not dominated by any member of , i.e.
Given the set of options , MO aims to find the nondominated set of options , called the Pareto set. The set of evaluations is called Pareto front.
MO BO have only be developed for standard (nonpreferential) BO, where multiobjectives can directly be evaluated. Many approaches rely on scalarisation to transform the MO problem into a singleobjective one, like ParEGO [knowles2006parego] and TSTCH [paria2020flexible]
(which randomly scalarize the objectives and use Expected Improvement and, respectively, Thompson Sampling).
[keane2006statistical] derived an expected improvement criterion with respect to multiple objectives. [ponweiser2008multiobjective] proposed an hypervolumebased infill criterion, where the improvements are measured in terms of hypervolume (of the Pareto front) increase. Other acquisition functions have been proposed in [emmerich2011hypervolume, picheny2015multiobjective, hernandez2016predictive, wada2019bayesian, belakaria2019max]. The most used acquisition function for MO BO is expected hypervolume improvement. In fact, maximizing the hypervolume has been shown to produce very accurate estimates
[zitzler2003performance, couckuyt2014fast, hupkens2015faster, emmerich2016multicriteria, yang2017computing, yang2019multi] of the Pareto front.2.4 Choice function
Individuals are often confronted with the situation of choosing between several options (alternatives). These alternatives can be goods that are going to be purchased, candidates in elections, food etc.
We model options, that an agent has to choose, as realvalued vectors and identify the sets of options as finite subsets of . Let denote the set of all such finite subsets of .
Definition 3.
A choice function is a setvalued operator on sets of options. More precisely, it is a map such that, for any set of options , the corresponding value of is a subset of (see for instance [aleskerov2007utility]).
The interpretation of choice function is as follows. For a given option set , the statement that an option is rejected from (that is, ) means that there is at least one option that an agent strictly prefers over . The set of rejected options is denoted by and is equal to . Therefore choice functions represent nonbinary choice models, so they are more general than preferences.
It is important to stress again that the statement implies there is at least one option that an agent strictly prefers over . However, the agent is not required to tell us which option(s) in they strictly prefer to . This makes choice functions a very easytouse tool to express choices. As depicted in Figure 2, the agent needs to tell us only the options they selected (in green) without providing any justification for their choices (we do not know which option in the green set dominates ).
By following this interpretation, the set can also be seen as the nondominated set in the Pareto sense for some latent function. In other words, let us assume that there is a latent vector function , for some dimension , which embeds the options into a space . The choice set can then be represented through a Pareto set of nondominated options. For example, in Fig. 2, . This approach was proposed in [pfannschmidt2020learning]
to learn choice functions. In particular, to learn the latent vector function, the authors devise a differentiable loss function based on a hinge loss. Furthermore, they add two additional terms to the loss function: (i) an
regularization term; (ii) a multidimensional scaling (MDS) loss to ensure that options close to each other in the inputs space will also be close in the embedding space. This loss function is then used to learn a (deep) multilayer perceptron to represent the embedding.
2.5 Contributions
In this work, we devise a novel multiobjective PBO based on choice functions. We follow the interpretation of choice functions as set function that select nondominated sets for an unknown latent function. First we derive a Bayesian framework to learn the function from a dataset of observed choices. This framework is based on a Gaussian Process prior on the unknown latent function vector.^{2}^{2}2Compared to the approach proposed in [pfannschmidt2020learning], the GPbased model is more sound – no multidimensional scaling is necessary – and it is a generative model. We then build an acquisition function to select the best next options to evaluate. We compare this method against an oracle that knows the true value of the latent functions and we show that, by only working with choice function evaluations, we converge to the same results.
3 Bayesian learning of Choice functions
In this work we consider options and, for , we model each latent function in the vector as an independent GP [rasmussen2006gaussian]:
(1) 
Each GP is fully specified by its kernel function , which specifies the covariance of the latent function between any two points. In all experiments in this paper, the GP kernel is Matern 3/2 [rasmussen2006gaussian].
3.1 Likelihood for general Choice functions
Having defined the prior on , we can now focus on the likelihood. We propose a new likelihood to model the observed choices of the agent. Given a set of observed choices , we are interested in learning a Paretoembedding coherent with this data in the sense that , where denotes the Pareto nondominated options in .
Assume that and let be the subset of indices of the options in , let be equal to and let the vector of dimensions of the latent space. Based on Definition 1, the choice of the agent expressed via implies that:
(2)  
(3) 
These equation express the conditions in Definition 1. Condition (2) means that, for each option , it’s not true ( stands for logical negation) that all options in are worse than , i.e. there is at least an option in which is better than . Condition (3) means that, for each option in , there is no better option in . This requires that the latent functions values of the options should be consistent with the choice function implied relations. Given , the likelihood function is one when (2)(3) hold and zero otherwise.
In practice not all choices might be coherent and we can treat this case by considering that the latent vector function is corrupted by a Gaussian noise with zero mean vector and covariance .^{3}^{3}3
We assume the noise variance is the same in each dimension but this can easily be relaxed.
Then we require conditions (2) and (3) to only hold probabilistically. This leads to the following likelihood function for the pair :(4) 
where denotes the th component of the vector and is the indicator function of the set . We now provide two results which allows us to simplify (4). We first compute the integral in the third product in (4).
Lemma 4.
(5)  
All proofs are in the supplementary material. We now focus on the first integral in (4), which can be simplified as follows.
Lemma 5.
(6)  
Note that eq. (6) is an expectation (with respect to ) of a product of Gaussian CDFs whose argument only depends on . We can thus write the above multidimensional integral as a sum of products of univariate integrals which can be computed efficiently, for instance by using a Gaussian quadrature rule.
Therefore, the likelihood of the choices given the latent vector function can then be written as follows.
Theorem 6.
The likelihood is
(7) 
with
(8)  
3.1.1 Likelihood for batch preference learning
In case (the latent dimension is one), we have that . This means the agent always selects a single best option. In this case, the likelihood (8) simplifies to
(9)  
The above likelihood is equal to the batch likelihood derived in [siivola2020preferential, Eq.3] and reduces to the likelihood derived in [ChuGhahramani_preference2005] when (that is the batch 2 case, i.e., ). This shows that the likelihood in (8) encompasses batch preferencebased models.
3.2 Posterior
The posterior probability of
is(10) 
where the prior over the component of is defined in (1), the likelihood is defined in (8
) and the probability of the evidence is
. The posterior is intractable because it is neither Gaussian nor Skew Gaussian Process (SkewGP) distributed. In this paper we propose an approximation schema for the posterior similar to the one proposed in [benavoli2020preferential, benavoli2021]. In [benavoli2020preferential], an analytical formulation of the posterior is available, the marginal likelihood is approximated with a lower bound and inferences are computed with an efficient rejectionfree slice sampler [gessner2019integrals]. In [benavoli2020preferential, benavoli2021]such approximation schema showed better performance in active learning and BO tasks than LP and EP. Here we do not have an analytical formulation for the posterior therefore we use a Variational (ADVI) approximation
[kucukelbir2015automatic]of the posterior to learn the hyperparameters
of the kernel and, then, for fixed hyperparameters, we compute the posterior of via elliptical slice sampling (ess) [pmlrv9murray10a].^{4}^{4}4We implemented our model in PyMC3 [salvatier2016probabilistic], which provides implementations of ADVI and ess. Details about number of iterations and tuning are reported in the supplementary.3.3 Prediction and Inferences
Let be a set including test points and . The conditional predictive distribution is Gaussian and, therefore,
(11) 
can be easily approximated as a sum using the samples from . In choice function learning, we are interested in the inference:
(12) 
which returns the posterior probability that the agent chooses the options from the set of options . Given a finite set (that is ), we can use (12) to compute the probability that a subset of is the set of nondominated options. This provides an estimate of the Paretoset based on the learned GPbased latent model.
4 Latent dimension selection
In the previous section, we provided a Bayesian model for learning a choice function. The model is conditional on the predefined latent dimension (that is, the dimension of the vector of the latent functions ). Although, it is sometimes reasonable to assume the number of criteria defining the choice function to be known (and so the dimension ), it is crucial to develop a statistical method to select . We propose a forward selection method. We start learning the model and we increase the dimension in a stepwise manner (so learning ) until some model selection criterion is optimised. Criteria like AIC and BIC are inappropriate for the proposed GPbased choice function model, since its nonparametric nature implies that the number of parameters increases also with the size of the data (as ). We propose to use instead the Pareto Smoothed Importance sampling LeaveOneOut crossvalidation (PSISLOO) [vehtari2017practical]. Exact crossvalidation requires refitting the model with different training sets. Instead, PSISloo can be computed easily using the samples from the posterior.
We define the Bayesian LOO estimate of outofsample predictive fit for the model in (10):
(13) 
where , ,
(14) 
As derived in [gelfand1992model], we can evaluate (14) using the samples from the full posterior, that is for .^{5}^{5}5We compute these samples using elliptical slice sampling. We first define the importance weights:
and then approximate (14) as:
(15) 
It can be noticed that (15) is a function of only, which can easily be computed from the posterior samples. Unfortunately, a direct use of (15) induces instability because the importance weights can have high variance. To address this issue, [vehtari2017practical] applies a simple smoothing procedure to the importance weights using a Pareto distribution (see [vehtari2017practical] for details). In Section 6.3, we will show that the proposed PSISLOObased forward procedure works well in practice.
5 Choicebased Bayesian Optimisation
In the previous sections, we have introduced a GPbased model to learn latent choice functions from choice data. We will now focus on the acquisition component of Bayesian optimization. In choicebased BO, we never observe the actual values of the functions. The data is , where is the set of the training inputs (options), is a subset of and is the choiceset for the given options . We denote the Paretoset estimated using the GPbased model as . In choicebased BO, the objective is to seek a new input point . Since can only be queried via a choice function, this is obtained by optimizing w.r.t. an acquisition function , where is the current (estimated) Paretoset. We define the acquisition function , with the aim to find a point that dominates the points in . That is, given the set of options , we aim to find such that .
The acquisition function must also consider the tradeoff between exploration and exploitation. Therefore, we propose an acquisition function which is equal to the % (in the experiments we use ) Upper Credible Bound (UCB) of with , and .
Note that the requirement for our acquisition function is strong. We could also define with different objectives in mind. For example we could seek to find a point which allows to reject at least one option in . We opted for UCB over the probability because it leads to a fast to evaluate acquisition function. In particular we only need to compute one probability for each new function evaluation. In future work we will study alternative approaches and the tradeoff between more costly acquisition function evaluations and faster convergence.
After computing the maximum of the the acquisition function, denoted with , consistently with the definition of the acquisition function, we should query the agent to express their choice among the set of options in . However, can be a very large set and human cognitive ability cannot efficiently compare more than five options. Therefore, by using the GPbased latent model, we select four options in which have the highest probability of being dominated by and query the agent on a five options set.^{6}^{6}6Details about the procedure we use to select these 4 options are reported in the supplementary material.
6 Numerical experiments
First, we assume (that is we assume that the latent dimension is known) and evaluate the performance of our algorithm on (1) the tasks of learning choice functions; (2) the use of choice functions in multiobjective BO. Second, we evaluate the latent dimension selection procedure discussed in Section 4 on simulated and real datasets.
6.1 Choice functions learning
Toy experiment
We develop a simple toy experiment as a controlled setting for our initial assessment. We consider the bidimensional vector function with .
We use to define a choice function. For instance, consider the set of options , given that
we have that and . In fact, one can notice that dominates on both the objectives, and and are incomparable. We sample inputs at random in and, using the above approach, we generate

random subsets of the 200 points each one of size (respectively ) and computed the corresponding choice pairs based on ;

random subsets each one of size (respectively ) and computed the corresponding choice pairs based on ;
for a total of four different datasets. Fixing the latent dimension , we then compute the posterior means and credible intervals of the latent functions learned using the model introduced in Section 3.2. The four posterior plots are shown in Figure 3. By comparing the 1st with the 3rd plot and the 2nd with the 4th plot, it can be noticed how the posterior means become more accurate (and the credible interval smaller) at the increase of the size dataset (from N=50 to N=150 choicesets). By comparing the 1st with the 2nd plot and the 3rd with the 4th plot, it is evident that estimating the latent function becomes more complex at the increase of . The reason is not difficult to understand. Given , includes the set of rejected options. These are options that are dominated by (at least) one of the options in , but we do not know which one(s). This uncertainty increases with the size of , which makes the estimation problem more difficult.
Considering the same inputs generated at random in , we add Gaussian noise to with and generate two new training datasets with and (i) ; (ii) . We aim to compute the predictive accuracy of the GPbased model for inferring from the set of options on additional unseen pairs . We denote the model learned using the and dataset respectively by ChoiceGP100 and ChoiceGP300. We compare their accuracy with that of an independent GP regression model which has direct access to : that is we modelled each dimension of as an independent GP and compute the posterior of the components of using GP regression from data pairs and, respectively, . We refer to this model as OracleGP: “Oracle” because it can directly query . The accuracy is:
ChoiceGP100  ChoiceGP300  OracleGP 
0.54  0.72  0.77 
averaged over repetitions of the above data generation process. It can be noticed that the increase of , the accuracy of ChoiceGP gets closer to that of OracleGP. This confirms the goodness of the learning framework developed in this paper.
Realdatasets
We now focus on four benchmark datasets for multioutput regression problems. Table 1 displays the characteristics of the considered datasets.
Dataset  #Instances  #Attributes  #Outputs 

enb  768  6  2 
jura  359  6  2 
realestate  414  5  2 
slump  103  7  2 
More details on the used datasets are in the supplementary material. By using 5fold crossvalidation, we divide the dataset in training and testing pairs. The target values in the training set are used to generate choice functions based pairs with and (i) ; (ii) . From the test dataset, we generated pairs. As before we denote the model learned using the and dataset respectively by ChoiceGP1 and ChoiceGP3 and compare their accuracy against that of OracleGP (learned on the training dataset by independent GP regression). The accuracy is:
ChoiceGP100  ChoiceGP300  OracleGP  

enb  0.74  0.77  0.77 
jura  0.44  0.47  0.53 
realestate  0.50  0.60  0.64 
slump  0.26  0.39  0.45 
As before, it can be noticed that the increase of , the accuracy of ChoiceGP gets closer to that of OracleGP.
6.2 Bayesian Optimisation
We have considered for five standard multiobjective benchmark functions: BraninCurrin (, ), ZDT1 (, ), ZDT2 (, ), DTLZ1 (, ), Kursawe (, ) and VehicleSafety^{7}^{7}7The problem of determining the thickness of five reinforced components of a vehicle’s frontal frame [yang2005metamodeling]. This problem was previously considered as benchmark in [daulton2020differentiable]. (, ). These are minimization problems, which we converted into maximizations
so that the acquisition function in Section 5 is welldefined. We compare the ChoiceGP BO (with ) approach proposed in this paper against ParEGO.^{8}^{8}8We use the BoTorch implementation [balandat2020botorch].
For ParEGO, we assume the algorithm can query directly and, therefore, we refer to it as OracleParEGO.^{9}^{9}9The most recent MO BO approaches mentioned in Section 1 outperform ParEGO. We use ParEGO only as an Oracle reference.
Conversely, ChoiceGP BO can only query via choice functions. We select and use UCB as acquisition function for ChoiceGP BO. We also consider a quasirandom baseline that selects candidates from a Sobol sequence denoted as ‘Sobol‘. We evaluate optimization performance on the five benchmark problems in
terms of loghypervolume difference, which is defined as the difference between the hypervolume of
the true Pareto front^{10}^{10}10This known for the six benchmarks. and the hypervolume of the approximate Pareto front based on the observed data .
Each experiment starts with initial (randomly selected) input points which are used to initialise OracleParEGO. We generate pairs of size by randomly selecting subsets of these points. These choices are used to initialise
ChoiceGP BO. A total budget of iterations are run for both the algorithms. Further, each experiment is repeated 15 times with different initialization. In these experiments we optimize the kernel hyperparameters
by maximising the marginal likelihood for OracleParEGO and its variational approximation for ChoiceGP.
Figure 4 reports the performance of the
three methods. Focusing on BraninCurrin, DTLZ1, Kursawe, and VehicleSafety, it can be noticed how ChoiceGP BO convergences to the performance of the OracleParEGO at the increase of the number of iterations. The convergence is clearly slower because ChoiceGP BO uses qualitative data (choice functions) while OracleParEGO uses quantitative data (it has directly access to ). However, the overall performance shows that the proposed approach is very effective. In DLTZ1 and DLTZ2, ChoiceGP BO outperforms OracleParEGO. The bad performance of OracleParEGO is due to the used acquisition function, which does not correctly balance exploitationexploration in these two benchmarks. Instead, the UCB acquisition function for ChoiceGP BO works well in all the benchmarks.
Computational complexity: The simulations were performed in a standard laptop. On average, the time to learn the surrogate model and optimise the acquisition function goes from 30s () to 180s ().
6.3 Unknown latent dimension
We assume that is unknown and evaluate the latent dimension selection procedure proposed in Section 4. First, we consider the onedimensional and so . We generate training datasets with and sizes and, respectively, and 10 test datasets with size . The following table reports the PSISLOO (averaged over the five repetitions) and the average accuracy on the test set for four ChoiceGP models with latent dimension .
N=30  N=300  

PSISLOO  acc. test  PSISLOO  acc. test  
1 (10/10)  10  0.75  75  0.93 
2  35  0.64  165  0.91 
3  44  0.64  333  0.86 
4  69  0.62  388  0.84 
By using PSISLOO (computed on the training dataset) as latent dimension selection criterion, we were able to correctly select the true latent dimension in all the repetitions (10 out of 10). The selected model has also the highest accuracy in the test set. We now focus on the bidimensional vector function and consider three different sizes for the training dataset .
N=30  N=50  N=300  

PSISLOO  acc. test  PSISLOO  acc. test  PSISLOO  acc. test  
1  56  0.20  89  0.23  493  0.30 
2  39  0.32  47  0.51  236  0.72 
3  39  0.32  49  0.49  269  0.65 
4  42  0.30  53  0.43  277  0.64 
For , has the highest PSISLOO in 4/10 cases, in 4/10 cases and in 2/10 cases. For , has the highest PSISLOO in 6/10 cases and in 4/10 cases. For , the PSISLOO selects in 10/10 cases. Note that, is never selected. Considering that the models are nested , this shows that the selection procedure works well even with small datasets by never selecting a latent dimension that is smaller than the actual one. Moreover, PSISLOO is able to select the correct dimension () at the increase of . In the supplementary material, we have reported a similar analysis for the datasets in Table 1 confirming that the procedure also works for realdatasets.
7 Conclusions
We have developed a Bayesian method to learn choice functions from data and applied to choice function based Bayesian Optimisation (BO). As future work, we plan to develop strategies to speed up the learning process by exploring more efficient ways to express the likelihood. We also intend to explore different acquisition functions for choice function BO.