When the Posterior Is Not Gaussian

Gaussian process regression owes its tidy formulas to one fact: a Gaussian prior combined with Gaussian observation noise gives a Gaussian posterior. The comparison models of Chapter 16 break that fact. Their likelihood is a probit or logistic curve, not a Gaussian, and the posterior over the utility is no longer Gaussian. Everything downstream, the predictive mean and band, the probability of the next answer, and the acquisition functions of Chapter 19, needs that posterior.

This chapter works almost entirely with the smallest case that shows the difficulty: one utility difference, a Gaussian prior, and a few comparisons. Its exact posterior can be computed on a grid and drawn, so each approximation can be checked against the truth. We then move to two and more latent values, where the exact answer turns out to be skewed only in the directions that comparisons touch, and end with the evidence on how much the choice of approximation matters in practice.

17.1 Non-Gaussian likelihoods #

Recall why regression was easy. With a prior $\vf \sim \N(\mathbf{0}, \mK)$ and observations $\vy = \vf + \bm{\varepsilon}$ with Gaussian noise, Bayes' rule multiplies two Gaussian densities in $\vf$. The product of Gaussian densities is again a Gaussian density up to a constant (Section 4.6), so the posterior has a closed form and so does its normalizing constant, the marginal likelihood. Prior and likelihood form a conjugate pair.

A comparison does not. Write $\vf$ for the utilities at the inputs that have been compared, and suppose $m$ answers have been recorded, the $k$-th saying that input $v_k$ was preferred to input $u_k$. With the probit model of Equation (16.2), Bayes' rule gives

The prior is Gaussian, but each likelihood factor is an S-shaped curve in one direction of $\vf$, the direction of the difference $f_{v_k} - f_{u_k}$. The product is not Gaussian, and $Z$ is the probability that a Gaussian vector lands in a region bounded by $m$ soft walls, an $m$-dimensional integral with no closed form in general (Section 5.7). Gaussian process classification, in which each input carries a yes-or-no label $y_i = \pm 1$ instead of a number, has the same structure, with one factor $\Phi(y_i f(\vx_i))$ per labeled input (Rasmussen and Williams, 2006, ch. 3), so the methods of this chapter were developed and tested there first.

17.1.1 The smallest case #

Strip the problem down to one number. Let $\Delta = g(A) - g(B)$ be the utility difference between two options, as in Section 16.6, with a Gaussian prior $\Delta \sim \N(0, v_0)$. Here $v_0$ is a variance; Equation (16.6) described the same kind of belief by its standard deviation $v$, so $v_0 = v^2$. Suppose a person says once that $A$ is better. Write $s = \sqrt{2}\,\sigma$ for the noise on the difference. The posterior is

The factor 2 is $1/Z$: by Equation (16.6) with a belief centered at zero, the prior probability of the answer is $\Phi(0) = 1/2$. A Gaussian density multiplied by a normal distribution function is called a skew-normal density, a family introduced by Azzalini (1985). It keeps the upper tail of the prior and cuts away the lower one, softly when the noise is large and sharply when it is small. Its mean has a closed form.

With $v_0 = 1$, the mean grows from about 0.46 at $\sigma = 1$ to $\sqrt{2/\pi} \approx 0.80$ as the noise vanishes, and the standard deviation falls from about 0.89 to 0.60 (Exercise 17.1). The mode, the peak of the density, behaves differently: it solves $\Delta/v_0 = \phi(\Delta/s)/(s\,\Phi(\Delta/s))$, and as $s$ shrinks it slides toward zero, to 0.05 at $\sigma = 0.01$. In the noise-free limit the posterior is the prior's upper half, a half-normal, whose peak sits at the cut while its mass lies well to the right of it.

Some things to try:

• Shrink the noise. Drag $\sigma$ toward 0.01. The likelihood becomes a step, the posterior becomes the half-normal, and the mode and the mean separate: the mode goes to the cut, the mean stays near 0.80.
• Raise the noise. At $\sigma = 2$ the likelihood is a gentle slope, the posterior is a slightly shifted bell, and mode and mean nearly coincide.
• Let both win. With three wins for $A$ and three for $B$ the posterior is pinned from both sides and is close to a Gaussian centered at zero. Skew comes from one-sided evidence.
• Let $A$ keep winning. Each further win pushes the mass up, but with small noise the left edge stays sharp: after a run of identical answers the evidence says "$A$ is better" more firmly, not "by how much".

