Predictive Entropy Search for Efficient Global Optimization of Black-box Functions

Introduction

Bayesian optimization techniques form a successful approach for optimizing black-box functions . The goal of these methods is to find the global maximizer of a nonlinear and generally non-convex function $f$ whose derivatives are unavailable. Furthermore, the evaluations of $f$ are usually corrupted by noise and the process that queries $f$ can be computationally or economically very expensive. To address these challenges, Bayesian optimization devotes additional effort to modeling the unknown function $f$ and its behavior. These additional computations aim to minimize the number of evaluations that are needed to find the global optima.

Optimization problems are widespread in science and engineering and as a result so are Bayesian approaches to this problem. Bayesian optimization has successfully been used in robotics to adjust the parameters of a robot’s controller to maximize gait speed and smoothness as well as parameter tuning for computer graphics . Another example application in drug discovery is to find the chemical derivative of a particular molecule that best treats a given disease . Finally, Bayesian optimization can also be used to find optimal hyper-parameter values for statistical and machine learning techniques .

We take a Bayesian approach to the problem described above and use a probabilistic model for the latent function $f$ to guide the search and to select $\widetilde{\mathbf{x}}_{N}$. In this work we use a zero-mean Gaussian process (GP) prior for $f$ . This prior is specified by a positive-definite kernel function $k(\mathbf{x},\mathbf{x}^{\prime})$. Given any finite collection of points $\{\mathbf{x}_{1},\dots,\mathbf{x}_{n}\}$, the values of $f$ at these points are jointly zero-mean Gaussian with covariance matrix $\mathbf{K}_{n}$, where $[\mathbf{K}_{n}]_{ij}=k(\mathbf{x}_{i},\mathbf{x}_{j})$. For the Gaussian likelihood described above, the vector of concatenated observations $\mathbf{y}_{n}$ is also jointly Gaussian with zero-mean. Therefore, at any location $\mathbf{x}$, the latent function $f(\mathbf{x})$ conditioned on past observations $\mathcal{D}_{n}$ is then Gaussian with marginal mean $\mu_{n}(\mathbf{x})$ and variance $v_{n}(\mathbf{x})$ given by

where $\mathbf{k}_{n}(\mathbf{x})$ is a vector of cross-covariance terms between $\mathbf{x}$ and $\{\mathbf{x}_{1},\dots,\mathbf{x}_{n}\}$.

Bayesian optimization techniques use the above predictive distribution $p(f(\mathbf{x})|\mathcal{D}_{n})$ to guide the search for the global maximizer $\mathbf{x}_{\star}$. In particular, $p(f(\mathbf{x})|\mathcal{D}_{n})$ is used during the computation of an acquisition function $\alpha_{n}(\mathbf{x})$ that is optimized at each iteration to determine the next evaluation location $\mathbf{x}_{n+1}$. This process is shown in Algorithm 1. Intuitively, the acquisition function $\alpha_{n}(\mathbf{x})$ should be high in areas where the maxima is most likely to lie given the current data. However, $\alpha_{n}(\mathbf{x})$ should also encourage exploration of the search space to guarantee that the recommendation $\widetilde{\mathbf{x}}_{N}$ is a global optimum of $f$, not just a global optimum of the posterior mean. Several acquisition functions have been proposed in the literature. Some examples are the probability of improvement , the expected improvement or upper confidence bounds . Alternatively, one can combine several of these acquisition functions .

The acquisition functions described above are based on optimistic estimates of the latent function $f$ which implicitly trade off between exploiting the posterior mean and exploring based on the uncertainty. We instead follow the approach described in and aim to maximize the expected gain of information on the posterior distribution of the global maximizer $\mathbf{x}_{\star}$. In Section 2 we derive a rearrangement of this acquisition function and a corresponding approximation that we call Predictive Entropy Search (PES). PES is more accurate than the approximation used in . In Section 3 we empirically evaluate this claim on both synthetic and real-world problems and show that this leads to real gains in performance.

