Observation Models, Surrogates, and Inference

Every preferential Bayesian optimization (PBO) system rests on three modeling choices, which Part IV introduced one at a time. The observation model says how a comparison, a ranking, or a choice arises from the latent utility: in Chapter 16, a probit or logistic link applied to a utility difference. The surrogate is the prior over the utility: in Chapter 18, a Gaussian process. Posterior inference combines the two and hands the result to the acquisition function: in Chapter 17, usually the Laplace approximation.

From 2017 to September 2026 the default combination barely changed. A Gaussian process prior, a probit link, and the Laplace approximation, which is the model of Chu and Ghahramani (2005), remain what BoTorch's PairwiseGP implements (Meta Platforms, Inc., 2026h). What changed is the scrutiny each layer received, and this chapter follows the three layers in order, ending with what the default implementation decides for you.

27.1 The baseline models #

Chu and Ghahramani (2005) placed a Gaussian process prior on the latent utility $f$ and used the likelihood

which is Thurstone's comparative judgment with Gaussian noise of scale $\sigma$ on each option (Section 16.3). They fitted the posterior with the Laplace approximation and chose hyperparameters by the approximate evidence. The pairwise model of Brochu et al. (2007) is the same Thurstone model, as Koyama et al. (2020) point out. Houlsby et al. (2011) took another route to the same place: a preference kernel that turns pairwise preference learning into Gaussian process classification over pairs.

