A Framework for Interdomain and Multioutput Gaussian Processes

One obstacle to the use of Gaussian processes (GPs) in large-scale problems, and as a component in deep learning system, is the need for bespoke derivations and implementations for small variations in the model or inference. In order to improve the utility of GPs we need a modular system that allows rapid implementation and testing, as seen in the neural network community. We present a mathematical and software framework for scalable approximate inference in GPs, which combines interdomain approximations and multiple outputs. Our framework, implemented in GPflow, provides a unified interface for many existing multioutput models, as well as more recent convolutional structures. This simplifies the creation of deep models with GPs, and we hope that this work will encourage more interest in this approach.


Introduction
Gaussian processes (GPs) [Rasmussen and Williams, 2006] are probability distributions over functions with many properties which make them convenient for use in Bayesian models. They can be used as both priors and approximate posteriors in a wide range of models where some kind of function needs to be learned, such as classification, regression, or state-space models. While GPs are usually formulated as having a single output, multiple interacting outputs can be considered as well [Álvarez et al., 2012]. This is crucial for building deep Gaussian process models where we want to be able to learn more complex intermediate mappings than just simple 1D warpings.
The most well-known drawback of GPs is that exact implementations scale with the number of training points N as O(N 3 ) in time. Multioutput GPs (MOGPs) have the additional problem of also needing to consider the covariances between all P outputs, leading to O(N 3 P 3 ) scaling. Fortunately, the development of variational inducing variable approximations has reduced the cost of expensive matrix manipulations to O(M 3 ), where M N , and has allowed for minibatch training [Titsias, 2009;Hensman et al., 2013]. These techniques have also been extended for use in MOGPs [Álvarez et al., 2010] and provided large improvements. Using interdomain approximations [Lázaro-Gredilla and Figueiras-Vidal, 2009], a wide range of efficient inference schemes can be implemented for MOGPs, which can reduce the computational complexity to O(M 3 P ).
In this work, we give a comprehensive introduction to sparse variational Gaussian process models, interdomain approximations, and multioutput Gaussian processes and show how they fit together in a unified description. We then present a unified software framework for specifying interdomain approximations and multioutput models, which we have implemented in GPflow . The flexible specification of interdomain approximations is particularly important for allowing the implementation of a wide range of efficient inference schemes for MOGPs, but our framework applies equally to single-output GPs. Our immediate practical goal is to flexibly allow convolutional Gaussian process  structure in deep GPs [Damianou and Lawrence, 2013;Salimbeni and Deisenroth, 2017], which requires image-to-image GPs [Blomqvist et al., 2018;Dutordoir et al., 2020]. More generally, we hope that by removing obstacles to implementing and using multioutput GPs, more researchers will experiment with architectues that contain them, and that more practitioners will apply them to problems that can benefit from them.
X ∈ R N ×D . A Gaussian process is defined by the property that its function values at arbitrary locations (here denoted X) are jointly Gaussian distributed as We usually take µ(·) = 0 without loss of generality. Throughout this document, we will manipulate a whole function f (·) as an infinitely long vector. Instead of having a finite number of elements addressed by a discrete index, we have an element for each continuous-valued input. From this view, f (x) picks a single element out of the infinitely long vector f (·), while f (X) picks a subset of points. We will abuse notation by denoting the Gaussian process distribution over f (·) as a density p(f (·)). Strictly speaking, it is not possible to write a density over an infinitely long vector. However, finite-dimensional marginals of the GP are well defined by densities (eq. 1), so we consider p(f (·)) to be an object that gives us the correct finite-dimensional Gaussian distribution when we marginalise out all but a finite set of points. 1 We denote this process as where f X (·) denotes the infinitely long vector f (·) excluding f (X). Since we work with finite datasets, we will always evaluate the GP at a finite number of locations, and therefore only need to consider a finite-dimensional marginal. We will also often factorise p(f (·)) into finite and infinite parts, e.g.
This splits the distribution over the infinitely long f (·) into the finite-dimensional Gaussian p(f (X)) (eq. 1), and the infinitely long conditional p(f X (·) | f (X)) that includes all function values except f (X). The conditional p(f X (·) | f (X)) is given by where [k f · ] n = C f [f (x n ), f (·)] = k(x n , ·). This expression is identical to the conditional of the whole process p(f (·) | f (X)). In this latter case, the random variable f (X) is included in f (·) on the left hand side of the conditioning and has a deterministic (delta function) relationship to the conditioned value. We explicitly denote f X (·) in eq. 4 so we do not double-count random variables on the right hand side of the equality.

Exact inference
We are interested in Bayesian models that use GPs as prior distributions. Common uses are regression, classification, or unsupervised learning (e.g. density estimation). Initially, our running example will be regression with arbitrary likelihoods. In our setup, we model N observations y = {y n } N n=1 that are produced by evaluating the GP at input locations x n . We denote our model as f (·) ∼ GP(0, k(·, · )) , y n | f (·), x n ∼ p(y n | f (x n )) , where we assume that we can evaluate the likelihood density p(y n | f (x n )).
When performing inference, we are firstly interested in the posterior for making predictions, and secondly in the marginal likelihood for learning hyperparameters. One of the early reasons for interest in GPs was that closed-form solutions for these quantities exist when the likelihood is Gaussian, p(y n | f (x n )) = N (y n ; f (x n ), σ 2 ). We can obtain the posterior for any predictive point using Bayes' rule, and the marginal likelihood through marginalising over the prior: p(f (·) | y, X) = p(f (·)) N n=1 p(y n | f (x n )) p(y | X) = p(f X (·) | f (X)) p(f (X)) N n=1 p(y n | f (x n )) p(y | X) = GP k T f · K ff + σ 2 I −1 y, k(·, · ) − k T f · K ff + σ 2 I −1 k f · , p(y | X) = N y; 0, K ff + σ 2 I , where [k f · ] n = C f [f (x n ), f (·)] = k(x n , ·). To find the predictive distribution for a specific set of inputs X * , we need to marginalise the posterior GP from eq. 6 to remove all points but X * : p(f (X * ) | y, X) = p(f X (·) | f (X)) p(f (X)) N n=1 p(y n | f (x n )) p(y | X) d{f X * (·)} = p(f (X * ) | f (X)) p(f (X)) N n=1 p(y n | f (x n )) p(y | X) df (X).
In GPflow, this model is implemented as GPR (Gaussian Process Regression). Note that you must include the likelihood variance when predicting actual observations, rather than the unobserved function values [Rasmussen and Williams, 2006]. Despite having closed-form solutions to the quantities of interest, actual computation using these equations quickly becomes intractable for even moderately sized datasets. The cause is the O N 3 computational cost for computing the inverse and determinant of K. The source of this intractability is unusual. Usually, Bayesian models have tractable priors and likelihoods, and any intractability comes from the fact that there is no closed-form solution for the normalising constant of the posterior (the marginal likelihood). The computational expense comes from needing to use methods like MCMC that avoid dealing with normalising constants. In GPs, however, it is the prior p(f (X)) that is intractable, as this density already contains the costly inversion and determinant calculations.
Non-Gaussian likelihoods present an additional challenge, due to the loss of conjugacy, which allowed the closed-form computation of eqs. 6 and 7. To summarise, any approximate method for inference needs to deal with 1) the cost of computing the prior, and therefore also the posterior, 2) likelihoods that are not conjugate to the Gaussian distribution. We will address this in the following section.