This lopsided shape is what every method below has to summarize, usually with a single Gaussian.

17.2 The Laplace approximation #

The simplest summary puts a Gaussian at the posterior's peak and gives it the peak's curvature. A bell shape is determined by where it is highest and how fast it falls off from there; if the posterior is close to a bell, those two facts recover it. This is the Laplace approximation, named after Laplace's method for approximating integrals, and Tierney and Kadane (1986) showed how far it goes for Bayesian computation: approximate posterior moments and marginal densities need only a maximization and the curvature at the maximum.

Write $\Psi(\vf) = \log p(\D \mid \vf) + \log \N(\vf; \mathbf{0}, \mK)$ for the log of the unnormalized posterior, and $\hat\vf$ for its maximizer, the mode or maximum a posteriori estimate.

Two things are needed: the mode and the curvature there. For the probit likelihood the log likelihood is concave, because $\log\Phi$ is concave, and the log prior is a concave quadratic, so $\Psi$ has a single maximum and Newton's method finds it quickly. Newton's method replaces the function by its quadratic approximation at the current point and jumps to that quadratic's maximum.

Starting from $\vf = \mathbf{0}$ and repeating Equation (17.4), halving the step whenever $\Psi$ fails to increase, converges in a handful of iterations. For classification, where $\mW$ is diagonal, Rasmussen and Williams (2006) give a numerically stable version as their Algorithm 3.1. For comparisons $\mW$ is not diagonal; Section 18.2 works out its structure.

The same expansion approximates the normalizing constant $Z$, the marginal likelihood used to fit hyperparameters (Section 9.3). Integrating the Gaussian of step 4 gives

the Laplace evidence that Chu and Ghahramani (2005) used for preferences (their equation 12) and that BoTorch maximizes to fit its preference model (Section 18.6).

17.2.1 Where it goes wrong #

The figure below adds the Laplace Gaussian to the exact posterior.

At moderate noise the blue curve sits close to the exact one. As the noise shrinks it fails in a specific way. At $\sigma = 0.01$ the exact posterior has mean 0.80, standard deviation 0.60, and gives $B$ a probability of 0.004 of being better. The Laplace Gaussian is centered at the mode, 0.05, with standard deviation 0.27, and gives $B$ a probability of 0.43: after a single nearly noise-free answer it is almost as unsure of the order as before (Exercise 17.2). The mode sits at the edge of the mass, and a Gaussian centered there must spill over the edge.

The general lesson was established for classification. Kuss and Rasmussen (2005) compared Laplace's method and expectation propagation with long sampling runs on binary Gaussian process classifiers and found that Laplace's method "systematically underestimates the mean", so that the approximate posterior over latent functions has "too small amplitude" and its predictive probabilities are over-conservative, although the sign of the latent function is mostly right; they concluded that it is "so inaccurate that we advise against its use, especially when predictive probabilities are to be taken seriously." The same pattern appears for comparisons: with nearly noise-free duels, the mode can be very far from the mean (Takeno et al., 2023).

Laplace's method survives because it is fast and simple. In robot reward learning from pairwise preferences, Bıyık et al. (2020) described expectation propagation as "more accurate than Laplace approximation" but "slower in practice", and chose Laplace for its computational efficiency. It is also the default in BoTorch (Section 18.6).

17.3 Expectation propagation #

The Laplace approximation looks at one point of the posterior. A better Gaussian would match the posterior's mean and variance, its mass rather than its peak. Computing those moments for the full posterior is as hard as the original problem, but computing them for a posterior with only one non-Gaussian factor is easy. Expectation propagation (EP) builds the approximation from such one-factor problems (Minka, 2001).