The paper that named the field used neither the probit link nor, as far as its text says, the Laplace approximation. González et al. (2017) modeled the probability that $\vx$ wins a duel against $\vx'$ with a Gaussian process classifier on the dueling space $\X \times \X$, with a Bernoulli likelihood and the latent joint reward squashed by a logistic function, assuming the latent function has the form $f([\vx, \vx']) = g(\vx') - g(\vx)$; their text says only that the posterior is intractable and needs approximations, without naming one. They listed the doubled input dimension as a limitation, and Nguyen et al. (2021) later criticized the model for working in twice the dimension of the objective, for not modeling the objective directly, and for not modeling ties. Because the skew Gaussian process theorem of Section 27.4 concerns the probit likelihood, it does not apply directly to this original model (inference).

Since then the two links have been used side by side. BoTorch offers both: PairwiseProbitLikelihood with $\Prob(v \succ u) = \Phi((f(v) - f(u))/\sqrt{2})$ and PairwiseLogitLikelihood with $\operatorname{sigmoid}(f(v) - f(u))$ (Meta Platforms, Inc., 2026g). The logistic link appears in the noise analysis of qEUBO (Astudillo et al., 2023), in POP-BO (Xu et al., 2024b), and in MR-LPF (Kayal et al., 2025). For two options, a softmax over utilities is exactly the logistic Bradley-Terry model. Both links are random utility models (Section 16.5): the probit one comes from Gaussian noise on each option, and the logistic one from Gumbel noise, a correspondence Kayal et al. use when they discuss lower bounds.

27.2 Extending the observation model #

A person's answer to "which is better?" can carry more, or different, information than one bit. Since 2017, researchers have extended the likelihood in eight directions. Table 27.1 summarizes them; the paragraphs after it give the details that matter for judging each one.

Choices from sets, and rankings. The sequential line search of Koyama et al. (2017) records the person's pick along a slider as a preference for the chosen point over both the current best and the point of highest expected improvement, using the Bradley-Terry-Luce model, the extension of Bradley-Terry to a choice among $m$ options; the Sequential Gallery models pairwise data with the Thurstone model and choices among $m$ options with the Bradley-Terry-Luce form (Koyama et al., 2020). Siivola et al. (2021) give a likelihood for arbitrary parallel feedback on two or more points, focus on the batch winner, and argue that full rankings are laborious for people. Nguyen et al. (2021) build a surrogate inspired by the multinomial logit and its ranking extension, interpretable as Gaussian process regression with i.i.d. Gumbel noise, which handles top-$k$ rankings and ties. Benavoli et al. (2023) handle set-valued choices (picking, say, three of five options shown, when none of the three is clearly better than the others) with several latent utilities; their abstract reports only simulations.

Ties and indifference. Bıyık et al. (2019) added an "About Equal" answer with a minimum perceivable difference $\delta \ge 0$: $\Prob(\text{About Equal}) = (e^{2\delta} - 1)\,\Prob(\text{choose }a)\,\Prob(\text{choose }b)$, which reduces to the softmax model at $\delta = 0$. Their model is a linear reward, not a Gaussian process, and comes with simulations and a user study. C-GLISp encodes better, worse, or similar as constraints on a radial basis function surrogate (Zhu et al., 2022). The tutorial of Benavoli and Azzimonti lists nine likelihoods, among them a just-noticeable-difference model that adds statements of indistinguishability, $\lvert u(\vx) - u(\vx')\rvert/\delta \le 1$, going back to Luce's 1956 notion of a discrimination threshold. Erarslan et al. (2025) represent indifference with a just-noticeable-difference threshold in an extended Thurstone model and report clear gains when 10% to 20% of comparisons are indifferent. The field has converged on two treatments, a threshold parameter in a Thurstone or logit model, or ties in a multinomial logit (inference). We found no Gaussian process PBO paper with an "abstain" or "skip" outcome distinct from a tie, the answer a person gives when they cannot compare the options at all.

Ordinal labels, strength, and confidence. ROIAL (Li et al., 2021) combines pairwise preferences with ordinal labels (very bad, bad, neutral, good), arguing that an $r$-level ordinal query yields at most $\log_2 r$ bits against one bit for a preference; its participants were 3 non-disabled people tuning 4 exoskeleton gait parameters. Wu et al. (2025a) add to the probit preference likelihood a Likert confidence likelihood with learnable cut points and a lapse rate; on human comparisons of robot gaits, they report that models trained with the confidence ratings "consistently achieve lower Brier scores and higher F1 scores", that is, smaller squared errors of the predicted probabilities and more accurate predicted choices.

Response times. How long a person takes to answer says something about how different the options felt. Shvartsman et al. (2024) approximated the drift-diffusion likelihood with a moment-matched skewed three-parameter distribution so that response times can enter a variational Gaussian process, and describe this as the first Gaussian process model with a joint likelihood of drift-diffusion response times and choices. On one participant's judgments of 1,225 pairs of quadruped robot gaits and on a visual psychophysics data set, the joint drift-diffusion models predicted held-out choices significantly better than the choice-only model, especially with small training sets, while a simpler stacking approach, which gives the choice model the logarithm of the response time as an extra input, did not improve on choice alone; on a recommender-system data set every joint model improved, stacking most. Li et al. (2024a) used the EZ-diffusion model for fixed-budget best-arm identification with linear utilities and proved that, for strong preferences, response times add information about preference strength.

Failure, validity, and crashes. Some experiments produce no output at all. Validity labels and feasibility surrogates attach a second label to each sample, and CrashPBO (Menn et al., 2026b) treats a crash report as a second kind of outcome, reduces crashes by 63% on synthetic benchmarks, and was validated on three robot platforms.

Intransitivity. Chau et al. (2022) placed a Gaussian process on a skew-symmetric preference function $g(\vx, \vx')$ with a kernel they call the generalized preferential kernel, which can represent cycles. On real data the new model's accuracy against the Chu-Ghahramani model was 0.78 against 0.51 for chameleon contests, 0.83 against 0.80 for flat lizards, 0.59 against 0.51 for NFL games from 2000 to 2018 (where linear pairwise logistic regression did best, at 0.65), and 0.74 against 0.66 for an arXiv citation graph. The authors conclude that their findings support the conjecture that violations of rankability, the assumption that one utility orders all items, are ubiquitous in real preference data. These data sets are ones where intransitivity is expected (contests and sports), the linear model won on the only human competition data, and none is a single user's design preference, so carrying the conclusion over to typical PBO settings lacks support (inference).

Heteroscedastic noise. Some comparisons are harder than others. Sinaga et al. (2026) ask the user for reliable "anchor" points and turn them into an input-dependent noise map with kernel density estimation, used in the probit likelihood $\Phi\big((f(\vx) - f(\vx'))/\sqrt{\sigma^2(\vx) + \sigma^2(\vx')}\big)$ or a logistic likelihood with input-dependent temperature; their consistency analysis holds only under an idealized i.i.d. anchor model, and the evaluation uses the Sushi and Candy data without live users. Fauvel and Chalk (2021) argue that acquisition functions for binary and preferential Bayesian optimization must separate epistemic from aleatoric uncertainty, since only the former should drive exploration.

Several users, and drift. Models for many users share low-dimensional structure across them, and crowdGPPL (Simpson and Gurevych, 2020) scales a combination of matrix factorization and Gaussian processes to thousands of users and items; Section 20.5 describes these models. We found no PBO or Gaussian process preference paper that models a utility changing over time; the finite-arm theory that does exist is in Section 29.10.

What the human evidence supports. For ties there is one user study with a linear model; for ordinal labels, three participants; for confidence, one offline robot-gait data set and two small user studies; for response times, three data sets, few participants, and offline fits; for crashes, three robot platforms; for heteroscedastic noise, two benchmarks converted from rating data; for intransitivity, no design preferences at all. Most evaluations compare model fits offline rather than closing the loop with people (inference). Of the signals that cost the person nothing extra, response times and confidence are the two that have improved models on real human data, which makes them the best-supported upgrades; recording both in every session costs nothing and lets them be tested in closed-loop PBO later (inference).

27.3 Surrogates beyond the Gaussian process #

The Gaussian process utility remains the workhorse, but four families of alternatives appeared, each for a specific reason.

Gaussian process variants. Bıyık et al. (2020) fixed $f(\bar{\Psi}) = 0$ at a reference point through the kernel, removing the additive non-identifiability of Section 18.4 in the prior. The skew Gaussian process line was built by Benavoli, Azzimonti, and Piga: a 2020 classification paper introduced the skew Gaussian process and proved it conjugate to the probit likelihood (Benavoli et al., 2020), a 2021 paper extended conjugacy to normal and affine probit likelihoods and their products (Benavoli et al., 2021a), and Benavoli and Azzimonti (2024) showed that a Gaussian process with linear inequality constraints is a skew Gaussian process, using it for monotone preference learning.

Radial basis functions, neural networks, and trees. GLISp (Bemporad and Piga, 2021) fits a radial basis function surrogate by linear or quadratic programming so that it satisfies the observed preferences as far as possible, and reports that, with the same number of comparisons, it usually gets closer to the optimum than PBO at lower computational cost; GLISp-r (Previtali et al., 2023) adds a proof of global convergence. These surrogates are not probabilistic, so the acquisition functions of Chapter 19 do not apply to them. Neural dueling bandits (Verma et al., 2025) estimate rewards from preference feedback with a neural network and prove sublinear regret, on synthetic data, and Wang et al. (2025a) replace the Gaussian process utility with an ensemble of monotone neural networks in preference exploration. DT-PBO (Leenders et al., 2025) builds shallow decision trees directly from comparisons, with a Laplace approximation in each leaf, and its authors report convergence competitive with Gaussian process PBO on 8 benchmark functions, particularly on rugged ones.

Amortized surrogates. PABBO (Zhang et al., 2025a) starts from the observation that the non-conjugate likelihood makes every PBO step expensive. It meta-learns surrogate and acquisition together in a transformer neural process, trained with reinforcement learning plus an auxiliary preference-prediction loss on synthetic tasks such as Gaussian process samples. Its abstract claims it is "several orders of magnitude faster than the usual Gaussian process-based strategies"; the paper reports a 12-fold average speedup over the fastest Gaussian process strategy on five problems and about 10-fold over noisy batch expected improvement on Gaussian process tasks. Its default evaluation configuration uses noise-free comparisons (Zhang, 2025). The authors list as limitations that the query set grows with dimension, each dimension needs its own model, how to build the pretraining data is open, and the person making comparisons is not modeled.

Priors across users, and generative latent spaces. Cross-user structure has been implemented through collaborative models such as crowdGPPL; through transfer at the acquisition level, as in Meta-PO (Li et al., 2025a), which stores earlier users' sessions as separate Gaussian processes and transfers them with a rank-weighted acquisition function (Section 32.5 reports its user study); and through guidance models trained on population data, such as the deep encoder of Granley et al. (2023). Several systems also optimize in the latent space of a generative model, such as FontCraft (Tatsukawa et al., 2025) in a font-style space and GimmBO (Liu et al., 2026b) over the merging weights of 20 to 30 diffusion adapters; none of these papers models how the geometry of the latent space relates to the utility's lengthscale (inference).

Overall, the debate about surrogates since 2020 has been mainly about the quality of inference in the same model, Gaussian approximation against the exact skew Gaussian process, not about the model class. Changes of model class were driven by structure (monotonicity, intransitivity, several utilities) or speed (PABBO, GLISp), not by evidence that a Gaussian process utility prior is wrong for human preferences (inference).

27.4 Posterior inference #

Chapter 17 introduced the approximations, and Section 17.7 reported the main evidence on how well they do; this section adds the details behind that evidence and the choices the libraries make. Table 27.2 lists where each method is used.

The Laplace approximation is the library default, mostly for speed. Bıyık et al. (2020) wrote of expectation propagation that "while it is more accurate than Laplace approximation, it is slower in practice", and chose Laplace; projective PBO (Mikkola et al., 2020) uses it "for the sake of simplicity". Koyama et al. use maximum a posteriori estimates with a tight log-normal prior on each hyperparameter, written $\mathrm{LN}(\mu, 0.01)$ in the paper, with $\mu = 0.5$ for every lengthscale and $\mu = 0.2$ for the amplitude, which nearly fixes the lengthscale (Koyama et al., 2020). Expectation propagation is the inference method in Siivola et al.'s experiments on real data, and the preferential Gaussian process sampler of optuna-dashboard uses it to fit its hyperparameters (Optuna developers, 2026c).

The exact posterior is a skew Gaussian process. Benavoli et al. (2021c) proved this in their PBO paper (Theorem 29.1) and sampled the posterior with LinESS, rejection-free elliptical slice sampling for linearly truncated Gaussians, at a cost they give as $O(n^3)$ time and $O(n^2)$ memory, the same bottleneck (a Cholesky factorization) as an ordinary Gaussian process; they report consistently better convergence and computation time than Laplace-based PBO.

How wrong the Gaussian approximations are. Takeno et al. (2023) took Gibbs sampling as ground truth (10,000 samples after 1,000 burn-in, thinned by 10), with an RBF kernel, noise variance $10^{-4}$, and uniformly random duels. Beyond the findings of Section 17.7, that the Laplace mode can be very far from the mean at small noise and that expectation propagation is accurate for means and credible intervals but distorts duel probabilities, they report that on the Ackley function expectation propagation underestimates probabilities near 0 or 1 and does not preserve the true ordering, which affects acquisition functions that depend on the joint distribution of two points, and that Gibbs sampling mixed faster than LinESS. They consider unrealistic the setting in which Benavoli et al. compared expectation propagation with the skew Gaussian process, where all inputs in one interval lose and all inputs in another win; an earlier study of binary Gaussian process classification had found expectation propagation accurate (Kuss and Rasmussen, 2005). And they reject Laplace for the posterior yet use the Laplace evidence for hyperparameters, so the accuracy needed for hyperparameters and for acquisition is being treated separately.

Takeno et al.'s strongest evidence comes from noise variance $10^{-4}$, where the posterior is nearly a truncated Gaussian; with noisier comparisons, closer to people, the skew should weaken. As of September 2026, no paper measures how large the errors of Laplace and expectation propagation are at realistic human noise levels, or replicates Takeno et al.'s calibration on real human comparisons. The experiment is cheap to describe: take comparisons from real people, run exact sampling as ground truth, and measure how far each approximation's predicted duel probabilities fall from it.

27.5 Inside the default implementation #

Most people who run PBO run BoTorch's PairwiseGP, so its defaults are, in effect, the field's defaults. Table 18.2 listed the choices inside version 0.18.1 (Meta Platforms, Inc., 2026h; Meta Platforms, Inc., 2026g). Three of them shape what the model can learn. The probit likelihood is $\Phi((f(v) - f(u))/\sqrt{2})$, which fixes the noise at 1, so the output scale (BoTorch's name for the kernel's amplitude $\sigma_f^2$) acts as the inverse of the noise level: it and the lengthscale must be learned jointly from comparisons alone, with no separate noise parameter, under a smoothed box prior on $[0.01, 100]$ and a constraint of $[0.005, 200]$ that a source comment calls a rule of thumb to keep estimates away from scales that saturate $\Phi$. The probit argument is clipped to $[-3, 3]$, which caps the likelihood of any single comparison at about 0.9987 during fitting, an implicit robustness to mislabeled answers (Section 18.6). And the lengthscale prior is Gamma(2.4, 2.7), initialized at its mode, in every dimension.

The changelog adds the history (Meta Platforms, Inc., 2026e). Version 0.9.0 (August 2023) added a pairwise Bayesian active learning by disagreement acquisition function; 0.10.0 (February 2024) added qEUBO; 0.12.0 (September 2024) switched most models to lengthscale priors that scale with dimension but explicitly excluded PairwiseGP; 0.18.0 (June 2026) fixed PairwiseGP's state handling in evaluation and cross-validation; and in 0.18.1 (June 2026) PairwiseGP still uses Gamma(2.4, 2.7). The dimension-scaled priors came from a study of scalar Bayesian optimization by Hvarfner et al. (2024), the subject of Chapter 30. The preferential sampler in optuna-dashboard, based on Takeno et al., uses a Matérn 3/2 kernel with one lengthscale per dimension and a Gamma(5, 10) prior (mode 0.4), also independent of dimension (Optuna developers, 2026c).

What the fixed prior does in high dimension can be computed (inference). The mode of Gamma(2.4, 2.7) is $(2.4 - 1)/2.7 \approx 0.52$ in every dimension, while the dimension-scaled prior's mode on the unit cube grows from about 0.65 at 10 dimensions to about 2.05 at 100. Since random points in $[0, 1]^d$ lie about $\sqrt{d/6}$ apart (Section 30.1.2 derives this), the default preference model in high dimension works where kernel values between typical points are near zero, and every comparison is nearly uninformative about every other point. This has not been verified directly on PairwiseGP.

27.6 Comparison graphs and conditioning #

The numerical trouble reported in 2026 has a linear-algebra core that Section 18.5 derived: the Hessian $\mW$ of the negative log-likelihood is the weighted graph Laplacian of the comparison graph, whose nodes are compared inputs and whose edges are answered pairs, and it has one zero eigenvalue for every connected component. The Laplace posterior precision $\mK^{-1} + \mW$ is always full rank, because the prior contributes $\mK^{-1}$, so the practical problem is not singularity but poor conditioning when the prior is weak along the shift directions, that is, when lengthscales are long or the output scale is large. The same fact has a statistical face, which Figure 27.1 shows: a shift direction that no comparison constrains is a relative utility that no answer has measured, and only the prior can say anything about it.

Shao et al. (2026), a 2026 preprint, report that in the pipeline of PairwiseGP, the Laplace approximation, and EUBO, the new pairs EUBO selects share no candidate with earlier queries, so each pair forms an isolated component of the comparison graph and the likelihood Hessian becomes rank deficient. They call the deficiency structural, one that "cannot be resolved by changing the surrogate modeling approach", and note that existing remedies either force comparisons to stay connected, wasting query budget, or apply uniform regularization that perturbs well-constrained directions. Their correction, whose size Section 18.5 reports, adds to the Hessian's diagonal a term scaled by the prior uncertainty, only during model fitting; an adaptive version applies it only when the surrogate is confident, and the benchmarks include a 16-dimensional plasma-medicine controller. Pukdee et al. (2026) give conditions under which the Bradley-Terry model still recovers the conditional preference distribution when the data violate it, and show that the margin and the connectivity of the comparison graph govern sample efficiency. Lemma 1 of local PBO (Menn et al., 2026a) shows that gradient and Hessian estimates taken from the Laplace posterior inherit the bias of its mode; the authors add that the Hessian estimate can be especially sensitive to the lengthscale and to the conditioning of the kernel matrix, and their experiments bound the lengthscale to $[0.05, 0.5]$ on the unit cube.

How hyperparameters are learned in practice. The approaches we found range from a weak Gamma prior with the Laplace evidence (BoTorch) and the exact skew Gaussian process evidence (Benavoli et al. 2021) to a tight log-normal prior that nearly fixes the lengthscale (Koyama et al.) and box constraints (local PBO). Whether any of them can identify lengthscales from the tens to one or two hundred comparisons a person gives is unstudied, and Section 30.3.3 states the gap and the simulation that would close it.

27.7 Settled, contested, missing #

Further reading #

• Benavoli and Azzimonti (2026a), the tutorial by Benavoli and Azzimonti, is the most complete single treatment of Gaussian process models for preferences and choices, including nine likelihoods.
• Benavoli et al. (2021c) and Takeno et al. (2023) are the two sides of the inference debate; read them together.
• The source of PairwiseGP (Meta Platforms, Inc., 2026h) is short and worth reading before trusting its defaults.
• Shvartsman et al. (2024) is the clearest example of adding a free behavioral signal, response time, to the likelihood.