Gaussian Process Preference Learning

Chapter 16 turned an answer to "which do you prefer?" into a likelihood: the probit choice model, under which the probability of preferring one option grows with the difference in their utilities. Chapter 17 showed how to handle the non-Gaussian posterior that such a likelihood creates. This chapter joins the two with the Gaussian process prior of Chapter 8, which is what lets a few dozen answers about a few dozen designs say something about every design, including those never shown.

The result is the model of Chu and Ghahramani (2005). It is the default model behind preferential Bayesian optimization (PBO) two decades later, and it is what runs inside the figures of Chapter 19. We set it up, derive how to fit it, use it to predict, and then ask two questions that the regression model never raised: what comparisons cannot tell us at all, and how the pattern of which pairs were compared shapes what they can.

18.1 The model #

A person has a latent utility $g$ over a domain $\X$ of designs: the color of a poster, the parameters of an exoskeleton controller, the settings of a photo filter. We place a Gaussian process prior on it,

so that before any answer, the utilities of any finite set of designs are jointly Gaussian with covariances given by the kernel $k$ (Chapter 7). The kernel encodes the one assumption that makes learning from few answers possible: designs that are close in $\X$ have similar utilities.

The data are $m$ answers among $n$ distinct designs $\vx_1, \dots, \vx_n$. The $k$-th answer says that design $v_k$ was preferred to design $u_k$, where $v_k$ and $u_k$ are indices into the list of designs. A design can appear in many answers, so $n$ is at most $2m$ and often much smaller. Collect the utilities at the designs in a vector $\vf = (g(\vx_1), \dots, g(\vx_n))^\T$, whose prior is $\N(\mathbf{0}, \mK)$ with $[\mK]_{ij} = k(\vx_i, \vx_j)$. Each answer follows Thurstone's Case V (Equation (16.2)), and answers are independent given the utilities:

Chu and Ghahramani (2005) derived this likelihood exactly as Section 16.3 did: start from an ideal judge who prefers $v_k$ whenever $g(\vx_{v_k}) \ge g(\vx_{u_k})$, then contaminate each latent value with independent Gaussian noise of variance $\sigma^2$. The noise scale $\sigma$ is a hyperparameter, alongside the kernel's lengthscale and amplitude.

The posterior over $\vf$ is Equation (17.1): the Gaussian prior times Equation (18.1), normalized. It is not Gaussian, so Chu and Ghahramani (2005) replaced it by the Laplace approximation of Section 17.2. Predictions at new designs then follow from the Gaussian process exactly as in regression, because given $\vf$ the utility anywhere else is conditionally Gaussian.

The same model can be seen from a second angle. The likelihood depends only on differences, so one can define a function of pairs, $h(\vx, \vx') = g(\vx) - g(\vx')$, and treat each answer as a binary label on a pair. Because $h$ is a linear transformation of $g$, it is itself a Gaussian process, with covariance

by expanding the covariance of two differences term by term. Houlsby et al. (2011) (a preprint) derived this preference kernel and concluded that the model of Chu and Ghahramani "is equivalent to GPC with a particular class of kernels", Gaussian process classification on pairs. Every tool for classification, including the approximations of Chapter 17, therefore applies.

Preference learning of this kind reached people early. Brochu et al. (2007) used the same Thurstone model to help people find material appearances for computer graphics by choosing between rendered examples, and it has since been used, with the Laplace approximation, to learn which exoskeleton gaits users prefer from their comparisons and ordinal labels (Li et al., 2021). Chapter 24 and Chapter 25 work through two such problems end to end.

18.2 Fitting it #

Fitting the model means finding the posterior mode $\hat\vf$ and the curvature there. Equation (17.4) gave Newton's step for any likelihood; what it needs from the likelihood is its gradient and its negative Hessian $\mW$. For comparisons both have a structure that explains much of the rest of this chapter, so we derive them in full. Write $s = \sqrt{2}\,\sigma$, and for the $k$-th answer write $z_k = (f_{v_k} - f_{u_k})/s$ and $\mathbf{a}_k = \mathbf{e}_{v_k} - \mathbf{e}_{u_k}$, the vector with $+1$ at the winner, $-1$ at the loser, and 0 elsewhere, so that $z_k = \mathbf{a}_k^\T\vf/s$.

Look at what Equation (18.2) says entry by entry. The gradient pushes each design up by $\lambda(z_k)/s$ for every answer it won and down for every answer it lost. The push is large when the answer was surprising: $\lambda(z)$ grows roughly like $-z$ when $z$ is very negative, that is, when the current utilities say the loser should have won, and it vanishes when the answer was already expected. The matrix $\mW$ has, on its diagonal, the total weight of the answers each design took part in, and off the diagonal, minus the total weight of the answers between two designs. That is the weighted graph Laplacian of the comparison graph, whose nodes are designs and whose edges are answered pairs. Section 18.5 returns to what that means.