EP replaces each likelihood factor $t_k(\vf)$, called a site, by an unnormalized Gaussian $\tilde t_k(\vf)$ in the same direction, so that the approximate posterior $q(\vf) \propto \N(\vf; \mathbf{0}, \mK) \prod_k \tilde t_k(\vf)$ is Gaussian. It then refines one site at a time.

The tilted distribution in step 3 has one non-Gaussian factor, which acts along one direction, so its moments reduce to a one-dimensional problem. For a probit site the answer is in closed form. Let the cavity give the difference $\Delta$ the Gaussian $\N(\mu, v)$, with $v$ a variance like $v_0$, and let the site be $\Phi(y\, \Delta / s)$ with $y = +1$ if $A$ won and $-1$ if $B$ won.

These are the same quantities, with $s = 1$, as in the expectation propagation algorithm for probit classification (Rasmussen and Williams, 2006, ch. 3). Each site update is a rank-one change to the covariance, which costs $O(n^2)$ for $n$ latent values, so a sweep over $m$ sites costs $O(mn^2)$, comparable to a Newton iteration when $m$ is of the order of $n$.

With a single comparison there is a single site, the tilted distribution is the exact posterior, and EP returns its exact mean and variance: at $\sigma = 0.03$, mean 0.80 and standard deviation 0.60. That is far better than Laplace's 0.12 and 0.32. Matching moments has a cost of its own, though. A Gaussian with the right mean and variance still puts mass where the exact posterior has none: EP gives $B$ a probability of 0.093 of being better, against an exact value of 0.013. With five wins at $\sigma = 0.1$, several sites act on the same direction and EP is no longer exact (mean 0.93 and standard deviation 0.45, against 0.90 and 0.58).

That small misplacement matters for comparisons specifically. Takeno et al. (2023) took a long run of a sampling method as ground truth (10,000 samples from the Gibbs sampler of Section 17.5, after discarding the first 1,000 and keeping every tenth), and found that EP estimates means and credible intervals very accurately but, for a duel already observed with $\vx_w$ beating $\vx_l$, overestimates the probability that $f(\vx_w) \le f(\vx_l)$. On the Ackley function, a standard test function, it underestimated duel probabilities near 0 or 1, and its estimates did not keep the true ordering between pairs, which can change which pair an acquisition function picks.

For classification the verdict on EP is favorable. Kuss and Rasmussen found its predictive probabilities and marginal likelihood estimates very close to those of long sampling runs (Kuss and Rasmussen, 2005), and a broader comparison of approximations for binary Gaussian process classification concluded that "the Expectation Propagation algorithm is almost always the method of choice unless the computational budget is very tight" (Nickisch and Rasmussen, 2008). Unlike Newton's method on a concave function, EP is not guaranteed to improve an objective at every step, so implementations damp the site updates when they oscillate.

17.4 Variational inference #

A third route turns approximation into optimization. Pick a family of simple distributions, Gaussians here, and find the member closest to the posterior, measuring closeness by the Kullback-Leibler divergence $\KL(q \,\|\, p)$ (Section 6.2). The divergence to the posterior cannot be computed directly, because it contains the unknown $\log Z$, but it can be minimized anyway.

$\mathcal{L}(q)$ is the evidence lower bound (ELBO). For a Gaussian $q$ and a probit likelihood, every term is cheap: the divergence between two Gaussians has a closed form, and the expected log likelihood is a sum of one-dimensional integrals, one per comparison, each over the Gaussian that $q$ assigns to a difference, so standard gradient-based optimizers can maximize it.

The direction of the divergence decides the character of the fit. $\KL(q \,\|\, p)$ averages $\log(q/p)$ under $q$, so it is enormous wherever $q$ puts mass and $p$ has almost none. The optimal $q$ therefore stays inside the posterior's support and tends to be too narrow. Expectation propagation's local moment update instead chooses the Gaussian closest to the tilted distribution in the other direction, $\KL(\hat p_k \,\|\, q)$, which punishes $q$ for missing mass and tends to make it too wide.

At $\sigma = 0.01$ the variational Gaussian has mean 0.88 and standard deviation 0.30, against the exact 0.80 and 0.60, and gives $B$ a probability of 0.002 of being better, close to the exact 0.004. It respects the hard edge that EP crosses, but it halves the uncertainty. Each method errs where its criterion is blind.