Predictive entropy search

We propose to follow the information-theoretic method for active data collection described in . We are interested in maximizing information about the location $\mathbf{x}_{\star}$ of the global maximum, whose posterior distribution is $p(\mathbf{x}_{\star}|\mathcal{D}_{n})$. Our current information about $\mathbf{x}_{\star}$ can be measured in terms of the negative differential entropy of $p(\mathbf{x}_{\star}|\mathcal{D}_{n})$. Therefore, our strategy is to select $\mathbf{x}_{n+1}$ which maximizes the expected reduction in this quantity. The corresponding acquisition function is

where $p(y|\mathcal{D}_{n},\mathbf{x},\mathbf{x}_{\star})$ is the posterior predictive distribution for $y$ given the observed data $\mathcal{D}_{n}$ and the location of the global maximizer of $f$. Intuitively, conditioning on the location $\mathbf{x}_{\star}$ pushes the posterior predictions up in locations around $\mathbf{x}_{\star}$ and down in regions away from $\mathbf{x}_{\star}$. Note that, unlike the previous formulation, this objective is based on the entropies of predictive distributions, which are analytic or can be easily approximated, rather than on the entropies of distributions on $\mathbf{x}_{\star}$ whose approximation is more challenging.

In this section we show how to approximately sample from the conditional distribution of the global maximizer $\mathbf{x}_{\star}$ given the observed data $\mathcal{D}_{n}$, that is,

Given a shift-invariant kernel $k$, Bochner’s theorem asserts the existence of its Fourier dual $s(\mathbf{w})$, which is equal to the spectral density of $k$. Letting $p(\mathbf{w})=s(\mathbf{w})/\alpha$ be the associated normalized density, we can write the kernel as the expectation

where $b\sim\mathcal{U}[0,2\pi]$. Let $\bm{\phi}(\mathbf{x})=\sqrt{2\alpha/m}\cos(\mathbf{W}\mathbf{x}+\mathbf{b})$ denote an $m$-dimensional feature mapping where $\mathbf{W}$ and $\mathbf{b}$ consist of $m$ stacked samples from $p(\mathbf{w},b)$. The kernel $k$ can then be approximated by the inner product of these features, $k(\mathbf{x},\mathbf{x}^{\prime})\approx\bm{\phi}(\mathbf{x})^{\mathsf{T}}\bm{\phi}(\mathbf{x}^{\prime})$. This approach was used by as an approximation method in the context of kernel methods. The feature mapping $\bm{\phi}(\mathbf{x})$ allows us to approximate the Gaussian process prior for $f$ with a linear model $f(\mathbf{x})=\bm{\phi}(\mathbf{x})^{\mathsf{T}}\bm{\theta}$ where $\bm{\theta}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$ is a standard Gaussian. By conditioning on $\mathcal{D}_{n}$, the posterior for $\bm{\theta}$ is also multivariate Gaussian, $\bm{\theta}|\mathcal{D}_{n}\sim\mathcal{N}(\mathbf{A}^{-1}\bm{\Phi}^{\mathsf{T}}\mathbf{y}_{n},\sigma^{2}\mathbf{A}^{-1})$ where $\mathbf{A}=\bm{\Phi}^{\mathsf{T}}\bm{\Phi}+\sigma^{2}\mathbf{I}$ and $\bm{\Phi}^{\mathsf{T}}=[\bm{\phi}(\mathbf{x}_{1})\dots\bm{\phi}(\mathbf{x}_{n})]$.