Two consequences follow immediately. Every $\mathbf{a}_k$ sums to zero, so the gradient sums to zero, $\mathbf{1}^\T\nabla\ell = 0$, and $\mW\mathbf{1} = \mathbf{0}$: raising every utility by the same amount changes no answer's probability. And the weights depend on $\vf$, so $\mW$ is recomputed at every Newton step.

With these two ingredients, the fit is the algorithm below. It is what fitPreference in this book's figure code does.

At the start, every $z_k = 0$, every answer is a coin flip under the model, and every weight is $\lambda(0)^2/s^2 = (2/\pi)/s^2$. The first step therefore treats all answers alike, and later steps reweight them by how surprising each turns out to be. In the figures of this chapter the fit converges in five to thirteen steps. Each step solves an $n \times n$ system, so the cost is $O(n^3)$ per step, the same as regression.

At the mode the gradient of the log posterior vanishes: $\nabla\ell(\hat\vf) - \mK^{-1}\hat\vf = \mathbf{0}$, so

which is equation (11) of Chu and Ghahramani (2005). The vector $\hat{\bm\beta}$ plays the role that $\bm\alpha = \mK^{-1}\vy$ played in regression (Equation (8.4)): one weight per design, positive for designs that won more than the model expected and negative for those that lost.

The same computations give the Laplace evidence Equation (17.5), which Chu and Ghahramani (2005) maximized over the kernel hyperparameters and $\sigma$ (their equation 12); BoTorch fits its preference model the same way (Section 18.6).

import numpy as np
from scipy.stats import norm

def terms(f, duels, s):
    """Gradient and negative Hessian of the pairwise probit log-likelihood."""
    n = len(f)
    grad, W = np.zeros(n), np.zeros((n, n))
    for i, j in duels:                        # design i was preferred to j
        z = (f[i] - f[j]) / s
        lam = np.exp(norm.logpdf(z) - norm.logcdf(z))   # phi(z) / Phi(z)
        grad[i] += lam / s
        grad[j] -= lam / s
        w = lam * (z + lam) / s**2
        W[i, i] += w; W[j, j] += w; W[i, j] -= w; W[j, i] -= w
    return grad, W

def fit_preference(K, duels, sigma=0.1, iters=50):
    n, s = len(K), np.sqrt(2) * sigma
    f = np.zeros(n)
    for _ in range(iters):
        grad, W = terms(f, duels, s)
        f_new = K @ np.linalg.solve(np.eye(n) + W @ K, W @ f + grad)
        done = np.max(np.abs(f_new - f)) < 1e-8
        f = f_new
        if done:
            break
    grad, W = terms(f, duels, s)              # beta and W at the mode
    return f, grad, W

18.3 Predicting preferences #

With the mode and the curvature in hand, the Laplace posterior over $\vf$ is Gaussian, $\N(\hat\vf, (\mK^{-1} + \hat\mW)^{-1})$, and predicting the utility at new designs is Gaussian conditioning, as in Section 8.1. For a new design $\vx$, write $\vk(\vx)$ for its prior covariances with the compared designs.

These match equations (17) and (18) of Chu and Ghahramani (2005) and Rasmussen and Williams's predictive equations for classification (Rasmussen and Williams, 2006, ch. 3), with the comparison $\hat\mW$ in place of the diagonal one. The second form in step 4 matters: $\hat\mW$ is a graph Laplacian, which always has a zero eigenvalue, so $\hat\mW^{-1}$ does not exist.

The quantity a preferential optimizer most often needs is the probability that a person will prefer one new design $\vx$ to another $\vx'$. The two latent values are jointly Gaussian under the approximation, with means $\mu, \mu'$, variances $\sigma^2(\vx), \sigma^2(\vx')$, and covariance $c$ from the same formula with $\vk(\vx')$ on the right. Their difference has mean $\mu - \mu'$ and variance $\sigma^2(\vx) + \sigma^2(\vx') - 2c$, and Equation (16.6) gives

equation (19) of Chu and Ghahramani (2005). The covariance $c$ matters. Two nearby designs have strongly correlated latent values, so their difference is much better known than either value, and the model can be confident about their order even when it is unsure how good either one is.

The figure below runs the model on the running objective of this book, the two-bump function of Chapter 11. A simulated person whose utility is that function answers each comparison with a little probit noise; you choose the pairs. Above the plot, the comparison graph draws one arc per answer.

Some things to try:

• Ask ten random pairs, twice. The figure starts with six answers, too few to find either bump. After 20 more the mean follows both bumps and places its maximum at $x = 0.72$, next to the true one at 0.73. Every comparison moved the curve only at its two ends and, through the kernel, their neighborhoods.
• Watch the band. It narrows near compared designs but stays wide everywhere, about $\pm 1$ even after 20 answers. That is not a failure; Section 18.4 explains it.
• Read the last prediction. Among the first answers the probabilities scatter, and some winners had been predicted to lose (a probability below one half); later most winners are predicted at 0.8 or above. A surprising answer moves the curve the most, because its $\lambda(z)$ is large.
• Change the lengthscale. At $\ell = 0.03$ each answer is a local bump, and the curve falls back toward zero wherever nothing nearby was compared. At $\ell = 0.5$ the curve is too stiff to hold two bumps, and its maximum lands in the wrong place.
• Press Best against random repeatedly. Every pair now includes the current best design, and the graph becomes a star. The region around the best is learned well and the rest poorly, a pattern that returns with the champion-and-challenger rules of Section 19.3 and with EUBO's collapse toward the current best (Section 19.6).

18.3.1 In two and more dimensions #

Nothing in Equation (18.1) to Equation (18.5) depends on the input dimension. Designs enter only through the kernel, as distances. The figure below runs the same model on a two-dimensional design space, with a simulated person whose utility is the Branin function, a standard test function with three equally good maxima along a curved valley, rescaled to the unit square with larger values better.

Some things to try:

• Compare the left and middle maps. With 20 random pairs among 17 designs the posterior mean recovers the valley's broad shape. In the default draw it orders 89% of all pairs of compared designs correctly, and its rank agreement with the hidden utility over the whole square is 0.87; press New draw a few times to see that agreement range from about 0.4 to 0.9.
• Look at the right map. Uncertainty is lowest where designs were compared and highest in the corners no comparison reached.
• Switch to Star. Every answer involves the same central design. Its comparisons with the others are learned, but the others are never compared with one another, and the model does worst here: averaged over 30 draws it orders 81% of the pairs correctly, with a rank agreement of 0.63, against 85% and 0.71 for random pairs.
• Switch to Isolated pairs. Twenty answers now touch 40 designs in 20 separate components. Each answer orders two designs and says nothing directly about how they compare with the other pairs, so the offsets between pairs come only from the kernel. With a kernel this smooth, that works: over 30 draws the model still orders 82% of the pairs correctly, and because 40 designs cover the square better than 17, its rank agreement over the square is slightly higher (0.79).
• Shorten the lengthscale to 0.08. The kernel now links designs less. The isolated pairs drop to 76% of pairs ordered correctly, the random pairs to 82%: the less the prior ties designs together, the more the answers must do it themselves.

How far does this go in more dimensions? The same code answers that. Table 18.1 reports an illustrative computation: 40 random designs in the unit cube, $m$ random pairs among them answered by a simulated person with noise $\sigma = 0.1$, an RBF kernel with lengthscale 0.52 (the mode of BoTorch's default prior, Section 18.6), and 20 repetitions. The utility is the Hartmann function in 3 and 6 dimensions, two standard benchmarks, and the 6-dimensional one hidden in 20 dimensions, where 14 inputs do not matter. The column "kernel value" is the average prior correlation between two of the designs.

Two patterns stand out. The model's grasp of new designs collapses with dimension: in 20 dimensions, random designs are so far apart in the kernel's eyes (an average kernel value of 0.006) that each answer informs almost nothing beyond its own two designs. Meanwhile its ranking of the designs it has seen keeps improving with more answers in every dimension, because those are ordered by direct comparisons. Making the lengthscale grow with the number of inputs $d$, to $0.3\sqrt{d}$, raised the average kernel value in 20 dimensions to 0.41 but improved the rank correlation on new designs only to 0.19, 0.30, and 0.28: the problem there is the 14 inputs that do not matter, which a kernel with one lengthscale shared by all inputs (an isotropic kernel) cannot ignore. These are one illustrative setup, not a benchmark (inference). They suggest why PBO in many dimensions depends on the kernel's structure, such as one lengthscale per input (Section 9.2), and on choosing pairs well (Chapter 19); Chapter 30 reports the research.

18.4 What can be identified #

The band in Figure 18.1 stayed wide after 20 answers, while the band of a regression model with 20 observations would have pinched at each one. The difference is not the approximation. It is what comparisons can measure.

Shift. Every answer's probability depends on a difference of utilities, so adding the same constant $c$ to every utility, $\vf \mapsto \vf + c\mathbf{1}$, leaves the likelihood unchanged. In Equation (18.2) this appears as $\mathbf{1}^\T\nabla\ell = 0$ and $\mW\mathbf{1} = \mathbf{0}$: no answer pushes along the direction $\mathbf{1}$, and no answer adds curvature along it. Only the prior says anything about the overall level, so the posterior keeps the prior's uncertainty about it, and that shared uncertainty is what keeps the band wide. The figure below shows the same model with the band redrawn for $g(\vx)$ minus its average over the domain, a quantity comparisons can measure.