Modular Gaussian Processes for Transfer Learning
Year: 2021

Abstract:
We present a framework for transfer learning based on modular variational Gaussian processes (GP). We develop a module-based method that having a dictionary of well fitted GPs, one could build ensemble GP models without revisiting any data. Each model is characterised by its hyperparameters, pseudo-inputs and their corresponding posterior densities. Our method avoids undesired data centralisation, reduces rising computational costs and allows the transfer of learned uncertainty metrics after training. We exploit the augmentation of high-dimensional integral operators based on the Kullback-Leibler divergence between stochastic processes to introduce an efficient lower bound under all the sparse variational GPs, with different complexity and even likelihood distribution. The method is also valid for multi-output GPs, learning correlations a posteriori between independent modules. Extensive results illustrate the usability of our framework in large-scale and multitask experiments, also compared with the exact inference methods in the literature.

heading: 1 Introduction
text: Imagine a supervised learning problem, for instance regression, where N data points are processed for training a model. At a later time, new data are observed, this time corresponding to a binary classification task, that we know are generated by the same phenomena, e.g. using a different sensor.
Having kept the observations from regression stored, a common approach would be to use them in combination with the classification dataset to generate a new model. However, this practice might be inconvenient because of i) the need of centralising data to train the model, ii) the rising data-dependent computational cost as the number of samples increases and iii) the obsolescence of fitted models, whose future usability is not guaranteed for the new set of observations. Looking at the deployment in large-scale scenarios, by any organization or use case, where data are ever changing, this solution becomes prohibitive. The main challenge is to incorporate unseen tasks as former models are recurrently discarded.
Alternatively, we propose a framework based on modules of Gaussian processes (GP) (Rasmussen and Williams, 2006). Considering the current example, the regression model (or module) is kept intact. Once new data arrives, one fits a meta-GP using the module, but without revisiting any sample. If no data are observed, combining multiple modules is also allowed — see Figure 1. Under this framework, a new family of module-driven GP models emerges.
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
Background. Much recent attention has been paid to the notion of efficiency, e.g. in the way of accessing to data. The mere limitation of data availability forces learning algorithms to derive new capabilities, such as, distributing the data for federated learning (Smith et al., 2017), recurrently observe streaming samples for continual learning (Goodfellow et al., 2014) and limiting data exchange for private-owned models (Peterson et al., 2019). A common theme is the idea of model memorising and recycling, that is, using the already fitted parameters for an additional task without revisiting any data. If we look to the probabilistic view of this idea, uncertainty is harder to be repurposed than parameters of functions. This is the point where Gaussian processes play their role.
The flexible nature of GPs for defining a prior distribution over non-linear function spaces has made them a suitable alternative in probabilistic regression and classification. However, GP models are not immune to settings where we need to adapt to irregular ways of observing data, e.g. noncentralised observations, asynchronous samples or missing inputs. Such settings, together with the well-known computational cost of GPs, typically O(N3), has motivated plenty of works focused on parallelising inference. Before the modern era of GP approximations, two seminal works distributed the computational load with local experts (Jacobs et al., 1991; Hinton, 2002). While the Bayesian committee machine (BCM) of Tresp (2000) focused on merging independently trained GP regression models on data partitions, the infinite mixture of GP experts (Rasmussen and Ghahramani, 2002) combined local GPs to gain expressiveness. Moreover, the emergence of large datasets, with size N>104, led to the introduction of variational methods (Titsias, 2009a) for scaling up GP inference. Two recent works that combined sparse GPs with distributed regression ideas are Gal et al. (2014); Deisenroth and Ng (2015). Their approaches are focused in exact GP regression and the distribution of computational load, respectively.
Contribution. We present a framework based on modular Gaussian processes. Each module is considered to be a sparse variational GP, containing only variational parameters and hyperparameters. Our main contribution is to keep such modules intact, to build meta-GPs without accessing any past observation. This makes the computational cost minimal, becoming a competitor against inference methods for large-datasets. The key principle is the augmentation of integrals, that can be understood in terms of the Kullback-Leibler (KL) divergence between stochastic processes (Matthews et al., 2016). We experimentally provide evidence of the benefits in different applied problems, given several GP architectures. The framework is valid for multi-output settings and heterogeneous likelihoods, that is, handling one data-type per output. Our idea follows the spirit of not being data-driven any longer, but module-driven instead, where the Python user specifies a list of models as input: models = {model1, model2, . . . , modelK}. This will be later called the dictionary of modules.