Let $\bm{\phi}^{(i)}$ and $\bm{\theta}^{(i)}$ be a random set of features and the corresponding posterior weights sampled both according to the generative process given above. They can then be used to construct the function $f^{(i)}(\mathbf{x})=\bm{\phi}^{(i)}(\mathbf{x})^{\mathsf{T}}\bm{\theta}^{(i)}$, which is an approximate posterior sample of $f$—albeit one with a finite parameterization. We can then maximize this function to obtain $\mathbf{x}_{\star}^{(i)}=\operatornamewithlimits{arg\,max}_{\mathbf{x}\in\mathcal{X}}f^{(i)}(\mathbf{x})$, which is approximately distributed according to $p(\mathbf{x}_{\star}|\mathcal{D}_{n})$. Note that for early iterations when $n<m$, we can efficiently sample $\bm{\theta}^{(i)}$ with cost $\mathcal{O}(n^{2}m)$ using the method described in Appendix B.2 of . This allows us to use a large number of features in $\bm{\phi}^{(i)}(\mathbf{x})$.

2 Approximating the predictive entropy

$f(\mathbf{x})$ is smaller than $f(\mathbf{x}_{\star})$. This simplified constraint only conditions on the given $\mathbf{x}$ rather than requiring $f(\mathbf{x}_{\star})\leq f(\mathbf{z})$ for all $\mathbf{z}\in\mathcal{X}$.

We incorporate these simplified constraints into $p(f(\mathbf{x})|\mathcal{D}_{n})$ to approximate $p(f(\mathbf{x})|\mathcal{D}_{n},\mathbf{x}_{\star})$. This is achieved by multiplying $p(f(\mathbf{x})|\mathcal{D}_{n})$ with specific factors that encode the above constraints. In what follows we briefly show how to construct these factors; more detail is given in Appendix B.

Implicitly we are assuming above that $\mathbf{f}$ depends on our observations and constraint C1.1, but is independent of C1.2 and C2 given $\mathbf{z}$. The computations necessary to obtain $\mathbf{m}_{\mathbf{f}}$ and $\mathbf{V}_{\mathbf{f}}$ are similar to those used above and in (1). The required quantities are similar to the ones used by EP to make predictions in the Gaussian process binary classifier . We can then incorporate C3 by multiplying $\mathcal{N}(\mathbf{f}|\mathbf{m}_{\mathbf{f}},\mathbf{V}_{\mathbf{f}})$ with a factor that guarantees $f(\mathbf{x})<f(\mathbf{x}_{\star})$. The predictive distribution for $f(\mathbf{x})$ given $\mathcal{D}_{n}$ and all the constraints can be approximated as

where $Z$ is a normalization constant. The variance of the right hand size of (8) is given by

3 Hyperparameter learning and the PES acquisition function

We now show how the previous approximations are integrated to compute the acquisition function used by predictive entropy search (PES). This acquisition function performs a formal treatment of the hyperparameters. Let $\bm{\psi}$ denote a vector of hyperparameters which includes any kernel parameters as well as the noise variance $\sigma^{2}$. Let $p(\bm{\psi}|\mathcal{D}_{n})\propto p(\bm{\psi})\,p(\mathcal{D}_{n}|\bm{\psi})$ denote the posterior distribution over these parameters where $p(\bm{\psi})$ is a hyperprior and $p(\mathcal{D}_{n}|\bm{\psi})$ is the GP marginal likelihood. For a fully Bayesian treatment of $\bm{\psi}$ we must marginalize the acquisition function (3) with respect to this posterior. The corresponding integral has no analytic expression and must be approximated using Monte Carlo. This approach is also taken in .

We draw $M$ samples $\{\bm{\psi}^{(i)}\}$ from $p(\bm{\psi}|\mathcal{D}_{n})$ using slice sampling . Let $\mathbf{x}_{\star}^{(i)}$ denote a sampled global maximizer drawn from $p(\mathbf{x}_{\star}|\mathcal{D}_{n},\bm{\psi}^{(i)})$ as described in Section 2.1. Furthermore, let $v_{n}^{(i)}(\mathbf{x})$ and $v_{n}^{(i)}(\mathbf{x}|\mathbf{x}_{\star}^{(i)})$ denote the predictive variances computed as described in Section 2.2 when the model hyperparameters are fixed to $\bm{\psi}^{(i)}$. We then write the marginalized acquisition function as