Approximate inference
Many approximations that address the two problems above have been proposed over the years. These can largely be categorised in terms of approximations to the model (see Quiñonero-Candela and Rasmussen [2005] for a review), and approximations to the posterior. When approximating Gaussian process models, we need to be careful that the approximation retains the desirable properties of the original model. Van der Wilk [2019, ch. 2] discussed desiderata that approximations to non-parametric models should satisfy, and argued that it was particularly important to maintain the uncertainty estimates afforded by using a non-parametric model. For this reason, we favour posterior approximations.
Broadly speaking, there are currently two general inference frameworks for posterior approximations to GPs that address the two problems above: expectation propagation (EP) variants [Minka, 2001;Csató and Opper, 2002;Bui et al., 2017], and variational inference (VI) methods [Blei et al., 2017;Titsias, 2009;Hensman et al., 2013;Matthews et al., 2016]. Both methods provide procedures for picking an approximate posterior from a family of tractable ones, and both methods use the same family of tractable approximate posteriors [Bui et al., 2017]. Currently, there is no consensus on which method is strictly better. In terms of pure predictive performance, neither method consistently outperforms the other. Bui et al. [2017] show that for approximations with few inducing points, interpolations between EP and VI often perform best, while Hensman et al. [2015b] provide favourable evidence for VI in classification. Bauer et al. [2016] and Van der Wilk [2019] compare EP (in its form as the FITC model approximation [Snelson and Ghahramani, 2006;Quiñonero-Candela and Rasmussen, 2005]) and VI more closely in terms of the quality of posterior approximation. They show that EP sometimes provides a posterior approximation with significantly different properties to the true posterior. On the other hand, VI was shown to approximate the true posterior in a predictable way. Moreover, VI can easily be extended to importance-weighted variational inference which has been shown to converge to the true marginal likelihood [Domke and Sheldon, 2018]. In this work, we focus on variational inference, mainly because of its practical convenience, and its strong theoretical guarantees including satisfying the desiderata from Van der Wilk [2019].
Variational inference works by minimising the KL divergence KL[q(f (·)) p(f (·) | y, X)] from the true posterior p(f (·) | y, X) to a distribution in the tractable family q(f (·)). While the KL divergence cannot directly be computed, we can minimise it by maximising a lower bound L to the log marginal likelihood, which has the KL as its gap: The lower bound L depends on the variational approximation q(f (·)), as well as the training data and any hyperparameters. Minimising the KL as a means to choosing the approximate posterior provides some important properties. Firstly, since the algorithm requires only optimisation, it is guaranteed to converge. Secondly, subject to optimisation difficulties, we can only improve the approximation by expanding the family of approximating distributions [Bauer et al., 2016;Matthews, 2017]. Finally, variational bounds are relatively easy to specify for a wide range of models, and very general algorithms exist to optimise them. These algorithms are based on stochastic evaluation of the bound [Hoffman et al., 2013] and its gradients through either score function estimators [Ranganath et al., 2014] or the reparameterisation trick [Kingma and Welling, 2013;Rezende et al., 2014;Titsias and Lázaro-Gredilla, 2014].

Inducing point posteriors
The variational lower bound L defined in eq. 9 depends on the family of approximate posteriors from which we will choose our approximation q(f (·)). We need a family that contains distributions that are convenient to manipulate mathematically, as well as computationally tractable. Additionally, we would like the family to be large enough to contain a distribution very close to the true posterior. We choose the family of Gaussian process posteriors introduced by Titsias [2009].
To construct a family of posteriors, we have to parameterise a distribution over functions, for which Gaussian processes are a convenient choice. We start by specifying the joint distribution over M function values at the input locations Z = {z m } M m=1 . We collect the function values f (Z) in the variable u ∈ R M , and specify a free mean and variance: Next, we specify the distribution for all other points f Z (·) using the prior conditioned on u.
Together, this gives a Gaussian process on f (·): where [k ·u ] m = k(·, z m ) and [K uu ] mm = k(z m , z m ). This family of posteriors is computationally tractable, as only a smaller M × M matrix K uu has to be inverted for making predictions.
To gain more intuition into this family of posteriors, we provide an alternative construction, similar to the unifying view from Bui et al. [2017]. We can consider the family of approximate posteriors from eq. 11 to be the collection of all the posteriors that we can get from observing M function values through an arbitrary Gaussian likelihood. To show this, we derive eq. 11 by setting up a regression problem akin to eq. 6 only with the likelihoodq(m | u) = N (m ; u, S ). We get where Z is the marginal likelihood for this surrogate regression problem that normalises q(f (·)). Following eq. 6, we obtain the approximate posterior from eq. 11 where m = m and S = (S −1 + K −1 uu ) −1 . This view provides a different way to intuitively understand the range of posteriors we allow. Following this interpretation, the variational method tries to find a set of M observations m which, together with the right likelihood variance S , give a posterior as similar as possible to the true posterior in terms of KL. This is often possible, particularly in the case where the data contain redundant information in a specific region of the input space. A concrete example of this is where we have many noisy observations in a concentrated region. In this situation, a few very precise observations can give the same posterior as many noisy observations. We can also immediately see that when M = N and a Gaussian likelihood is used in the model, the true posterior is contained within the family of distributions. Burt et al. [2019] show that in this case an arbitrarily good posterior can be found as N → ∞, with M N and the exact rate depending on the expected eigenvalue decay of K ff .
It is also worth pointing out that the approximate posterior has error bars that behave in the same way as the full model. That is, both the full model, and our approximate posterior predict with an infinite number of basis functions, allowing the posterior to be uncertain in regions unconstrained by the data. This is not the case for many parametric approximations (see Quiñonero-Candela and Rasmussen [2005] or Van der Wilk [2019, ch. 2]).
Finally, it is worth pointing out that the choice of this family of posteriors was not accidental. This family allows the KL divergence between the q(f (·)) and p(f (·) | y, X) to be finite. Many distributions over functions, including many Gaussian processes, would give an infinite KL. By choosing a family consisting of distributions that are all posteriors with respect to the prior of the model, we can guarantee a finite KL [Matthews, 2017].

Variational inference for Gaussian processes
Now that we have specified the family of approximate posteriors, we need to find the variational lower bound. Matthews et al. [2016] showed that, despite technically not being able to deal with densities over functions, variational inference in GP models proceeds very similarly to the well-understood parametric case. The key observation is that the KL divergence between the approximate distribution over functions (eq. 11) and the true posterior can be computed from the KL between a finite marginal distribution on function values: We can informally see that this makes sense, since we can include any other set of function values (denoted f (·)) without changing the total KL divergence: Substituting this KL divergence into eq. 9 and applying Bayes' rule to the true posterior gives the bound discussed by Hensman et al. [2013] and Hensman et al. [2015b]: The original variational GP method by Titsias [2009] is a special case for regression, where q(u) is analytically optimised, and substituted in. The bound takes the same form as familiar bounds for parametric models, with an expected log-likelihood term and a KL to the prior. Conveniently, the final bound is no more complex than that of parametric models, despite the technical difficulties of dealing with distributions over functions.
To evaluate the bound, we need to take expectations over the approximate posterior of a single function value, given by the univariate Gaussian distribution: To reduce the cost of summing over large numbers of data points, Hensman et al. [2013] suggested using minibatches to obtain an unbiased estimate of the sum. Other likelihoods can be handled by estimating the expected likelihood using Monte Carlo (unbiased), or Gauss-Hermite quadrature (small bias, but zero variance) [Hensman et al., 2015b].