Variational inference is the approximation that scales. With inducing points, a small set of pseudo-inputs that summarize the function (Titsias, 2009), and stochastic optimization, it handles thousands of comparisons and any likelihood whose expected log can be estimated. That is why it underlies much of the recent preference modeling: crowd preference learning with thousands of users and items (Simpson and Gurevych, 2020), top-$k$ rankings (Nguyen et al., 2021), choice functions (Benavoli et al., 2023), response times (Shvartsman et al., 2024), mixed likelihoods that combine comparisons with confidence ratings, in a 2025 preprint (Wu et al., 2025a), and the experiments of qEUBO (Astudillo et al., 2023). As of September 2026 we found no systematic assessment of its error for preference likelihoods comparable to those for Laplace and EP.

17.5 Sampling #

The three methods so far replace the posterior by a Gaussian. Markov chain Monte Carlo (MCMC) methods instead draw a sequence of samples whose distribution converges to the posterior itself. Any quantity of interest, the mean, a credible interval, the probability that one option beats another, is then estimated by averaging over the samples. The answer becomes exact as the number of samples grows, at the price of computation and of having to judge when a chain has run long enough.

For a Gaussian prior times a likelihood, elliptical slice sampling is the natural choice. Murray et al. (2010) designed it for "models with multivariate Gaussian priors": it has "simple, generic code", "no free parameters" to tune, and "works well for a variety of Gaussian process based models". Its idea is to move along an ellipse that passes through the current state and a fresh draw from the prior, so every proposal is already plausible under the prior and only the likelihood needs checking.

The shrinking bracket always contains $\theta = 0$, the current state, so the loop terminates. Comparisons add a refinement. Because the probit factor is the probability that a Gaussian noise variable falls below the utility difference, the exact posterior is the marginal of a Gaussian restricted to a region cut out by linear constraints. Benavoli et al. (2021c) sample it with LinESS, a rejection-free elliptical slice sampler for linearly truncated Gaussians, at a cost of $O(n^3)$ time and $O(n^2)$ memory in the number of comparisons. Takeno et al. (2023) used Gibbs sampling, which redraws one variable at a time from its distribution given all the others, for the same distribution and found it faster: about 0.54 against 1.61 seconds in their Table 1, though LinESS should win when the number of truncations far exceeds the dimension.

With 200 samples the histogram is ragged and, with the default seed, the mean is off by almost 0.2: successive samples of a Markov chain are correlated, so 200 of them carry less information than 200 independent draws. With 1,000 the error is under 0.05, and with 5,000 the estimates agree with the exact column to about two digits. Sampling is the only method in this chapter whose error can be driven to zero by spending more computation. That is why studies of the other methods use it as the reference.

17.6 An exact answer: the skew Gaussian process #

The one-comparison posterior Equation (17.2) has a name and a closed form. Does the general posterior Equation (17.1) have one too? It does, and the reason is the same random-utility trick used throughout Chapter 16: a probit factor is the probability that a hidden Gaussian variable is positive.

A Gaussian vector conditioned on a linear transformation of itself, plus independent Gaussian noise, being positive is a unified skew-normal distribution. Durante (2019) showed that the posterior of parametric probit regression belongs to this family for any Gaussian prior, which makes the family conjugate to the probit likelihood. For functions, Benavoli et al. (2021c) proved that "the true posterior distribution of the preference function is a Skew Gaussian Process (SkewGP), with highly skewed pairwise marginals", and argued from it that Laplace's method "usually provides a very poor approximation". Theorem 29.1 states the theorem with its parameters.

The derivation also says where the posterior is skewed. Every factor depends on $\vf$ only through the differences $\mathbf{D}\vf$. In any direction orthogonal to all the $\mathbf{a}_k$, the likelihood is constant and the posterior is the Gaussian prior conditioned on the differences. The skew lives in at most $m$ directions, and never in the direction that shifts every utility by the same amount, because no difference changes along it. The figure below shows the smallest instance: two utilities and comparisons between them.