Note that PES is effectively marginalizing the original acquisition function (2) over $p(\bm{\psi}|\mathcal{D}_{n})$. This is a significant advantage with respect to other methods that optimize the same information-theoretic acquisition function but do not marginalize over the hyper-parameters. For example, the approach of approximates (2) only for fixed $\bm{\psi}$. The resulting approximation is computationally very expensive and recomputing it to average over multiple samples from $p(\bm{\psi}|\mathcal{D}_{n})$ is infeasible in practice.

Experiments

Figure 1 shows the objective functions produced by RS, ES and PES for a particular $\mathcal{D}_{n}$ with 10 measurements whose locations are selected uniformly at random in $^{2}$. The locations of the collected measurements are displayed with an “x” in the plots. The particular objective function used to generate the measurements in $\mathcal{D}_{n}$ is displayed in the left part of Figure 2. The plots in Figure 1 show that the PES approximation to (2) is more similar to the ground truth given by RS than the approximation produced by ES. In this figure we also see a discrepancy between RS and PES at locations near $\mathbf{x}=(0.572,0.687)$. This difference is an artifact of the discretization used in RS. By zooming in and drawing many more samples we would see the same behavior in both plots.

In these experiments we compared the performance of PES with that of ES and expected improvement (EI) , a widely used acquisition function in the Bayesian optimization literature. We again assume that the optimal hyper-parameter values are known to all methods. Predictive performance is then measured in terms of the immediate regret (IR) $|f(\widetilde{\mathbf{x}}_{n})-f(\mathbf{x}_{\star})|$, where $\mathbf{x}_{\star}$ is the known location of the global maximum and $\widetilde{\mathbf{x}}_{n}$ is the recommendation of each algorithm had we stopped at step $n$—for all methods this is given by the maximizer of the posterior mean. The right plot in Figure 2 shows the decimal logarithm of the median of the IR obtained by each method across the 1000 different objective functions. Confidence bands equal to one standard deviation are obtained using the bootstrap method. Note that while averaging these results is also interesting, corresponding to the expected performance averaged over the prior, here we report the median IR because the empirical distribution of IR values is very heavy-tailed. In this case, the median is more representative of the exact location of the bulk of the data. These results indicate that the best method in this setting is PES, which significantly outperforms ES and EI. The plot also shows that in this case ES is significantly better than EI.

We perform another series of experiments in which we optimize well-known synthetic benchmark functions including a mixture of cosines and Branin-Hoo (both functions defined in $^{2}$) as well as the Hartmann-6 (defined in $^{6}$) . In all instances, we fix the measurement noise to $\sigma^{2}=10^{-3}$. For both PES and EI we marginalize the hyperparameters $\bm{\psi}$ using the approach described in Section 2.3. ES, by contrast, cannot average its approximation of (2) over the posterior on $\bm{\psi}$. Instead, ES works by fixing $\bm{\psi}$ to an estimate of its posterior mean (obtained using slice sampling) . To evaluate the gains produced by the fully Bayesian treatment of $\bm{\psi}$ in PES, we also compare with a version of PES (PES-NB) which performs the same non-Bayesian (NB) treatment of $\bm{\psi}$ as ES. In PES-NB we use a single fixed hyperparameter as in previous sections with value given by the posterior mean of $\bm{\psi}$. All the methods are initialized with three random measurements collected using latin hypercube sampling .