Interdomain Gaussian processes
When working with Gaussian processes, we are usually concerned with observations and predictions of the function f (·) itself. However, point observations may not be the only properties of a function that we are interested in, or that we can observe. For example, when modelling the trajectory of an object, its velocity can be observed as well as location, which gives additional observations of the function's derivative [O'Hagan, 1992;Rasmussen and Williams, 2006, ch. 9.4]. Alternatively, in Bayesian quadrature, we are interested in learning about the integral of a function from point observations [Rasmussen and Ghahramani, 2003]. These quantities rely on a transformation of f (·) to a different domain, and we will name them "interdomain" variables, following Lázaro-Gredilla and Figueiras-Vidal [2009]. We start by discussing how observations of these interdomain quantities can be incorporated into an exact posterior. However, we are primarily interested in using interdomain transformations to define inducing variables for specifying GP approximations, which we will discuss afterwards. Here, we give an informal overview of interdomain methods, with a more formal discussion of some difficulties surrounding measurability being left to Matthews et al. [2016] and Matthews [2017].

Interdomain observations
The earlier examples of derivative or integral observations can be seen as observations of f (·) ∈ R X after having been passed through some linear operator L : R X → R. In other words, L takes the entire function f (·) and returns a real value. We modify the likelihood from p(y n | f (·), x n ) to p(y n | L n f (·)), where L n is the linear operator for the nth observation. A common example for the linear operator L n is an integral operator, as used by Lázaro-Gredilla and Figueiras-Vidal [2009]: We can see that L n indeed maps the entire function f (·) to a single real value. A different example of a linear operator which we can express, is the derivative of the dth input dimension as In this case, the linear operator only depends on a small neighbourhood around x n . Since an operator is allowed to depend on the entire function f (·), it is also free to only depend on a part of it while fitting within the operator framework. In fact, even regular point observations fit within this framework, where the linear operator simply depends on a single point: We will now derive the posterior for these operator observations. To ease notation, we redefine the operator L to return all N quantities needed for the likelihood terms. This makes L : To find the posterior, we first need to find the joint distribution between variables on f (·), and Lf (·), which the likelihood depends on. Given that the dependence between f (·) and Lf (·) is linear, and that the prior on f (·) is Gaussian, the joint will be Gaussian as well. Keeping in mind the usual shorthand of f (·) referring to function values at arbitrary inputs, we obtain where the covariances k L· and K LL are defined as These covariances can be evaluated by using the definition of Lf (·), taking advantage of the linearity of L, and then taking expectations over f (·): where with L j we denote the linear operator that works on the · parameter of k(·, · ). A specific example is when all elements in L use the integral operator of eq. 18, which gives We now have all the required elements for obtaining the posterior for interdomain observations. Using the fact that p(f (·), Lf (·)) is jointly Gaussian, everything reduces to Gaussian conditioning as for eq. 6. We obtain the posterior 2 which differs only from the point observation posterior (eq. 6) through the use of the interdomain covariances k L· and K LL .

Interdomain inducing variables
Apart from being quantities of interest for observation or prediction, interdomain variables can also be used to specify approximate posteriors.

used interdomain variables tailored to a particular model construction, while Lázaro-Gredilla and Figueiras-Vidal
[2009] discussed interdomain approximations for a given squared exponential kernel. Both these early approaches used a FITC style approach, while Álvarez et al. [2010] was the first to use an interdomain posterior within the variational framework. The construction of approximate posteriors using interdomain variables follows the exact same procedure as using inducing points. In eq. 12 we showed that we specified our family of approximate posteriors by conditioning on M arbitrary point observations. By changing the covariances as we did in eq. 27, we can specify our approximate posterior through conditioning on interdomain variables instead. The final form of the approximate posterior (eq. 11) is changed only through K uu and k u· being evaluations of the interdomain covariances of eq. 22. The variational lower bound remains unchanged, with the KL[q(u) p(u)] term now being between the interdomain variables [Matthews et al., 2016;Matthews, 2017].
There are many possible interdomain inducing variables that can be used to construct approximate posteriors. In certain cases, an appropriate choice can make a crucial difference to the effectiveness of an approximation, particularly when interdomain variables can be tailored to a particular kernel of interest. Álvarez et al. [2010] construct a kernel through a deterministic transformation of a white noise process, and want to create an approximation by learning about the white noise process. Normal inducing points are completely ineffective at representing knowledge about white noise, but well-chosen interdomain variables can capture the useful low-frequency characteristics of a realisation. Similarly, computationally tractable convolutional GPs rely on tailored interdomain variables . Well-chosen interdomain variables can also be used to induce computationally favourable structure in K uu for general Matérn kernels .
In the following section we introduce multioutput GPs and discuss how interdomain approximations allow such models to scale to large numbers of data points as well as large numbers of outputs. In section 5 we then introduce our software framework that allows kernels and interdomain variables to be specified in a flexible manner, which allows all combinations to be specified with minimal code duplication.

Multioutput Gaussian processes
So far, we have discussed approximations for scalar GP models. Multioutput Gaussian processes (MOGPs) (see Álvarez et al. [2012] for a review) are powerful models that can represent useful correlations between related outputs. Taking these correlations into account allows the information from one output to be used in the prediction of another. This helps in multi-task learning, when predicting multiple outputs using a single model. Convolutional layers in deep convolutional Gaussian process models [Blomqvist et al., 2018;Dutordoir et al., 2020] also exhibit multioutput structure, where outputs for patches in different image locations are correlated. While the additional correlations are helpful for the model, they present an even larger computational challenge than single-output models. In this section, we will introduce MOGPs and discuss how we can do efficient inference in the variational inference framework, leading up to the discussion of our software framework in section 5.
We are interested in models that learn a multioutput function f (·) : X → R P . We take the input space X to be R D , but much carries over to the general case. We denote by f p (x) the pth output of f (·) at the input x. A prior over such functions is a multioutput Gaussian process [Álvarez et al., 2012] if the distribution of a vector of values is Gaussian distributed. Note that for each element we specify the input x n as well as which output p n we want to consider. As with single-output GPs, multioutput GPs are fully specified by their kernel (we assume zero mean functions). Kernels for multioutput GPs can be seen as matrix-valued (or operator-valued in general) where a covariance matrix for all outputs is returned [Micchelli and Pontil, 2005]. This view has allowed in-depth analysis of kernel methods for learning vector-and even function-valued functions [Kadri et al., 2016]. An alternative view is that the matrix-valued kernel obeys the same positive definiteness properties of a single-output kernel where the index of the output is treated simply as another input. In other words, we can think of multioutput kernels as functions on the input space X extended by the index of the output: This means that the collection on random variables f from eq. 28 will have the distribution This insight shows that when we observe an arbitrary collection of different output dimensions of f (·) for different inputs, which we will refer to as heterotopic data, we can simply use the same procedure and software as for single-output GPs. While this view is conceptually simple, it obscures structure that can be taken advantage of for convenient and efficient implementation.
The "output as an input" view is particularly ill-suited for dealing with data where we observe (nearly) all output dimensions of f (·) for each input (homotopic data), both in terms of mathematical notation and software implementation. In these situations, we will instead use to denote the stacked outputs for all N inputs x n in X. We can obtain the covariance of p(f (X)) by taking the (n, n )th P × P block of K to be k(x, x ), using the matrix-valued kernel of eq. 29: We will sometimes index into f (X) and K as multidimensional arrays, in their unstacked form. Specifically, we denote The true index given n and p can be found to be i = (n − 1)P + p. This notation provides several benefits. Firstly, when dealing with all outputs, it is inconvenient to explicitly denote the output dimensions as inputs. Secondly, always considering all outputs may introduce structure in K, which can be computationally exploited (we will see examples later). Finally, in software terms, we want to organise a collection of vector-valued outputs and their covariances in high-dimensional arrays, which makes their grouping explicit, which simplifies later use.
We will now discuss several examples of multioutput kernels. In order to take advantage of details that help efficient implementation, we take the view of evaluating a vector-valued function at particular inputs, rather than the "output as an input" view. This will then lead to a software framework which a) allows easy specification of the most efficient implementation, and b) does not compromise on unified software interfaces. Álvarez et al. [2012] discussed different choices for priors on the correlations between multiple outputs. Here we briefly review these as well as convolutional GPs.