heading: 2 Modular Gaussian processes
text: We consider supervised learning problems, where we have an input-output training dataset D = {xn, yn}Nn=1 with xn ∈ Rp. We assume i.i.d. outputs yn, that can be either continuous or discrete variables. For convenience, we will refer to the likelihood term p(y|θ) as p(y|f) where the generative parameters are linked via θ = f(x). We say that f(·) is a non-linear function drawn from a zeromean GP prior f ∼ GP(0, k(·, ·)), and k(·, ·) is the covariance function or kernel. Importantly, when non-Gaussian outputs are considered, the GP output function f(·) might need an extra deterministic mapping Φ(·) to be transformed to the appropriate parametric domain of θ. Data modules. The dataset D is assumed to be partitioned into an arbitrary number of K subsets or modules, that are observed and processed independently, that is, {D1, D2, . . . ,DK}. There is not any restriction on the number of modules, {Dk}Kk=1 do not need to have the same size, and we only restrict them to be Nk<N . However, since we deal with a huge number of observations, we still consider that Nk for all k ∈ {1, 2, . . . ,K} is sufficiently large for not accepting exact GP inference due to temporal and computational demands. Notice that k is an index while k(·, ·) is the kernel.

heading: 2.1 Sparse variational approximations for independent modules
text: For each data module, we adopt the sparse variational GP approach based on inducing variables, and the framework of Titsias (2009b), where posterior mean vectors and covariance matrices are computed as parameters. The use of auxiliary variables with approximate inference methods is widely known in the GP literature (Hensman et al., 2013, 2015; Matthews et al., 2016). In the context of K independent data modules, we define subsets of Mk Nk inducing inputs Zk = {zm}Mkm=1, where zm ∈ Rp.
Their non-linear function evaluations by f(·) are denoted as uk = [f(z1), f(z2), . . . , f(zMk)]>. As a first approach, we assume that f is stationary across all modules, being all uk ∀k ∈ {1, 2, . . . ,K} function evaluations of f(·). This will be later relaxed in the following section. To obtain the independent approximation to the posterior distribution p(f |Dk) of the GP, we introduce a kth variational distribution qk(f) for each data module Dk. In particular, this variational density factorises as qk(f) = p(f 6=uk |uk)qk(uk), with qk(uk) = N (uk|µk,Sk) and p(f6=uk |uk) being the standard conditional GP prior distribution given the hyperparameters ψk of the kth kernel. To fit the GP modules, and their variational distributions qk(uk), we build lower bounds Lk on the marginal log-likelihood (ELBO) of every dataset Dk. Then, we use gradient-based optimisation methods to maximise the K objective functions Lk, one per module, in an separate asynchronous manner. The ELBOs, for either regression or classification tasks are obtained as follows
Lk = Nk∑ n=1 Eqk(fn) [log p(yn|fn)]− KL[qk(uk)||pk(uk)], (1)
with pk(uk) = N (uk|0,Kkk), where Kkk ∈ RMk×Mk has entries k(zm, zm′) with zm, zm′ ∈ Zk and conditioned to certain kernel hyperparameters ψk that we also aim to optimise. The variable fn corresponds to f(xn) and the marginal posterior comes from qk(fn) = ∫ p(fn|uk)qk(uk)duk. In practice, the distributed bounds Lk are identical to the one presented in Hensman et al. (2015) and can be combined with stochastic variational inference (Hoffman et al., 2013; Hensman et al., 2013).
Dictionary of modules. The principal goal here is to obtain a dictionary, containing the already fitted GP modules, for their later use. Such dictionary consists, for instance, of a list of module objects {M1,M2, . . . ,MK} without any specific order, where eachMk = {φk,ψk,Zk,L∗k}, being φk the corresponding variational parameters µk and Sk. Notice that we also include L∗k in addition to the fitted inducing points and hyperparameters. This is the ELBO obtained during optimization. It will play a key role later for the module-driven learning process.