The plots in Figure 3 show the median IR obtained by each method on each function across 250 random initializations. Overall, PES is better than PES-NB and ES. Furthermore, PES-NB is also significantly better than ES in most of the cases. These results show that the fully Bayesian treatment of $\bm{\psi}$ in PES is advantageous and that PES can produce better approximations than ES. Note that PES performs better than EI in the Branin and cosines functions, while EI is significantly better on the Hartmann problem. This appears to be due to the fact that entropy-based strategies explore more aggressively which in higher-dimensional spaces takes more iterations. The Hartmann problem, however, is a relatively simple problem and as a result the comparatively more greedy behavior of EI does not result in significant adverse consequences. Note that the synthetic functions optimized in the previous experiment were much more multimodal that the ones considered here.

We finally optimize different real-world cost functions. The first one (NNet) returns the predictive accuracy of a neural network on a random train/test partition of the Boston Housing dataset . The variables to optimize are the weight-decay parameter and the number of training iterations for the neural network. The second function (Hydrogen) returns the amount of hydrogen production of a particular bacteria in terms of the PH and Nitrogen levels of the growth medium . The third one (Portfolio) returns the ratio of the mean and the standard deviation (the Sharpe ratio) of the 1-year ahead returns generated by simulations from a multivariate time-series model that is adjusted to the daily returns of stocks AXP, BA and HD. The time-series model is formed by univariate GARCH models connected with a Student’s $t$ copula . These three functions (NNet, Hydrogen and Portfolio) have as domain $^{2}$. Furthermore, in these examples, the ground truth function that we want to optimize is unknown and is only available through noisy measurements. To obtain a ground truth, we approximate each cost function as the predictive distribution of a GP that is adjusted to data sampled from the original function (1000 uniform samples for NNet and Portfolio and all the available data for Hydrogen ). Finally, we also consider another real-world function that returns the walking speed of a bipedal robot . This function is defined in $^{8}$ and its inputs are the parameters of the robot’s controller. In this case the ground truth function is noiseless and can be exactly evaluated through expensive numerical simulation. We consider two versions of this problem (Walker A) with zero-mean, additive noise of $\sigma=0.01$ and (Walker B) with $\sigma=0.1$.

Figure 4 shows the median IR values obtained by each method on each function across 250 random initializations, except in Hydrogen where we used 500 due to its higher level of noise. Overall, PES, ES and PES-NB perform similarly in NNet, Hydrogen and Portfolio. EI performs rather poorly in these first three functions. This method seems to make excessively greedy decisions and fails to explore the search space enough. This strategy seems to be advantageous in Walker A, where EI obtains the best results. By contrast, PES, ES and PES-NB tend to explore more in this latter dataset. This leads to worse results than those of EI. Nevertheless, PES is significantly better than PES-NB and ES in both Walker datasets and better than EI in the noisier Walker B. In this case, the fully Bayesian treatment of hyper-parameters performed by PES produces improvements in performance.

Conclusions

We have proposed a novel information-theoretic approach for Bayesian optimization. Our method, predictive entropy search (PES), greedily maximizes the amount of one-step information on the location $\mathbf{x}_{\star}$ of the global maximum using its posterior differential entropy. Since this objective function is intractable, PES approximates the original objective using a reparameterization that measures entropy in the posterior predictive distribution of the function evaluations. PES produces more accurate approximations than Entropy Search (ES), a method based on the original, non-transformed acquisition function. Furthermore, PES can easily marginalize its approximation with respect to the posterior distribution of its hyper-parameters, while ES cannot. Experiments with synthetic and real-world functions show that PES often outperforms ES in terms of immediate regret. In these experiments, we also observe that PES often produces better results than expected improvement (EI), a popular heuristic for Bayesian optimization. EI often seems to make excessively greedy decisions, while PES tends to explore more. As a result, EI seems to perform better for simple objective functions while often getting stuck with noisier objectives or for functions with many modes.

References

Appendix A Details on approximating GP sample paths

In this section we give further details about the approach used in Section 2.1 to approximate a GP using random features. These random features can be used to approximate sample paths from the GP posterior. By optimizing these sample paths we obtain posterior samples over the global maxima $\mathbf{x}_{\star}$. We derive in more detail the kernel approximation from (5). Formally, the theorem of states