Linear model of coregionalization
A simple way of introducing correlations in the outputs is to construct our multioutput function f (·) from a linear transformation W ∈ R P ×L of L independent functions g (·) as: with f (x) ∈ R P and g(x) ∈ R L . This is known as the Linear Model of Coregionalization (LMC) [Journel and Huijbregts, 1978], and implies the covariance This construction can be viewed as a form of input-dependent PCA or Factor Analysis, where we have a linear relationship between a latent and observed space. We can choose L < P to obtain a low-dimensional explanation.

Convolution processes
The LMC only allows for a limited kind of correlations across outputs. It cannot capture, for example, outputs which are simply delayed versions of each other, i.e. f p (x) = f p (x + δ). Such more complicated relationships between outputs can still be expressed linearly by constructing f p (x) through a convolution of some latent process g(·). The resulting convolution process (CP) can exhibit relationships like time-lags, and more general linear dependence on past observations. Such models have been investigated in different forms over the years [Higdon, 2002;Boyle and Frean, 2005;. We follow the description by Álvarez et al. [2010], which constructs f (·) from a convolution of g(·) as which, when taking the prior on g(·) as before, results in the covariance G(·) is usually parameterised in a way that makes the integral tractable, and adds a number of parameters that balances flexibility against susceptibility to overfitting.

Image convolutional Gaussian processes
Gaussian processes with convolutional structure  have been proposed for learning functions on images. They work by learning functions on a patch-by-patch basis, after which all patch responses are combined to form a single response for the image. This is similar to convolutional neural networks, with an exact correspondence in a particular infinite limit [van der Wilk, 2019]. The construction starts with a patch response function g(·) : R w·h → R, operating on image patches of size w × h. An image response is obtained by summing the response to all P patches where x [p] is the pth patch in the image x, and P = (W − w + 1) × (H − h + 1) is the number of patches for a W ×H image. The construction can be phrased similarly to a convolution processes (section 4.1.2) by taking G(x, z) = P p=1 δ x [p] − z , and integrating over the space of patches rather than images.
This construction is easily adapted to create an image-to-image MOGP, by simply returning all patch responses, rather than summing them. We obtain the resulting function f (x) : R D → R P which has one output per patch. This construction was discussed by Dutordoir et al. [2020] to create a convolutional layer for use in a deep convolutional Gaussian process model. Since we construct all outputs using the same patch response function g(·), the correlation of all outputs is given by the covariance function of g(·):

Approximate inference
Exact inference with MOGPs is complicated by the same issues that were discussed in section 2.2. The computational scaling in particular is worse, since if we observe P outputs for all N inputs, we will obtain a P N × P N covariance matrix, resulting in O(N 3 P 3 ) time scaling. In the same way as for single-output GPs, these problems can be addressed in the variational framework (section 2.3).
In the following, we will consider general likelihoods for MOGP priors of the form p(y n | f (x n )), where the vector-valued observation y n can depend on all outputs of f (·) for the input x n . A simple but useful example to keep in mind is where we use a correlated Gaussian as the likelihood, e.g.: To find the variational bound for MOGP models, we have to take this modified likelihood into account, as well as any changes that come from using multioutput variational posteriors. Luckily, much of the single-output case carries through. As before, we define the approximate posterior through the prior conditioned on the inducing observations u (i.e. p(f (·) | u)), together with a free-from density q(u): Denoting the prior conditonal p(f (·) | u) is similar to the single-output case, but requires slightly more sophisticated notation. To obtain a similar expression, we will make the following conventions: • In multioutput scenarios, we will useM to refer to the number of elements in u, i.e. u ∈ RM . This is an important distinction for the interdomain inducing variables.
• We take the kernel k(x, x ) to be a matrix-valued kernel, as in eq. 29, i.e. k(x, x ) ∈ R P ×P .
• The covariance matrices K uf (x) and K uf (X) are obtained by the same stacking as eq. 34. This makes K uf (x) ∈ RM ×P , and K uf (X) ∈ RM ×N P .
• We index into these matrices in the same way as in eq. 35.
These conventions allow us to write the approximate posterior at x as We can now write the multioutput variational lower bound with a variational expectation over the P -dimensional q(f (x)) to take the multivariate nature of the new likelihood into account: This bound is written for homotopic data, as it assumes that the likelihood depends on all P dimensions of f (x). Heterotopic can be dealt with in the same way by simply marginalising out unobserved variables. For example, if we only observe a subset of outputs described by the index set S n at the input x n with the correlated Gaussian likelihood from eq. 42, we can integrate out all elements in y n that are not in S n . The likelihood will then also only depend on a subset of the outputs in S n . This only changes the variational expectation in the ELBO: Integrating out the unused random variables is done by selecting the relevant rows and columns of K T uf and Σ.
In the discussion above, we have used the vector of inducing outputs u without defining exactly which random variables they correspond to. This general discussion only assumes that we have a covariance of u as given by K uu . Creating efficient multioutput methods depends strongly on the choice of u, as this can induce structure in K uu which can make it easier to invert. In the next section, we will discuss various choices of u and their computational benefits.

Example: Linear Model of Coregionalization
In the previous section, we stated the variational bound for multioutput GP models, and that its computational efficiency depended on the structure in K uu . Here, we discuss variants of sparse approximations to the Linear Model of Coregionalization (section 4.1.1). The main design choice in the approximation is the definition of u, which determines the matrices K uu and K uf . We will show how different choices lead to approximations with different computational characteristics, both in terms of the ease with which K uu can be inverted, as well as how many kernel elements need to be computed. This section discusses examples in mathematical terms, while section 6.5 delves into the implementation.

Inducing points
We start with by introducing two formulations of inducing points for MOGPs. In the same way as in single-output GPs, these inducing points are general, and work with every multioutput kernel.
We can choose our inducing variables in an unstructured way, by specifyingM input-output pairs {z m , p m } using the "output as an input" formulation. This gives the following u and corresponding covariance matrices: Alternatively, we can define inducing points for MOGPs as specifying M inducing inputs, and using all corresponding outputs as inducing variables: In this case, we haveM = M P inducing variables for M inducing inputs. This "vector-valued inducing point" approach defines more inducing variables for fewer inducing inputs, which reduces the number of parameters to be optimised. This may be beneficial even when dealing with heterotopic data, particularly if the input distribution is the same for all outputs. Both these methods yield dense K uu and K uf (x) matrices for which no additional simplifying assumptions can be made. However, this method does immediately work for all multioutput kernels. We include the vector-valued inducing point approach in GPflow (section 6.5.1) as the most basic and general method.

Latent inducing points
For the LMC, approximate inference can be made more efficient by a more judicious choice of inducing variables. In section 4.1.1, we constructed the LMC from L independent latent processes g(·). If we choose evaluations of g(·) as our inducing variables, we can take advantage of the prior independence of g (·) over . We collect all outputs of g(·) for M inputs in Z, to obtainM = M L inducing variables. Thanks to the prior independence, K uu will have a block diagonal structure, meaning that the largest matrix that has to be inverted only has size M × M . We specify the approximation as We impose an additional mean-field restriction on q(u) by choosing S = blockdiag L =1 [S ], where S ∈ R M ×M . This gains extra computational efficiency at the cost of some accuracy. We can now rewrite eq. 44 to take advantage of this block-diagonal structure. We write K −1 uu for the inverse of the th M × M block of K uu .
where we implicitly sum over repeated indices (Einstein notation). Furthermore, by noting that k(x, x) pp = W p k (x, x)W p , we can write eqs. 51 and 52 as a linear transformation of the predictive distributions of g(·). If we denote the predictive density at the input x as q(g(x)) = N g(x), µ g , Σ g (note that Σ g is diagonal), we have

Latent inducing points for the Intrinsic Model of Coregionalisation
We can further gain computational benefits by constraining the kernels of each g (·) to be equal, which is known as the Intrinsic Model of Coregionalisation (IMC). In terms of modelling quality, this constraint may be beneficial or a hindrance depending on the characteristics of the dataset. However, it is much faster than the more flexible model of the previous section: instead of L separate inversions, only a single M × M inversion has to be performed. Our predictive mean and covariance become We can again rewrite these terms as in eq. 53, only in this case the computation of µ g and Σ g would require only a single inversion.

Computational complexity
The choice of inducing variable has a large impact on the computational complexity of the method. We compare the methods of the previous three sections in table 1. For the methods in sections 4.3.2 and 4.3.3, we have two ways to compute the required quantities: 1) first calculate K uf (X) followed by a summation over mm p, 2) first sum over mm to obtain µ g and Σ g , and then sum over . While method 2 is more efficient than method 1, this choice influences the complexity of implementation, which we discuss in section 6.5. In this section, we will present our unified implementation of interdomain and multioutput Gaussian processes, following the mathematical derivation that we described so far. Our framework is implemented as a large extension to GPflow . We start with describing the desiderata which guided our design choices, followed by a description of the code architecture. We will specifically refer to locations in the code where key items are implemented, so the reader can gain an idea for how other kernels and interdomain variables can be implemented.