heading: 2.2 Module-driven lower bound
text: Ideally, to obtain an inference solution given the GP modules included in the dictionary, the resulting posterior density should be valid for all data modules {Dk}Kk=1. This is only possible if we consider the entire dataset D in a maximum likelihood criterion, that is, using the global log-evidence. Specifically, our goal is now to obtain an approximate posterior distribution q(f) ≈ p(f |D), by maximising a lower bound LM under the log-marginal density log p(D). This bound should not revisit any data but only the objects in the dictionary of modules, that are models. Notice that the GP model will be no longer data-driven, but module-driven instead. This is, we build a new model from models, what we call a GP meta-model or meta-GP.
We begin by considering the full posterior distribution of the stochastic process, similarly as Burt et al. (2019) does for obtaining an upper bound on the Kullback-Leibler (KL) divergence. The idea is to use the large-dimensional integral operators introduced by Matthews et al. (2016) in the context of variational inference, and previously by Seeger (2002) for standard GP error bounds. The use of the large-dimensional integrals is equivalent to an augment-and-reduce strategy (Ruiz et al., 2018). It consists of two steps: i) we augment the model to be conditioned on the high-dimensional index set of the stochastic process and ii) we apply properties of Gaussian marginals to reduce the integral operators to a finite amount of GP function values of interest. Similar strategies have been previously used in continual learning for GPs (Bui et al., 2017; Moreno-Muñoz et al., 2019).
Global evidence objective. The construction considered is as follows. We first denote y as all the output targets {yn}Nn=1 in the dataset D and f+ as the augmented large-dimensional GP. Notice that f+ sufficiently large and it contains the function values taken by f(·) at both {xn}Nn=1 and {Zk}Kk=1 for all data modules. We consider it finite. The log-marginal expression is therefore
log p(y) = log p(y1,y2, . . . ,yK) = log ∫ p(y, f+)df+, (2)
where yk = {yn}Nkn=1 are the output data-points already used for training each GP moduleMk. The joint distribution in (2) factorises according to p(y|f+)p(f+), where the l.h.s. term is the augmented likelihood distribution and the r.h.s. term would correspond to the GP prior over the large index set of the stochastic process. Then, we introduce the global variational distribution q(u∗) = N (u∗|µ∗,S∗)
that we aim to fit by maximising a lower bound under log p(y) given the dictionary of modules. The variables u∗ correspond to function values of f(·) given the new subset of inducing inputs Z∗ = {zm}Mm=1, where M is the free-complexity degree of the global variational distribution. To derive the bound, we exploit the reparametrisation introduced by Gal et al. (2014) for distributing the computational load of the expectation term with GP latent variable models (GPLVM) (Lawrence, 2005). This is based on the decoupling of data-points conditioned on the inducing inputs u∗. Applying Jensen’s inequality, it is obtained as
log p(y) = log ∫∫ q(u∗)p(f+6=u∗ |u∗)p(y|f+) p(u∗)
q(u∗) df+6=u∗du∗ ≥ Eq(u∗) [ Ep(f+ 6=u∗ |u∗) [log p(y|f+)] + log p(u∗)
q(u∗)
] , (3)
where we applied the properties of Gaussian conditionals to factorise the GP prior as p(f+) = p(f+6=u∗ |u∗)p(u∗). Notice that f+ = {f+6=u∗ ,u∗}. Here, the last prior distribution is p(u∗) = N (u∗|0,K∗∗) where [K∗∗]m,n := k(zm, zn), with zm, zn ∈ Z∗. We also consider k(·, ·) conditioned to the global kernel hyperparameters ψ∗ that we also aim to estimate. The double expectation in (3) comes from the factorization of the large-dimension integral given f+ and the application of Jensen’s inequality twice. Its derivation is in the appendix.
Likelihood reconstruction. The augmented likelihood distribution p(y|f+) is perhaps the most important part of the derivation. It allows us to apply conditional independence (CI) between the data modules Dk, and particularly, the output variables yk. This gives a factorized term that we will later use for introducing the GP modules Mk stored in our dictionary, that is, log p(y|f+) = ∑K k=1 log p(yk|f+). To avoid revisiting each kth likelihood term, and hence, evaluating the corresponding data which might be unavailable, we want to estimate the likelihood distributions. In particular, we use Bayes theorem conditioned to the large-dimensional augmentation f+. Bayes theorem indicates that each kth likelihood can be approximated as
p(yk|f+) ≈ Zk qk(f+)
pk(f+) = Zk
(((( ((p(f+6=uk |uk)qk(uk) ((( (((p(f+6=uk |uk)pk(uk) = Zk qk(uk) pk(uk) , (4)
since qk(f+) ≈ p(f+|yk). Here, Zk refers to the normalization constant or evidence of Dk, which is still difficult to compute. We also assumed that the augmented approximate density factorises as qk(f+) = p(f+6=uk |uk)qk(uk), using the structure of the GP prior pk(f+), as in the proposition of Titsias (2009a). Similar approximations based on the stochastic process conditionals were previously used in Bui et al. (2017) and Matthews et al. (2016), with emphasis on the theoretical consistency of augmentation. Importantly, all parameters and variables needed for the approximation in (4) are given, and stored in the kth moduleMk. Then, the nested conditional expectation in (3) turns out to be
Ep(f+6=u∗ |u∗) [log p(y|f+)] ≈ K∑
k=1
Ep(uk|u∗) [ logZk qk(uk)
pk(uk)
] ,
where we applied properties of Gaussian marginals to reduce the large-dimensional expectation. For instance, the integral ∫ p(f+)df+6=uk is analogous to ∫ p(f+6=uk ,uk)df+6=uk=p(uk) via marginalisation.
Variational contrastive expectations. The introduction of K expectation terms over the log-ratios given the likelihood approximations in (4), leads to particular advantages. Having a nested integration in (3), first over u∗ given the variational density q(u∗), and second over uk for the log-ratio qk(uk)/pk(uk), we can exploit the GP predictive equation to write down
K∑ k=1 Eq(u∗) [ Ep(uk|u∗) [ logZk qk(uk) p(uk) ]] = K∑ k=1 EqC(uk) [ logZk qk(uk) pk(uk) ] , (5)
where we obtained qC(uk) via the integral qC(uk) = ∫ q(u∗)p(uk|u∗)du∗, that coincides with the approximate predictive GP posterior. The computation can be obtained analytically for each kth subset uk using the following expression
qC(uk) = N (uk|K>∗kK−1∗∗µ∗,Kkk + K>∗kK−1∗∗ (S∗ −K∗∗)K−1∗∗K∗k),
where, once again, φ∗ = {µ∗,S∗} are the global variational parameters that we aim to learn. One important detail of the sum of expectations in (5) is that it works as an average contrastive indicator that measures how well the global q(u∗) is being fitted to the modular approximations qk(uk). Without the need of revisiting any data, the GP predictive qC(uk) is playing a different role in contrast with the usual one. Typically, we assume the approximate posterior fixed and fitted, and we evaluate its performance on some test data points. In this case, it goes in the opposite way, the approximate variational distribution is unfixed, and it is instead evaluated over the inducing inputs Zk provided by each kth moduleMk. Module-driven lower bound. We are now able to simplify the initial bound in (3) by substituting the first term with the contrastive expectations presented in (5). This substitution gives us an initial version of the lower bound under the log-marginal likelihood of the GP, where no data intervene,
log p(y) ≥ K∑
k=1
logZk + K∑ k=1 EqC(uk) [log qk(uk)− log pk(uk)]− KL [q(u∗)||p(u∗)] .
However, since normalization constants {Zk}Kk=1 are difficult to compute, this immediately implies that the bound is not well defined without them. We therefore, make use of the information provided by the modules. One can show that using the alternative construction of variational lower bounds on sparse GP models, the constant satisfies logZk = L∗k + KL[q(f,uk)||p(f,uk|yk)], where L∗k is obtained from Eq. (1). Since we know the KL divergence is greater or equal to zero, and each L∗k per module is a lower bound of logZk, we have that ∑K k=1 L∗k ≤ ∑K k=1 logZk. This allows us to write down a closed-form bound LM ≤ log p(y) ,
LM = K∑
k=1
L∗k + K∑
k=1
EqC(uk) [log qk(uk)− log p(uk)]− KL [q(u∗)||p(u∗)] . (6)
The maximisation of (6) is w.r.t. the parameters φ∗, the hyperparameters ψ∗ and Z∗. To assure the positive-definitiness of variational covariance matrices {Sk}Kk=1 and S∗ on both local and global cases, we consider that they all factorize according to the Cholesky decomposition S = LL>. We can then use unconstrained optimization to find optimal values for the lower-triangular matrices L.
A priori, the module-driven bound is agnostic with respect to the likelihood model chosen. There is a general derivation in Matthews et al. (2016) of how stochastic processes and their integral operators are affected by projection functions, that is, different linking mappings of the function f(·) to the parameters θ. In such cases, the local lower bounds Lk in (1) might include expectation terms that are intractable. Since we build the framework to accept any possible data-type, we propose to solve the integrals via Gaussian-Hermite quadratures as in Hensman et al. (2015); Saul et al. (2016), and if this is not possible, an alternative would be to apply Monte-Carlo methods.

heading: 2.3 Transfer learning with multi-output modules
text: The arguments used so far are based on the idea that there exists a stationary GP across all modules, being uk ∀{1, 2, . . . ,K} values of the same function f(·). This viewpoint can be satisfied if we deal with data from the same domain — see Figure 2 for illustrative examples with MNIST images.
However, regarding the type of outputs that can be modelled, it becomes difficult when we have modules of heterogeneous data from the same phenomena, e.g. a mix of continuous, binary or discrete variables with different likelihood models. In such cases, we know that independent modules in our dictionaryM might be strongly correlated, but not parametrized by one single-output GP. As a result, we relax the stationary assumption over f to accept an additional set of independent latent functions V = {vq(·)}Qq=1. This is, we assume the output functions {fk(·)}Kk=1 per moduleMk ∈ M to be linear combinations of V . Based on Moreno-Muñoz et al. (2018), this is equivalent to introducing a multi-output Gaussian process (MOGP) prior (Alvarez et al., 2012) in the formulation.
Multi-parameter modules. Given Q latent functions V , each latent function vq(x) is assumed to be drawn from an independent GP prior, such that vq(·) ∼ GP(0, kq(·, ·)), where kq can be any valid covariance function, and the zero mean is assumed for simplicity. Since we adopt the linear model of coregionalisation (LMC) (Journel and Huijbregts, 1978), each output function fk per moduleMk is then given as fk(x) = ∑Q q=1 ∑Rq i=1 a i k,qv i q(x), where viq(x) are i.i.d. samples from each latent qth GP prior and aik,q ∈ R are their associated coefficients for i = 1, 2, . . . , Rq samples. Notice that the variational parameters per moduleMk, are given in terms of uk, the evaluation of fk(·) given the index set
Zk. Notice that the output functions {fk}Kk=1 are CI, given V . This leads to heterogeneous likelihoods that we can also factorise. The augmented distribution becomes log p(y|F+) = ∑K k=1 log p(yk|fk+), where we denote F+ = {fk+}Kk=1. Augmentation of multiple likelihoods. We follow the same formalism for latent functions V used in the single-output scenario, where we introduced the inducing inputs Z∗ = {zm}Mm=1. Sparse approximations in the context of MOGP have been previously explored in Alvarez and Lawrence (2009); Álvarez et al. (2010). Additionally, we define Q variational distributions, such that q(v∗q) = N (v∗q|µ∗q,S∗q), where v∗q = [vq(z1), vq(z2), . . . , vq(zM )]>. This gives us a way to construct approximations to multi-output likelihood terms p(yk|fk+), and particularly, we use the augmented joint distribution p(F+,V+), that factorises p(F+,V+) = p(F+,V+6=∗|V∗)p(V∗), where p(V∗) = ∏Q q=1 p(v∗q), GP priors per qth latent function are p(v∗q) = N (0,K∗q∗q) and matrices [K∗q∗q]m,n := k(zm, zn), with zm, zn ∈ Z∗. Since we can obtain closed-form conditionals between latent variables v∗q and uk from modules using the cross-covariance matrices of the MOGP, we are able to rewrite our multi-output module-driven bound as
LMOM = K∑
k=1
L∗k + K∑
k=1
EqC(uk) [log qk(uk)− log pk(uk)]− Q∑
q=1
KL [q(v∗q)||p(v∗q)] , (7)
where we also rewrote the contrastive density qC(uk) = ∫ p(uk|v∗)q(v∗)dv∗, with v∗ = {v∗q}Qq=1
and q(v∗) = ∏Q q=1 q(v∗q). This results in closed-form solution for qC(uk), which is in the appendix.
Computational cost and connections. The computational cost of training modules is O(NkM2k ), while the meta-GP reduces to O(( ∑ k Mk)M
2) and O(M2) in training and memory, respectively. If we consider the multi-output approach, the cost only increases to O(( ∑ k Mk)QM
2). Note that learning the meta-GP is no longer data-driven. The methods for distributed GPs typically need O( ∑
k N 2 k ) for global prediction. We also find a link between the module-driven bound in (6) and the
underlying idea in Tresp (2000); Deisenroth and Ng (2015). To approximate the predictive posterior, such methods combine local estimates divided by the GP prior. This is analogous to (6), where we formulate in the logarithmic plane and the variational inference setup.
(∗) Respectively, Gaussian and non-Gaussian regression (N & non-N REG), classification (CLASS), heterogeneous (HET) and storage (ST).

heading: 3 Related work
text: In terms of distributed inference for scaling up computation, that is, the delivery of tasks across parallel nodes, we are similar to Gal et al. (2014). However, their approach distributes calculus operations, i.e. product of vectors or matrix inversions that are later centralized, but no usable models under our definition of module. This is, we obtain a meta-model from models that are pre-trained and non-obsolete. Alternatively, if we look to the property of having nodes that contain usable GP models (Table 1), we are similar to Deisenroth and Ng (2015); Cao and Fleet (2014); Neumann et al. (2009) and Tresp (2000), with the difference that we introduce variational approximation methods for non-Gaussian likelihoods. An important detail is that the idea of exploiting properties of full stochastic processes (Matthews et al., 2016) for substituting likelihood terms in a general bound has been previously considered in Bui et al. (2017), ending in the derivation of expectation-propagation (EP) methods for streaming inference in GPs. There is also the method of Bui et al. (2018) for both federated and continual learning, but focused both on EP and the Bayesian approach of Nguyen et al. (2018). A short analysis of their application to GPs is included for continual learning settings but far from the large-scale scope of our paper. Moreover, the spirit of using inducing-points as pseudo-approximations of local subsets of data is shared with Bui and Turner (2014), that comments its potential application to distributed setups. More oriented to dynamical modular models, we find the work by Velychko et al. (2018), whose factorisation across tasks is similar to Ng and Deisenroth (2014) but closer to state-space methods. It is also important to mention that we are different from Duvenaud et al. (2013), because we work in the concept of models, not a dictionary of kernels as they do. We also find Mallasto and Feragen (2017), where they incorporate uncertainty in a large population analysis from the optimal transport perspective, using GPs and 2-Wasserstein metrics.

heading: 4 Experiments
text: In this section, we evaluate the performance of the module-based GP framework (MODULARGP) for multiple architectures, scenarios and data access settings. To illustrate its usability, we present results in three different learning scenarios: i) regression, ii) classification and iii) multi-output learning where we also consider heterogeneous likelihoods. All experiments are numbered from one to six in roman characters. Performance metrics are given in terms of the negative log-predictive density (NLPD), root mean-square error (RMSE) and mean-absolute error (MAE). For standard optimization, we used the Adam algorithm (Kingma and Ba, 2015). We provide Pytorch code that allows to easily learn the meta models from GP modules.1 It also includes the baseline methods used. Details about strategies for initialization and optimization are provided in the appendix. We also remark that data are never revisited and their presence in the meta-GP plots is just for clarity.
GPoE 2.79± 0.56 2.43± 0.52 1.96± 0.48 2.73± 0.72 2.19± 0.91 1.71± 0.84 2.72± 0.52 1.98± 0.61 1.32± 0.05 RBCM 2.96± 0.51 2.49± 0.51 2.02± 0.46 3.03± 0.86 2.51± 1.12 1.99± 1.04 2.56± 0.06 1.82± 0.02 1.37± 0.03 ModularGP 2.71± 0.11 1.56± 0.04 0.97± 0.05 2.89± 0.07 1.73± 0.01 1.23± 0.02 2.87± 0.09 1.87± 0.07 1.34± 0.09 Acronyms: BCM (Tresp, 2000), PoE (Ng and Deisenroth, 2014), GPoE (Cao and Fleet, 2014) and RBCM (Deisenroth and Ng, 2015).
4.1 Regression. In our first experiments for sparse variational GP regression, we provide both qualitative and quantitative results about the performance of the modular GP framework. (i) Toy concatenation: In Figure 4, we show three of five tasks united in a new GP model. Tasks are GP modules fitted independently with Nk=500 synthetic data points and Mk=15 inducing variables
1The code is publicly available in the repository: https://github.com/pmorenoz/ModularGP/.
per module. The ensemble fits a variational solution of dimension M=35. Notice that the variational meta-GP tends to match the uncertainty of the modules. (ii) Distributed GPs: We provide error metrics for the modular GP framework compared with the state-of-the-art models in Table 2. The training data is synthetic and generated as a combination of sin(·) functions (in the appendix). For the case with 10K observations, we used K=50 tasks with Nk=200 data-points and Mk=3 inducing variables in each GP module. The scenario for 100K is similar but divided into K=250 tasks with Nk=400. Our method obtains better results than the exact distributed solutions, as the meta-GP finds the average solution among all modular GPs. The baseline methods are based on a combination of solutions, if one is bad-fitted, it has a direct effect on the predictive performance. We also tested the data with the inference setup of Gal et al. (2014), obtaining an NLPD of 2.58±0.11 with 250 nodes for 100K data. It is slightly better than ours and the baseline methods. Note that it only distributes the computation of matrix products and inversions, which is equivalent to the standard model in Hensman et al. (2013).
4.2 Classification. We adapted the modular GP framework to accept non-Gaussian likelihoods, and in particular, binary classification with Bernoulli distributions. We used the sigmoid mapping to link the GP function and the probit parameters. Figure 2 shows an illustrative pixel-wise MNIST {0, 1} experiment inspired in Van der Wilk et al. (2017), where {0, 1} meta-GPs are free of data revisiting. (iv) Banana dataset: We used the popular dataset in sparse GP classification for testing our method with M=25. We obtained a test NLPD= 7.21±0.04, while the baseline variational GP test NLPD was 7.29±7.85×10−4. The improvement is understandable as the total number of inducing points used, including the ones in modules (C), is higher in the modular GP scenario.
4.3 Multi-output learning. We tested the performance of multi-output meta-GPs on large datasets. Error metrics are averaged across all output functions {fk}Kk=1 and we are interested in the difference
of predictive performance between the outputs of the meta-GP vs. modules and the standard MOGP model (Alvarez and Lawrence, 2009). Links to datasets are provided in the appendix. (v) Airline delays (US flight data): We took data of US airlines from 2008 (1.5M), where 200K samples where used for testing. The output target is the flight delay and inputs are 6 dimensional. Our goal is to analyse if having GP regression modules per day (K=366) or per month (K=12), one can obtain a meta-MOGP linking one-module per output. The results are shown in Table 3, where we see that the relative error difference (∆) per output is small. This is, the meta-MOGP prediction error only decreases around 1% approx. w.r.t. the GP modules and the variational MOGP baseline. (vi) London household: Based on Hensman et al. (2013), we obtained the register of properties sold in the Greater London County during 2017. The inputs are longitude-latitude coordinates and we use an heterogeneous model with an heteroscedastic Gaussian and the two Bernoulli likelihood distributions. We trained GP modules on the Gaussian samples and one of the binary data channels (prices and type of contracts). The scheme is provided in Figure 3. Later, we obtained a meta-MOGP from the two modules and new binary data register (type of house). The predictive curves of the metaMOGP for the binary module are provided in Figure 6. We obtained NLPD=4.18±0.06 in regression and NLPD=4.78±0.03 in the classification module. The meta-MOGP obtained NLPDf1=5.49±0.17, NLPDf2=4.64±0.06 and NLPDf3=4.51±0.11 in the Gaussian and two Bernoulli tasks respectively. The performance has a slight decrease on the regression task, while improves in the binary output w.r.t. the module. The dataset size is N=20K. We used Q=3 and M=Q×16 inducing inputs.

heading: 5 Conclusion
text: We introduced a new framework for building meta-models from independently trained GP modules. Our main contribution is to keep modules intact based on their parameters, avoid their obsolescence and mix them to form new usable tool without revisiting any data. The formulation is principled and allows for GP regression, classification and multi-output learning with heterogeneous likelihoods. We analysed its performance on synthetic and real data, and compared it with the state-of-the-art works. Experimental results show remarkable evidence that the method is robust, and successfully fits to the dictionary of modules. In future work, it would be interesting to extend the framework to include convolutional kernels (Van der Wilk et al., 2017) for large-scale image processing and functional regularisation (Titsias et al., 2020; Moreno-Muñoz et al., 2019) for continual learning applications. As a potential societal impact of our work, we argue that in the long-term, fitted models must be protected, as we do with data to the preserve privacy of individuals.

heading: Acknowledgements
text: The authors want to thank Daniel Hernández-Lobato for his constructive comments and Javier González for the useful discussion on this work during the PhD defense of PMM in Madrid last March. We also thank Søren Hauberg for the valuable feedback and the Section for Cognitive Systems at DTU for providing the computational resources. PMM has been supported by FPI grant BES2016-077626 and ERC funding under the EU’s Horizon 2020 research and innovation programme (grant agreement nº 757360). AAR acknowledges the grant TEC2017-92552-EXP (aMBITION) by Ministerio de Ciencia, Innovación y Universidades and the grants TEC2017-86921-C2-2-R (CAIMAN) and RTI2018-099655-B-I00 (CLARA) jointly with the European Comission (ERDF). He also acknowledges the grant Y2018/TCS-4705 (PRACTICO-CM) by the Comunidad de Madrid. MAA has been financed by the EPSRC Research Projects EP/T00343X/2 and EP/V029045/1.

Review summaries:
A novel approach to learn from multiple tasks (involving different partitions of predictor datasets and different types of target variables) that share the same type of input and domain. The method approximates and bounds the overall (cross-task) log-likelihood with ELBOs and posteriors of each task/module. Using inducing points posteriors, the data in each task can be summarized by inducing points of each task and is not required in later tasks. This is an appealing property in practice since it obviates the need to store all datasets of all tasks, and the training can be performed in a task-distributed fashion. 

Summary: the paper proposes a modular decomposition of the conventional GP process. The author claims that such novel decomposition can empirically improve the prediction accuracy. However, the paper lacks certain theoretical reasoning behind the claims. Many critical definitions are missing and the presentation quality of the paper is below the NeuRIPS standard. 

The authors propose a framework of reusing previously trained sparse variational GP models, without the need for revisiting the original data. Each variational sparse GP model is called a "module". The authors assume K modules, each of which is trained individually on a dataset $D_k$ (for $1 \leq k \leq K$) by means of ELBO maximisation (equation (1)). The trained modules are subsequently combined into a meta-GP, which is defined as a global variational approximation of a large-dimensional GP. This is achieved by deriving a global variational bound, where the global likelihood is approximated by utilising the local variational distribution in eq. (4). Eventually, the evidence lower bound of eq. (6) is produced, which is optimised with respect to the global variational parameters $\phi = \{\mu, S\}$, as well as the hyperparameters $\psi $ and the (global) inducing point locations $Z$. The framework is also extended to a multi-output context following MOGP. The experimental campaign demonstrates that the methodology proposed is highly competitive on a number of large-scale datasets.