A continuous, shift-invariant kernel is positive definite if and only if it is the Fourier transform of a non-negative, finite measure.

As a result given some kernel $k(\mathbf{x},\mathbf{x}^{\prime})=k(\mathbf{x}-\mathbf{x}^{\prime},\mathbf{0})$ there must exist an associated density $s(\mathbf{w})$, known as its spectral density, which is the Fourier dual of $k$. This can be written as

Further, we can treat this measure as a probability density $p(\mathbf{w})=s(\mathbf{w})/\alpha$ where $\alpha=\int s(\mathbf{w})\,d\mathbf{w}$ is the normalizing constant. Consequently, the kernel can be written as

We now briefly show the equivalence between a Bayesian linear model using random features $\bm{\phi}$ and a GP with kernel $k$. Consider now a linear model $f(\mathbf{x})=\phi(\mathbf{x})^{\mathsf{T}}\bm{\theta}$ where $\bm{\theta}\sim\mathcal{N}(\mathbf{0},\mathbf{I})$ has a standard Gaussian distribution and observations $\mathcal{D}_{n}=\{(\mathbf{x}_{i},y_{i})\}_{i\leq n}$ of the form $y_{i}\sim\mathcal{N}(f(\mathbf{x}_{i}),\sigma^{2})$. The posterior of $\bm{\theta}$ given $(\mathcal{D}_{n},\bm{\phi})$ is also be normal $\mathcal{N}(\mathbf{m},\mathbf{V})$ where

and where $[\bm{\Phi}]_{i}=\bm{\phi}(\mathbf{x}_{i})$ and $[\mathbf{y}]_{i}=y_{i}$ consist of the stacked features and observations respectively. We can also easily write the predictive distribution over $f$ evaluated at a test point $\mathbf{x}$, which is Gaussian distributed with mean and variance given by

By a simple application of the matrix-inversion lemma these quantities can be rewritten in terms which only make use of the inner products between features,

the expectations of which are equivalent to the kernel $k$ and we obtain the same expressions as that in (1).

Appendix B Details on approximating the predictive variance

We now provide further details on approximating the predictive variance $v_{n}(\mathbf{x}|\mathbf{x}_{\star})$ of inputs $\mathbf{x}$ given the position of the global optimizer $\mathbf{x}_{\star}$. In particular we include all steps omitted in the presentation of Section 2.2.

Here $\mathbf{c}$ contains the random variables that we will condition on in order to enforce constraint C1.1. Given the input locations $\mathbf{x}$ and $\mathbf{x}_{\star}$ we can construct a kernel matrix $\mathbf{K}$ containing the covariance evaluated on the stacked vector $[\mathbf{z};\mathbf{c}]$. We again refer to in constructing this matrix which includes derivative observations, the computations of which are tedious but not overly complicated. Note also that the portions of $\mathbf{K}$ which correspond to $y_{i}$ will have an additional $\sigma^{2}$ due to the observation noise. Next let $\mathbf{K}_{\mathbf{z}}$, $\mathbf{K}_{\mathbf{c}}$, and $\mathbf{K}_{\mathbf{z}\mathbf{c}}$ denote the corresponding diagonal and off-diagonal blocks of the kernel matrix. We can now condition on the observed values of $\mathbf{c}$ to write

where $\mathbf{m}_{0}=\mathbf{K}_{\mathbf{z}\mathbf{c}}\mathbf{K}^{-1}_{\mathbf{c}}\mathbf{c}$ and $\mathbf{V}_{0}=\mathbf{K}_{\mathbf{z}}-\mathbf{K}_{\mathbf{z}\mathbf{c}}\mathbf{K}^{-1}_{\mathbf{c}}\mathbf{K}_{\mathbf{z}\mathbf{c}}^{\mathsf{T}}$.

B.2 Incorporating the non-analytic latent constraints (C1.2 and C2)