Desiderata
Our contribution to GPflow maintains its limited scope of implementing only variational approximations to GP models. While many of the abstractions described here are useful for other methods (e.g. EP), we maintain the narrow focus to ensure high-quality code. Even with this limitation, the previous sections illustrated that a wide variety of models can be approximated with a wide variety of approximate posteriors. The goal of our software framework is to facilitate the implementation of as many combinations of model and approximation structure as possible, in a computationally efficient manner. We summarise this goal using three desiderata: 1) Efficiency: Our framework must allow the most computationally efficient code to be implemented, by exploiting any mathematical structure in the model and approximation.
2) Modularity: Using different combinations of modelling and approximation assumptions should require a minimum of effort and code-duplication.
3) Extensibility: Our framework should allow new model and approximation variations to be implemented without modifying core GPflow source, following the open/closed design principle [Martin, 1996].
The mathematical framework discussed in the previous sections requires us to make three choices that influence the overall method. The kernel and likelihood determine model properties, while the choice of interdomain variables u determines properties of the approximation. We will achieve modularity and extensibility by expressing the design choices as class hierarchies, which specify the most efficient code path through multiple dispatch.

GPflow model setup and examples
We focus our implementation around the variational bound from eqs. 16 and 45 [Hensman et al., 2013;Hensman et al., 2015b], as it works for general kernels, inducing variables, and likelihoods, in contrast to more specialised methods like Titsias [2009]. Our contribution extends the existing SVGP (Sparse Variational GP) implementation in GPflow . SVGP requires that three main components are specified: 1. the kernel, which determines the prior on f (·), 2. the inducing variables u, and 3. the likelihood.
Each is passed as an object to the SVGP, presenting a flexible interface for users to combine different kernels, inducing variables, and likelihoods. A model can be trained to the dataset (X, Y) by minimising its objective function with respect to its variational and hyperparameters, using: We now review several examples that show how our interface can be used, and that it allows for the creation of models that flexibly mix and match different components.

Single-output and inducing points
The most straightforward use case is simple single-output regression with inducing points on a dataset where X.shape == (N, D) and Y.shape == (N, 1). This model can be set up with listing 1. 1 kernel = gpflow . kernels . S q u a r e d E x p o n e n t i a l () 2 Z = X [: M , :]. copy () # Initialise inducing locations to the first M inputs in the dataset 3 i n duc ing _va ria bl e = gpflow . i n d u c i n g _ v a r i a b l e s . Induci ngPoin ts ( Z ) 4 likelihood = gpflow . likelihoods . Gaussian () 5 model = gpflow . models . SVGP ( kernel , likelihood , i n d u c i n g _ v a r i a b l e ) Listing 1: Set-up code for a single-output GP model, with an inducing point posterior.
Single-output and Fourier features The properties of the inducing variables greatly influence the computational properties of the method. Variational Fourier Features , for example, reduce the cost of inverting K uu , and allow precomputation of many costly quantities. We can use them by simply changing the inducing variable that is passed into SVGP in listing 1: where a set of N images is represented as a matrix X of shape N × (W · H). While using "default" InducingPoints is possible, van der Wilk et al. [2017] argue that tailoring the inducing variables to the convolutional structure of the kernel leads to a much more efficient method. We achieve this by instead using InducingPatches as inducing variables: A fully-working example can be found in the GPflow documentation [2019a].
Multioutput kernels Multioutput models also fit into this framework, and can be specified using a multioutput kernel. The usual InducingPoints can still be used to specify a posterior, as described in section 4.3.1, resulting in the code: 1 kernel_list = [ 2 gpflow . kernels . S q u a r e d E x p o n e n t i a l ( D ) for _ in range ( L ) 3 ] # specify a kernel for each output 4 W = np . random . randn (P , L ) # Linear mixing of the GP outputs 5 kernel = gpflow . kernels . L i n e a r C o r e g i o n a l i s a t i o n ( kernel_list , W = W ) 6 iv = gpflow . i nd u c i n g _ v a r i a b l e s . In ducin gPoint s ( Z ) Alternatively, we can use inducing points defined in g(·) to take advantage of the most efficient code path (section 4.3.2): 1 from gpflow . in d u c i n g _ v a r i a b l e s import InducingPoints , S h a r e d I n d e p e n d e n t I n d u c i n g V a r i a b l e s 2 iv = S h a r e d I n d e p e n d e n t I n d u c i n g V a r i a b l e s ( Induci ngPoin ts ( Z ) ) See the GPflow documentation [2019b] for more examples on multi-output GPs in GPflow.