The additional constraints C1.2 and C2 can be introduced explicitly as in (6), which takes the form of a single Gaussian factor and $d+1$ non-Gaussian factors

We approximate this distribution using a single multivariate Gaussian $q(\mathbf{z})$ where each non-Gaussian factor is replaced by a Gaussian approximation $\widetilde{t}_{i}(z_{i})=\mathcal{N}(z_{i};\widetilde{m}_{i},\widetilde{v}_{i})$ such that

where this approximation is parameterized by $\mathbf{m}=\mathbf{V}[\widetilde{\mathbf{V}}^{-1}\widetilde{\mathbf{m}}+\mathbf{V}_{0}^{-1}\mathbf{m}_{0}]$ and $\mathbf{V}=(\widetilde{\mathbf{V}}^{-1}+\mathbf{V}_{0}^{-1})^{-1}$. The parameters of the approximate factors are combined to form the vector $[\widetilde{\mathbf{m}}]_{i}=\widetilde{m}_{i}$ and the diagonal matrix $[\widetilde{\mathbf{V}}]_{ii}=\widetilde{v}_{i}$.

To compute the approximate factors we use expectation propagation (EP). EP is a procedure that starts from some initial values for the approximate factors $(\widetilde{m}_{i},\widetilde{v}_{i})$ and iteratively refines these quantities; here we initialize $\widetilde{m}_{i}=0$ and $\widetilde{v}_{i}=\infty$ which corresponds to $\mathbf{m}=\mathbf{m}_{0}$ and $\mathbf{V}=\mathbf{V}_{0}$. At each iteration, for every factor $i$, we remove the contribution of the $i$th approximate factor to form the cavity distribution $q_{\setminus i}(\mathbf{z})\propto q(\mathbf{z})/\widetilde{t}_{i}(z_{i})$. Given the independent factors we consider here we can focus on each individual component $q_{\setminus i}(z_{i})$ separately with mean and variance

For both sets of constraints used in this work the moments can easily be obtained by computing the normalizing constant $\bar{Z}_{i}=\int\mathcal{N}(z_{i}|\bar{m}_{i},\bar{v}_{i})\,t_{i}(z_{i})\,dz_{i}$ and using the following identities:

B.3 Incorporating the prediction constraint (C3)

Given some test input $\mathbf{x}$ we now turn to the problem of making predictions about $f(\mathbf{x})$. We again note that both the “prior” terms $\mathbf{m}_{0}$, $\mathbf{V}_{0}$ and the EP factors, $\widetilde{\mathbf{m}}$ and $\widetilde{\mathbf{V}}$, are independent of $\mathbf{x}$ and can be precomputed once for later use at prediction time.

Let $\mathbf{f}=[f(\mathbf{x});f(\mathbf{x}_{\star})]$ be a vector given by the concatenation of the latent function at $\mathbf{x}$ and $\mathbf{x}_{\star}$. The distribution for $\mathbf{f}$ given the first two constraints can be written as

By writing $p(\mathbf{f}|\mathbf{z},\mathbf{c})$ above we are assuming that $\mathbf{f}$ is independent of C1.2 and C2 given $\mathbf{z}$ and as a result the above is simply an integral over the product of two Gaussians. Let $\mathbf{K}_{\dagger}$ be the cross-covariance matrix evaluated between $\mathbf{f}$ and $[\mathbf{z};\mathbf{c}]$ and $\mathbf{K}_{\mathbf{f}}$ the covariance matrix associated with $\mathbf{f}$. The posterior above will then be Gaussian with mean and variance

where $\widetilde{\mathbf{W}}$ is a block-diagonal matrix where the first block is zero and the second is $\widetilde{\mathbf{V}}$ (note this matrix is also diagonal since $\widetilde{\mathbf{V}}$ is diagonal). Finally, these values can be plugged into (8–9) in order to arrive at $v_{n}(\mathbf{x}|\mathbf{x}_{\star})$.