Implementation of SVGP
The main role of SVGP is to specify the computation of the training objective, and to allow for future predictions once the model is trained. As stated earlier, we train our model by maximising the ELBO (eqs. 16 and 45), which we reproduce here: We need to implement SVGP and its components in a way that satisfies our desiderata of efficiency, modularity, and extensibility. The example use-cases in the previous section determine the interface of SVGP to ensure modularity. To satisfy extensibility, we need to ensure that SVGP does not rely on properties that are specific to a particular kernel or inducing variable, and that new behaviour can be implemented without modifications to GPflow. We will separately discuss each of the four components that form the ELBO: 1) the KL term, 2) the predictive distribution q(f (x)), 3) the expectation of the log-likelihood over the predictive distribution, and 4) the sum of the expectations over the log-likelihood. Although slightly hidden from the implementation of SVGP, we will also discuss how the covariance matrices K uu and K uf are computed, which the KL term and predictive distribution rely on.
We see a skeleton of the implementation of SVGP in listing 2. The interface of the constructor permits the usage described in the previous section, by accepting a kernel, inducing variable, and likelihood. Optionally, a user can specify a mean function, and a different parameterisation of the variational parameters (using whiten, see section 5.3.3). The core of SVGP computes the objective function in SVGP.elbo(), where we see the separation of the different terms of the ELBO: • KL divergence The first quantity we compute in elbo() is the KL[q(u) p(u)] term (eq. 56), i.e. the KL divergence between the approximate posterior and prior of the inducing outputs. This can be calculated in closed form as both distributions are Gaussian. The method that computes this quantity, kullback_leiblers.prior_kl() takes the inducing variable and kernel objects, in addition to the specification of the variational posterior through q_mu and q_sqrt. We pass the inducing variable and kernel in as objects so we can specialise the computational path to take advantage of any resulting structure in K uu .
• Predictive distribution The predictive distribution q(f (x)) is computed through a call to the function conditional() inside predict_f(). The method will return the mean fmean and the variance fvar of the posterior predictive GP evaluated at points Xnew, following the general equation of a conditional GP eq. 13. As with the KL divergence, we want to tailor the code path to the combination of inducing variable and kernel that is being used, so we pass these objects in to the function as parameters.
• Expected log-likelihood The expected log-likelihood is encapsulated in classes inheriting from Likelihood, which only requires the predictive distribution (fmean and fvar), and is therefore independent of any model choices once conditional() has been called.
• Minibatch sum In the last line of elbo() we compute the sum of the expected loglikelihoods of the data points in the batch and rescale its contribution to the overall ELBO to compensate for the minibatching. Listing 2: GPflow implementation of SVGP.
We will now delve into the implementation of these terms in more detail.

Multiple dispatch
Both for computing the KL divergence (prior_kl()) and for computing the predictive density (conditional()), the code path depends on the type of inducing variables and kernel we use. One possible approach for customising the code paths inside these functions would be to have branches for each pair that is implemented, e.g.: While this implementation would allow the most efficient code paths to be specified, it is not extensible. We want to be able to implement specific code for new combinations of kernel and inducing variable without having to edit core GPflow code.
To provide this functionality, we rely on multiple dispatch (we use the Python implementation by Rocklin [2014]). When a multiple dispatch function is called, the runtime types of its arguments are used to determine which specific code path to run. This is a runtime equivalent of function overloading in statically typed programming languages. Multiple dispatch allows new code paths to be registered as separate functions, which can be placed anywhere in code, therefore enabling extensibility. A code path for a new pair of types is registered using the dispatcher.register(*types) decorator. For example, in GPflow we register various implementations of conditional(): The covariance matrices used to calculate the conditional have the following shape : """ Most efficient routine to project L independent latent gps through a mixing matrix W .

25
The mixing matrix is a member of the ' LinearCoregionalization ' and has shape [P , L ].

26
The covariance matrices used to calculate the conditional have the following shape : When conditional() is called, the appropriately typed implementation from above is run. The inheritance structure of the inducing variable and kernel classes is respected, with the most specific type that has an implementation being preferred, and parent types providing default implementations. For example, if a new inducing variable NewInducingVariable inherits from InducingPoints, conditional() will call the second implementation in listing 3, unless a more specific implementation is registered with @conditional.register(NewInducingVariable, kernels. Kernel).
The key benefits of implementing SVGP in this way, is that the model is completely general to any kernel, inducing variable, and likelihood. Any code that is specific to a particular combination (e.g. for efficiency) can still be specified, but external to the model and even the GPflow package itself (e.g. GPflow documentation [2019c]). This achieves our desiderata of efficiency, and extensibility. In the next sections, we will describe the details of each component of the ELBO: covariances (section 5.3.2), prior_kl() (section 5.3.3), conditional() (section 5.3.4), and variational_expectations() (section 5.3.5).

Covariances K uu & K uf
Although not directly visible in SVGP, both prior_kl() and conditional() need to compute the covariances K uu and K uf . Since their expressions (eq. 22) depend on both the type of the kernel and inducing variables, we also implement their interface with multiple dispatch (listing 4).   Listing 4: Three definitions for interdomain covariances, using multiple dispatch for extensibility.

Kullback-Leibler divergences
To compute the KL divergence between the approximate posterior process q(f (·)) and the prior process p(f (·)), we only need to compute the finite-dimensional KL divergence at the inducing points [Matthews et al., 2016]: Due to the possible structure in K uu , we use multiple dispatch on the inducing variables and kernel. This gives a general interface for prior_kl() of: 1 def prior_kl ( inducing_variables , kernel , q_mu , q_sqrt , whiten = False ) : 2 ...
For unstructured K uu , GPflow provides a default gauss_kl() function to compute eq. 58. In the formulation above, it requires the K uu matrix, which is computed with the multiple dispatch function Kuu, together with the variational parameters m and S. We parameterise S through its lower-triangular Cholesky factorisation, i.e. S = L S L T S , which makes log |S| = 2 M m=1 log[L S ] mm , which is computationally cheaper.
We can alternatively compute the KL without direct access to K uu by whitening the inducing variables. In SVGP this is indicated by setting whiten=True. Instead of conditioning on u, we use  a transformed version v designed to have an identity covariance in the prior: where K uu = L uu L T uu . It is clear that this construction implies the correct distribution p(u) = N (0, K uu ). If we represent the variational posterior over v instead of u, then Since the prior over v is now a centered standard Normal, this simplifies the KL divergence: Whitening can make optimisation or MCMC sampling more efficient [Hensman et al., 2015a]. In the GPflow source code a whiten flag indicates whether we are passing (m v , S v ) or (m, S).

Conditionals
The computation of the approximate posterior distribution q(f (X)) is the most complex component to be specified. It is required for computing the expected log-likelihood (section 5.3.5), as well as for making predictions. In its most complete form, q(f (X)) is specified by covariances between every output for every input, giving N × P × N × P elements in total. Computing the expected log-likelihood never requires covariances between outputs for different inputs, and if the likelihood factorises over the output dimensions, only the P marginal variances are needed. Covariances between outputs for different inputs are generally only needed for plotting. This gives four different predictive covariances that could be computed, determined by whether covariances are needed for different inputs (full_cov={True|False}), and outputs (full_output_cov={True|False}). We bring all this functionality into the single function conditional(), which deals with computing q(f (x)) with all four output covariance shapes. From the point of use in SVGP, as described in section 5.2, we need a fixed signature for conditional(), which we specify as: The return values of conditional() are the mean and covariance of function values for the N input points in Xnew. The return shape of the covariance depends on the value of full_cov and full_output_cov, summarised in table 2. We choose the layout of the returned covariances based on what is most convenient for follow-on operations. Since these are usually matrix multiplications or Cholesky decompositions, the tailing dimensions always contain the matrix of covariances that are requested. In principle, conditional() only needs to implement eq. 44. However, different combinations of priors and inducing variables induce different structure in K uu and K uf , which can be exploited for computational efficiency. As described earlier, we will use multiple dispatch to allow the most efficient code path to be specified, while keeping a consistent interface. We dispatch on the specific inducing_variable and kernel that are passed to conditional(). In the following sections, we describe three different conditionals that are implemented in GPflow, each making a different assumption about K uu and K uf . In section 6, we discuss various models (including the LMC from section 4.3) that can be implemented efficiently using these and other custom conditionals.

Likelihood
The expected log-likelihood is computed by Likelihood.variational_expectations(Fmu, Fvar, Y), which implements the integral This function is separated from all other components as it depends only on the variational predictive distribution q(f (x n )). This is computed by conditional(), and represented through its mean µ n ∈ R P (passed as Fmu) and covariance Σ n ∈ R P ×P (passed as Fvar).
For certain likelihoods this expectation can be calculated in closed form. This is implemented by inheriting from the Likelihood base class, and implementing the closed-form expression in Likelihood.variational_expectations(Fmu, Fvar, Y). If this not possible, GPflow can fall back to either Gauss-Hermite quadrature or Monte Carlo estimation. If Σ n is diagonal (i.e. the outputs under the variational posterior are independent), or if the likelihood is independent over dimensions as log p(y n | f (x n )) = p log p(y np | f p (x n )), the multidimensional integral becomes equivalent to a sum over single dimensional integrals, which makes Gauss-Hermite or Monte Carlo estimation particularly effective. These can be used by inheriting from the classes Likelihood or MonteCarloLikelihood, respectively.

Implementation of example methods in our framework
In this section, we wish to demonstrate the efficiency, modularity, and extensibility of our framework by discussing several example methods. We pay particular attention to demonstrating the use of multiple dispatch, and show how the appropriate Kuu(), Kuf(), and conditional() functions can be implemented.

Single-output Gaussian processes
We begin with the most simple example of a single-output GP with an approximate posterior that induces no structure in K uu . Listing 1 sets up an example of such a model with an inducing point posterior, with Kuu() and Kuf() being defined in listing 4. To complete the implementation we need to define a version of conditional() and register it using @conditional.register(...). Since there is no structure in K uu , we need to define a version of conditional that simply implements eq. 17 in accordance with the shapes in table 2. We register the function using the following code: We register the function as a conditional with the dispatcher for general Kernel and Inducing Variable objects, using @conditional.register(...). This will make this code path the default for all kernels and inducing variables that inherit from the base classes Kernel and Inducing Variable, unless more specialised implementations are provided. The main work of implementing eq. 17 is done in base_conditional(). We pass in K uu and K uf as matrices, after they have been evaluated using the multiple dispatched functions Kuu() and Kuf(). To keep base_conditional() simple, it does not implement full_output_cov=True. To ensure that the returned shapes variances are consistent with other multioutput conditionals, we use the helper function _expand _independent_outputs().

Multiscale interdomain inducing variables
In section 3 we discussed the use of integral transformations of a function f (·) as inducing variables. Lázaro-Gredilla and Figueiras-Vidal [2009] proposed using w m (x) being a Gaussian, or a Gaussianwindowed sinusoid. When combined with the SquaredExponential kernel, the integrals for the required covariances (eqs. 23 and 24) have closed form expressions. We included the Gaussian integral transformation as the Multiscale inducing variable in GPflow. While these inducing variables do not provide special structure in K uu and K uf , they do allow the posterior to represent longer-range correlations with fewer inducing variables. In terms of code, this method can be implemented by simply specifying the appropriate covariances: Using the @method.register(...) line, we ensure that this code is specific to the combination of the Multiscale inducing variable, and the SquaredExponential kernel. Since Multiscale and SquaredExponential inherit from InducingPoints and Kernel respectively, the same version of conditional(...) will be called as in the previous section.

Image convolutional Gaussian processes
Similarly to Multiscale inducing variables, the inducing patches for the convolutional Gaussian process  can be implemented using by only specifying Kuu() and Kuf(), while relying on the default conditional. We define InducingPatch for specifying inducing variables in g(·) (see section 4. Image-to-image convolutional Gaussian processes [Dutordoir et al., 2020] also fit within this framework, but require custom conditionals to implement the covariances between the different outputs.  define inducing variables in the spectral domain of the process which lead to periodic functions for K uf . More importantly, this careful choice leads to a "low-rank plus diagonal" form of K uu , i.e.

Variational Fourier features
where α ∈ R M and β ∈ R M ×d , where d − 1 is the mean-squared differentiability of the functions in the prior. This method naturally fits into GPflow's interdomain and multiple dispatch framework. We can define a new type of feature, FourierFeatures, which inherits from InducingVariable. It would be mathematically possible and an immediate fit into the existing framework if we implement a new Kuu() method which computes the full (non-structured) K uu , and register it within the multiple-dispatch framework.
This approach, however, is sub-optimal and would not exploit any of the known structure (diagonal + low rank) of K uu . To take advantage of this, we need to route the model through more specialised code for both the KL divergence and the conditional. We let the covariance method for K uu return an object that captures the structure of diag(α) + ββ T (in this case, a TensorFlow LinearOperator). We then register a specialised conditional method that can deal with subsequent computations in the computationally most efficient way given its access to diag(α) and β.
The covariances can be implemented closely following the equations in Hensman et al., 2017, e.g. for the Matérn 1/2 kernel: Note that different orders of Matérn kernels have different structures in their K uu matrix, so we have to separately register one implementation for each combination of FourierFeatures inducing variable and one of the Matérn kernels. Similarly, we need to register separate implementations for K uf . The conditional, however, only needs to know that the solve operation is efficient. A sample implementation is given in the GPflow documentation [2019c].

Multioutput inducing points
When moving to specifying approximate posteriors for multioutput kernels, we have to deal with all the possible covariance structures described in table 2. We start with the inducing point method for multioutput GPs, which was described in section 4.3.1 and applied to the LMC. This method takes allM = M P output values from M inputs as inducing variables. The single-output conditional (section 6.1) does not compute the required covariances for full_output_cov=True, so we will need define a new conditional().
To implement the predictive equation (eq. 44), Kuu() needs to compute all M × P × M × P covariances, while Kuf() needs to compute the M × P × N × P covariances between u and the predictions. We let these functions simply compute and return the covariances defined in the kernel. Since this implementation of the InducingPoints inducing variable is general to all multioutput kernels, we register the implementation for all kernels that inherit from the multioutput kernel base class MultioutputKernel: We now require an implementation of conditional() that can use these differently shaped covariances. When full_cov and full_output_cov are equal, i.e. we only care about the marginals, or we want the complete joint over both the inputs and outputs, we can rely on the base_conditional() that we used earlier after a simple reshape. For the other cases, we implement fully_correlated_conditional(). As with the implementation of the covariances, we define this conditional() to be the default for all multioutput kernels when InducingPoints are used. For this example, we can assume that P is equal to L. In order to use this code with the LMC, all we have to do is define the covariance of the prior in an object that inherits from the multioutput base class MultioutputKernel:

Latent inducing points
In section 4.3.2, we discussed that placing the inducing variables in the process g(·) reduced the computational cost of the method.
The previous implementation focused on generality. Here we want to take advantage of additional structure. The IndependentLatent kernel allows inference with inducing variables Fallback {Shared|Separate}InducingVariables. These inducing variables give a block-diagonal structure in K uu , which is represented as a L × M × M array (section 4.3.2). To take advantage of this we define a custom conditional for these combinations of inducing variables and kernels. The exact details are abstracted away in independent_interdomain_conditional(*) which takes advantage of the block-diagonal structure. The previous implementation of LMC using the default MultioutputKernel and InducingPoints leads to a scaling of O P 3 M 3 + N P 3 M 2 . In this section, using the specialised codepath gives O LM 3 + N P LM 2 (see section 4.3.4).

Background
So far we have discussed a set of models for which all inputs {x n } N n=1 are known and fixed variables: they form the features of the model. This is the typical setting for most supervised learning algorithms. In certain applications, however, the GP models have to deal with inputs which are uncertain, and typically described by a distribution {p(x n )} N n=1 , rather than a single point. A wellknown example of this type of model is the GP-LVM [Lawrence, 2004], where we are interested in learning the joint posterior of two variables: the distribution of the input locations {q n (x n )} N n=1 of the GP function and the GP mapping f (·) itself. More recent models such as the deep Gaussian process (DGP) [Damianou and Lawrence, 2013] and the conditional density estimation model [Dutordoir et al., 2018] also require dealing with distributional or uncertain inputs.
Uncertain inputs are omnipresent in GP models and need to be accounted for in our software tool. In the following sections we discuss the optimisation bounds of the GP-LVM and DGP models, and show how they are implemented in GPflow. The focus of this section is to show that even when constructing these more complex models, we can still follow our desiderata of efficiency, modularity and extensibility. This leads to a software package that can be used to extend the functionality or modelling assumptions of DGPs by re-implementing only certain bits of the computational graph while staying inside the framework.

Latent Variable Model
The GP-LVM [Lawrence, 2004] is an unsupervised algorithm which learns for each of the given points in the dataset {y n } N n=1 an associated input distribution {q n (x n )} N n=1 , together with a GP mapping, so that y n = f (x n ) + noise. This is different from the typical regression or classification settings, where both the input and outputs are known and the only quantity to infer is the underlying mapping.
Deriving the bound for the GP-LVM is similar as in section 2.5: the approximate GP is defined by a set of inducing variables, which are learned by minimising the KL between the true and the approximate posterior. To learn the approximate posterior of the inputs in the GP-LVM we follow the same procedure, and define an approximate posterior over the uncertain inputs {q n (x n )} N n=1 . Following Bayes' rule, as in section 2.5 we obtain the lower bound (65) Most of the elements of this bound are already discussed in section 5.3, and stay unchanged in the case of uncertain inputs. However, the additional unknown variables and associated approximate posterior introduce two new components to the bound.
The first component is a KL term (for each datapoint) between the prior and approximate posterior of the uncertain inputs x n . As with the KL term on the inducing variables, given reasonable choices for the form of distributions this quantity can be computed in closed form.
The second component is an additional expectation of the data-fit term with respect to q(x). In general, solving this expectation analytically is not possible (as is the case for the inner expectation over q(f (x n )) for most likelihoods). To deal with this we resort to Monte Carlo and sample x n ∼ q(x n ) to estimate the quantity. Using re-parameterised samples [Kingma and Welling, 2013;Rezende et al., 2014] from the approximate posterior distribution allows us to get unbiased gradient estimates of the ELBO and optimise for the parameters of q(x), simultaneously with the other parameters as before.
For certain kernels and forms of q(x n ) the expectation can be computed analytically by making use of the so-called kernel expectations [Damianou et al., 2016]. While this approach adds less stochasticity to the bound, it increases the computational costs and limits the set of kernels that can be used.

Deep Gaussian processes
A deep Gaussian process consists of multiple GPs which are hierarchically chained to form a compositional function, f = f L • f L−1 • . . . • f 1 . At each level the input of the current function is given by the output of the previous one. The uncertain inputs in a DGP arise from the fact that the hidden state that is used as the input to the next GP is given by the output of a GP, which is itself a distribution.
As was the case for shallow GP models (see section 2.3 and section 4.2), we construct an approximate posterior GP by conditioning the prior GP on a set of inducing variables. The inducing variables u , for which we have a set for each layer, are then integrated out with respect to their posterior q(m , S ), leading to an approximate posterior of the form N (µ (h −1 ), Σ (h −1 , h −1 )) where the mean and the variance are given by The crucial difference between this approximate posterior and eq. 44 is that here the inputs to the mean and variance are the outputs of the previous layer h −1 , and thus uncertain variables. As was the case in the GP-LVM model, integrating out the inputs to account for their uncertainty is not possible due to the non-linear mapping of the inputs in the kernel. However, when appropriately propagating the uncertainty through the layers of the DGP, it is possible to compute a sample at each layer: Using the approximate posterior for every GP in the composition, we can lower-bound the log marginal likelihood to get the ELBO for a DGP: In this bound, compared to the shallow GP case (eq. 16), we are dealing with a variational expectation over the output of the last layer q(h L ) and a sum over the KL terms of the inducing variables of every layer. Evaluating the KLs in this bound is straightforward given the Gaussianity of both prior and approximate posterior of u . The required computation is identical to the one required in the shallow case. The expectation over the final layer is approximated using Monte Carlo. Following the procedure of eq. 69, we obtain a sample from q(h L ) which is used to get an unbiased estimate of the variational expectation, as described in Salimbeni and Deisenroth [2017].

Implementation
Implementing bounds for models with uncertain inputs, e.g. the GP-LVM (eq. 65) and DGP (eq. 70) is not significantly harder, as most of the required components are identical to the building blocks of the SVGP model discussed in section 5.3. The most notable difference between the standard SVGP model and models with uncertain inputs is that for the SVGP, we compute the mean and variance of q(f (x n )) at known locations {x n } N n=1 , whereas for models with uncertain inputs, we evaluate the posterior GP at samples from the source of uncertainty. This can be the distribution of the inputs q(x n ) in the case of the GP-LVM or the output distribution of the previous layer q(h ) for DGPs.
To accommodate this use-case, GPflow provides a method called sample_conditional, whose role is to compute a sample from the approximate posterior q(f (x)). The method has the same signature as conditional and in most cases will use conditional to compute the mean µ and covariance Σ of the approximate posterior, so that it can construct a sample µ + Σ 1 2 , with ∼ N (0, I). However, for certain types of kernels and inducing variables it is possible to compute a valid sample from the posterior in a more efficient way, which is why sample_conditional itself makes use of multiple dispatch. As was the case for conditional (see section 5.3.4 and section 5.3.1), thanks to multiple dispatch, specialised code can be executed depending on the types of the arguments to the method in order to compute the sample in the most efficient way.

Modularity and extensibility
The sample_conditional method computes a sample from the GP's approximate posterior, which can then be propagated through the layers of a DGP to obtain a valid sample of the composite function. A user can alter the behaviour of a layer in a DGP by defining a new kernel and inducing variable, together with the corresponding sample_conditional. This makes GPflow easy to extend and modular, as evidenced by e.g. Blomqvist et al. [2018] and Dutordoir et al. [2020] who implemented deep convolutional (multi-output) GPs following this approach, and Salimbeni et al. [2019] who implemented a DGP with latent variable layers.

Conclusion
In this work we presented the current state of Gaussian process models. We reviewed current modelling assumptions such as multi-output GPs and discussed methods for efficient inference: sparse variational inference and inter-domain inducing variables. We presented this in a unified mathematical framework, starting from single-output models up to multi-output inter-domain GPs. Throughout the text it becomes clear that efficient inference in GP models relies strongly on choosing useful inducing variables in a way that can take computational advantage of independence properties of the prior. To this end, we developed a software framework in GPflow that allows for the flexible specification of both multi-output priors and inducing variables.