Modular Gaussian Processes for Transfer Learning
Authors: Pablo Moreno-Muñoz
Antonio Artés
Mauricio A Álvarez
Venue: NeurIPS 2021 Poster (2021)
Keywords: Gaussian Processes

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, each model being characterised by its hyperparameters, pseudo-inputs and their corresponding posterior densities, one could build ensemble GP models without revisiting any data. 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 multi-task experiments, also compared with the exact inference methods in the literature.

Modular Gaussian Processes for Transfer Learning
Pablo Moreno-Muñoz" Antonio Artés-Rodriguez"
Mauricio A. Álvarez!
*Section for Cognitive Systems, Technical University of Denmark (DTU) 'Dept. of Signal Theory and Communications, Universidad Carlos III de Madrid, Spain Evidence- Behavior (cB2), Spain Dept. of Computer Science, University of Sheffield, UK pabmošdtu. dk, antonio@tsc. uc3m. es, mauricio. alvarez0sheffilld ac uk
Abstract
We present a framework for transfer learning based on modular variational Gaussian processes (GP). We develop a module-based method that having 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 nigh-dimensiona integral operators based on the Kullback-Leiblee divergence between stochastic processes to introduce 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 posteriori between independent modules. Extensive results illustrate the usability of our framework large-scale and multi- task experiments, also compared with the exact inference methods in the literature.
Introduction
Imagine supervised learning problem, for instance regression, where N data points are processed for training a model. At later time, new data are observed, this time corresponding to binary classification task, that 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 with the classification dataset to generate new model. However, this practice might be inconvenient because of i) the need of centralising data to train the model, ii) the rising date-deeendent computational cost as the number of samples increases and iii) the obsolescence of fitted models, whose future usability 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 to incorporate unseen tasks as former models recurrently discarded.
$PA
SPp
GPc
them in combination
gPp
unobserved
observed
Figure GP modules (A, B, C) are used for training (D) without revisiting data.
Alternatively, we propose 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 meta-GP using the module, but without revisiting any sample. If no data are observed, combining multiple modules is also allowed see Figure Under this framework, new family of module-driven GP models emerges.
35th Conference on Neural Information Processing Systems (NcurIPS 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- 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. non- centralised observations, asynchronous samples or missing inputs. Such settings, together with the well-known computational cost of GPs, typically ((N*), 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., [991] Hinton 2002). While the Bayesian committee machine (BCM) of Tresp (2000) focused on merging independently trained GP regression models 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>10*, 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 (0013). Their approaches 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 -eilback-Leiblerr (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.
2 Modular Gaussian processes
We consider supervised learning problems, where we have an input-output training dataset 'D- {I.,. yn NN with n,, E RD. 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)0) as p(vi.f) where the generative parameters are linked via = ((x)). We say that f(-) is a non-linear function drawn from a zero- mean GP prior f ~ GP ((. :), 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 O.
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, (Di, Dz...., Dx ). There is not any restriction on the number of modules, (Dx K {D&1 do not need to have the same size, and we only
restrict them to be N₂ <N. However, since we deal with a huge number of observations, we still consider that N₂ for all k E 21,2..... K} is sufficiently large for not accepting exact GP inference due to temporal and computational demands. Notice that & is an index while k(.,,) is the kernel.
2.1 Sparse variational approximations for independent modules
For each data module, we adopt the sparse variational GP approach based on inducing variables, and the framework Titsias (acbaas(200g)) (2009b), where posterior mean vectors and covariance matrices are computed as parameters. The use of auxiliary variables with approximate inference methods 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 M₁ N inducing inputs Z₁ {zm] Im where Zm E R²

---

Their non-ling non-linear function evaluations by f()) are denoted as u, = [f( f(2))..... f II As a first approach, we assume that f is stationary across all modules, being all un Ak {1, 21,2..... K}
function evaluations of ((.). This will be later relaxed in the following section.
To obtain the independent approximation to the posterior distribution p(f(Ds) of the GP, we introduce a kth variational distribution (f((f) for each data module Dr. In particular, this variational density factorises as qu (x( = p(J s [ux)q: (ux)) with qu (ux) - (uus Mss S&) and p(fzu (us) being the standard conditional GP prior distribution given the hyperparameters uhu of the kth kernel. To fit the GP modules, and their variational distributions qi (ux), we build lower bounds Ls on the marginal log-likelib regliiklihhod (ELBO) of every dataset Dr. Then, we use gradient- based optimisation methods to maximise the к objective functions Lt, one per module, in an separate asynchronous manner. The ELBOs, for either regression or classification tasks are obtained as follows
NL
Lx = Σ (f.) [log p(yn [f]]] - KL (ux) Ilpx (ux))
(I)
n=1
with ps(uss = =N(ux K where KAs E IR Me X has entries k( 2me) with Zm, Rm/ E Zu and conditioned to certain kernel hyperparameters the that we also aim to optimise. The variable f. corresponds to fxx) and the marginal posterior comes from qu (f,,)
p(f,, [e ) (ug )uy. In
practice, the distributed bounds Lx 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 dictionary, containing the already fitted GP modules, for their later use. Such dictionary consists, for instance, of a list of module objects [2,,.M....... Mx} without any specific order, where each M2 = { 1/ Zk, Cg),, being &x the corresponding variational parameters μ4, and Sμ.. Notice that we also include Lj. in addition to the fitted inducing points and hyperparameters. This is the ELBO obtained during optimization. will play key role later for the module-driven learning process.
2.2 Module-driven lower bound
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 (Dx KK This is only possible if we consider the entire dataset D in maximum likelihood criterion, that is, using the global log log-evidence. Specifically, our goal is now to obtain an approximate posterior distribution {f) = p(f(D), by maximising 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 longer data-driven, but module-driven instead. This is, we build 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 ctt al. (2019) does for obtaining an upper bound on the Kullback-Leibhe (KL) divergence. The idea is to use the dirgee-dmmnsional integral operators introduced by Matthews et al. (016)) in the context of variational inference, and previously by [Segerr (2002) for standard GP error bounds. The use of the large-dimen dimenssonal integrals is equivalent to augment- tinnent-and-eee 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 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 J in the dataset D and f4 the augmented nine-diimesiooall GP. Notice that f+ sufficiently large and it contains the function values taken by f(-) at both {I,))_ and { Zz K
for all data modules. We consider it finite. The log- marginal expression is therefore
log p(y) = log p(y), ......, yK) = log /
No are the output data-points already used for training each GP module Mr. The joint distribution in factorises according to p(y).+))(({+) where the l. 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.) = Nu, (ua. S.)
where yı: = [
3

---

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 m=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. (204)) 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 is inequality, it is obtained as
logr(y) = log /f
wu p(y)s+ ((ua)) plun) +ff+ du.
B.p Bepsstee [ogg(y]]+)] + log p(ux)
2 Ryy(
(3)
q( ((uu))
where we applied the properties of Gaussian conditionals to factorise the GP prior as p(f+)
p(f+/u. ur ).(.). Notice that f+ = {ftzu. ). Here, the last prior distribution is p(u,) =
Nu,, K..) where [m,, == k(zm; Zn), with 2m, 2n E Z.. We also consider k(;,)) condi- tioned the global kernel hyperparameters s/. that we also aim to estimate. The double expectation in B comes from the factorization of the carge-dimena Jensen's 's inequality twice. Its derivation is in the appendix.
integral given f+ and the application of
Likelihood reconstruction. The augmented likelihood distribution p(ylf+) is perhaps the most important part of the derivation. It allows the data modules DR. and particularly, the output variables yk- This gives a factorized term that we will later use for introducing the GP modules MM stored in our dictionary, that is, log f4)] = =k log f+)) To avoid revisiting each kth likelihood term, and hence, eval- uating 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 f4. Bayes theorem indicates that each kth likelihood can be approximated
apply conditional independence (CI) between
p(uxlf+) = Z' (e({+) = 2% p(f+ U1x (u))) px(f+)
plff Tee Ok (ux)
9k (UAA)
Zk
(4)
ps py(u)))
since qu (f+) s p(f+ (v). Here, Za refers to the normalization constant or evidence of Dz. which still difficult to compute. We also assumed that the augmented approximate density factorises as
(f+ = pffff((uu 14 (2 using the structure of the GP prior (f((+). ). as in the proposition of Titsias (009a)) Similar approximations based on the stochastic process conditionals were previously used in Bui et al. (017) and Matthews et al. (0066), 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 module Ma. Then, the nested conditional expectation in 3 turns out to
be
K
9k (Uk)
(pp(st (u.) [logp(y][+) R
Σ Ep(uz/us)
log Zk
Pl: Pa.(ua))
where we applied properties of Gaussian marginals to reduce the large- dinge-dimensional expectation. For instance, the integral J p(f+)df )df+ is analogous to f'f(ff Mff+/uu Un 4x y =p(ux) via marginal- isation.
Variational contrastive expectations. The introduction of K expectation terms over the log-ratios given the likelihood approximations in 4 leads to particular advantages. Having nested integration in B first over given the variational density ((u)), and second over Uy for the log-ratio qe (up)/pp (ux), we can exploit the GP predictive equation to write down
Σ Ey(uu)) (ppcuslu))
K
log Zn
Qk (uR)
]
Σ
K
Eye (nee) (vee)
[log Ze Ts (Ux))
(5)
p(u)))
pk (21z.
where we obtained ge (us) via the integral (c(ux) )==)= = { approximate predictive GP posterior. The computation can be obtained analytically for each kth subset un using the following expression
g(u. p(ux |ua )du., that coincides with the
qc (ux) Nuu
Ms KEN KIrK,, (S. K.. J.. K.t);

---

(A) Modular
(B) MNIST GP
Modialar
(D) MNIST (1) GP
1.0
66
10
60
-4.5
-10-25 nnD a5
L0
-83
5n 65 14
-1.0 -033
0.1
0.5
1.0
-10 -05
DV0
14
-10
-1.5
0.33
0.5
Eaput
input
icpus
21 input
iapus
Figure 2: Modular GPs for {0, 1} MNIST data samples. The meta- GPs (B-D) are built from fitted modules and do not revisit samples, only parameters and variables stored in each MM.
where, once again, Φ. (ue. ) are the global variational parameters that we aim to learn. One important detail of the sum of expectations (5) is that it works an average contrastive indicator that measures how well the global ((u.) is being fitted to the modular approximations ga (ux) Without the need of revisiting any data, the GP predictive (c(ss) is playing a different role 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 the opposite way, the approximate variational distribution is unfixed, and it is instead evaluated over the inducing inputs Z provided by each kth module Mr.
Module-driven lower bound. We are now able to simplify the initial bound in by substituting the first term with the contrastive expectations presented (5). This substitution gives us an initial version of the lower bound under the log- marriial likelihood of the GP, where no data intervene,
K
K
poy) N log Z +
Σ k=1
(ny) [log qu (ux) - lo (ux)) ) - KL [g(u.))lp(u..
However, since normalization constants { are difficult to compute, this immediately implies that the bound by the modules. One can show that using the alternative construction of variational lower bounds on sparse GP models, the constant satisfies log Zx - C + KL[q(f, u (ps, us (yk)), where Lz is obtained from Eq. I Since we know the KL divergence is greater or equal to zero, and each 2 per module is lower bound of log ZR. have that Ek= 2% < > Lam log Zμ. This allows us to write down a closed-form bound LM S logp(y).
well defined without them. We therefore, make use of the information provided
K
K
LM = Σ
C
+
E Enc( [log 9k (ux) (ogg(();) - KL ig (q(uu))lp((x))
(6)
Jmm
The maximisation of (6) is w.r.t. the parameters dh, the hyperparameters v/. and Z. To assure the positive-definitines. of variational covariance matrices [S 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 find optimal values for the lower-triangular matrices L.
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 0. In such cases, the local lower bounds C. in (I) 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 in Hensman al. (2015); (Suul et al. (2016), and if this is not possible, an alternative would be to apply Monte-Carte CCalo methods.
2.3 Transfer learning with multi-output modules
The arguments used so far are based on the idea that there exists stationary GP across all modules, being us 1 ... deal with data from the same domain see Figure 2 for illustrative examples with MNIST images.
2,.... K} values of the same function f(-). This viewpoint can be satisfied if we
5

---

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. mix of continuous, binary or discrete variables with different likelihood models. In such cases, we know that independent modules in our dictionary M might be strongly correlated, but not parametrized by one single-output GP. As result, we relax the stationary assumption over f to accept an additional set of independent latent functions
V= (ug(-)) This is, we assume the output functions (fx( ( per module M₁ E M to be linear combinations of V. Based on MorrnooMMuño et al. (2018), this is equivalent to introducing a multi-output Gaussian process (MOGP) prior (Alvarez et al. 20122 in the formulation.
Multi-parameter modules. Given Q latent functions V, each latent function
Regression GP
Classiication GP
New data
Py (2:) assumed to be drawn from an independent GP prior, such that (.).)
CP(O, ())),, where ki 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 1& per module Mg is then given as fx(()) = - E ak.g (xr), where oy are i..dd samples from each latent gth GP prior and ai.q E R are their associ- ated coefficients for = =, ........ Rg sam- ples. Notice that the variational parameters per module M&. are given in terms of utg. the evaluation of f.). given the index set Zz Notice that the output functions {Lk } are CI, given V. This leads to heterogeneous likelihoods that we can also factorise. The augmented distribution becomes log = ΣK log pog where we denote F+ (ffx- 1
M1
M2
2. 3/
fO
fz(()
fz())
Multi-output GP
Ensemble
Figure 3: Illustration of modular MOGP. Regression,
classification and data tasks are linked to outputs fx)()-
We follow the same formalism for latent functions V used in the single-output scenario, where introduced the inducing inputs Z. = (zm 1M 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 (vg)) ( plag, S.g), where Uag (v.y(21), (22 -e)))).---yy(// [2([ This gives way to construct approximations to multi-output likelihood terms p(ya: fu+), ). and particularly, we use the augmented joint distribution p(F+, V+), that factorises p(F+, V+) = pFF+, V+ Vv)p(V.), where
Augmentation of multiple likelihoods.
p(V.) = [[= p(vag), GP priors per qth latent function are p( (((gg) = WO, K and matrices IK. = k(zm, 2n), with 2m, Zn E Z.. Since we can obtain closed-form conditionals between latent variables Uny and from modules using the cross-covariance matrices of the MOGP, we are able to rewrite our multi-output module-driven bound as
K
K
MM
Σ
E (ny) (wx) [log 9u (un) - log Pa (ux))
Σ KL (g(
(7)
1=1
k=1
where we also rewrote the contrastive density (c(ux) ((uull=))(vv)dde with v, = fw g== and g(v.) = _ q(v.g). This results in closed-form solution for qe ge(ux) (ux), which is in the appendix.
Computational cost and connections. The computational cost of training modules is M()) while the meta-GP reduces to O M2) and ((M2) in training and memory, respectively. we consider the multi-outie multi-output approach, the cost only increases to O( MM)2) Note that learning the meta-GP is no longer data-driven. The methods for distributed GPs typically need OEE, N2) for global prediction. We also find a link between the module- driven bound in (6) and the underlying idea in Tresp Deisenroo 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.

---

Table I: Properties of distributed/modular GP models
MODEL
N REG.
non-N REG.
CLASS.
HRT.
INFERENCE
GNOOEE
FREE OF DATA ST.
Tresp (2000
x
X
K
Analytical
X
Ng and Deisenroli 2014
X
X
X
Analytical
X
Cao and Flact 2014
X
X
X
Analytical
X
Deisenroth and Ng 2015
X
X
X
Analytical
X
Gal at 2014
Variational
work
Variational
(*) Respectively. Gaussian and non-Gaussian regression (N non-N RBG), classification (CLASS), beterogencous (HET) and storage (ST).
3 Related work
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 meta-model from models that are pre-trained and non-obsolete. Alternatively, we look to the property of having nodes that contain usable GP models (Table M.. 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-Gane nau-saus likelihoods. An important detail is that the idea of exploiting properties of full stochastic processes (Matthews et al.] (016)) for substituting likelihood terms in general bound has been previously considered in Bui al. (2017), ending in the derivation of expectation expeecationn-poopaattie 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 iet 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-approximaatins 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.] (0188), whose factorisation across tasks is similar to Ng and Deisenroth (2014) but closer state-space methods. It 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 Mallasto and Feragen (2017), where they incorporate uncertainty in a large population analysis from the optimal transport perspective, using GPs and 2- 2-Wasserstein metrics.
4 Experiments
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 six in roman characters. Performance metrics are given in terms of the negative log-predictive density (NLPD), root mean-square error (RMSE) and mean-abss 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 It also includes the baseline methods used. Details about strategies for initialization and optimization are never revisited and their presence in the meta-GP GP plots is just for clarity.
provided in the appendix. We also remark that data
Table 2: Comparative error metrics for distributed GP models.
DAIA siZE -3
10x
100k.
IM
MODEL
NLPD
RMSE
MAB
NLPD
RMSE
MAE
NLPD
RMSE
MAE
BCM
2.90 0.94
11.94 18.80
2.05 1:
3.51 0.775
2.35 0.96
1.34 1 1.05
NA
NA
PoB
2.79 0.15
2.32 2.22
1.8620222
2.82 0.67
2.19 0.91
1.71
291
0.63
1.98
1.32 2: 0.85
GPoF
2.70+00.6
1.96 0.48
2.73 0.772
2.19-00.1
+-..4
2.72
0.52
1.98 0.61
1.32 0.05
RBCM
2.06 0.5
2.49 0.51
2046
3.08
2.51
2.56
0.0G
1.37 0.03
ModolarP
2.71 0.11
1 1.56
0.04
0.97
0.055
2.80:1..00 0.07
1.78
0.01
1.23 0.02
2:87 1 0.00
1.87
0.07
1.34 0.00
Acronyms: BCM (Tresp 2000
PoE
(Ng and Deiseuroth 2014),
GPoE (Cao and Flcet 20144 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 we show three of five tasks united in new GP model. Tasks are GP modules fitted independently with 55= synthetic data points and M;=15 inducing variables
'The code is publicly available in the repository: https: httptttwww..om//omoeee////ttuuta
7

---

Table 3: Comparative metrics of modular multi multi-output GPs for US-FLIGHT dataset.
PARTITION
MODULES (f)
2.36
0.18
1.48 d: 0.26
2.31
5: 0.24
2.03
3: 0.02
1.53 0.06
1.83
0.03
MODULARGP
1.58 1.14
MODELARGP (Q
2.38
+
2:38
0.113
2.38
MODULARGP (Q 4)
2.36 0.15
1.49 2 0.226
2.31
1 0.24
2.39 + 0.03
1:014
2.37 +
0.12
MOGP
2.49
1.51 1:
2.31
1:
1.61 1.112
2.23 1
-3)
2.39
1.11
2.18
MOGP (Q
2.37
0.17
1.51
+ 0.26
2:31
4 0.25
2.34 =: 0.28
1.63 0.11
2.14
0.13
PARTITION
DAYS
MONTHS
Ave. DIEFERENCE PER OUTPUT -
ANLPD
AMAE
ARMSE
ÁNLPO
AMAE
ARMSE
MODULARGP QQ 2] Modules
-3.91%
-0.64%
-0.17%
-17.69%
-1.41%
-22.73%
MODULABGP (Q 3) Modules
-0.51%
-0.99%
-0.16%
14..19%
1.41%
22.81%
MODULARGP (Q
Modules
+0.13%
-0.19%
-14.93%
-1.77%
-22.34%
MonmLARGPP (Q
-3.71%
-0.75%
-0.15%
-4.21%
-18.04%
MODULARGP (Q = 3) MOGP
-0.91%
-1.25%
-0.33%
-13.96%
-3.32%%
-15.669
MODULARGG (Q) 4) MOGP
-0.11%
-1.50%
-0.27%
-11.54%
-4.92%
-13.84%
(f) Modales metric reference.
module. The ensemble fits variational solution of dimension M-35. Notice that the variational meta-GP tends match the uncertainty of the modules.
(ii) Distributed GPs: We provide error metrics for the modular GP framework compared with the statc-of-ther thr---rrthe-rtttte models in Table The training data synthetic and generated as a combination of sin(.) functions (in the appendix). For the case with 10K observations. we used K==00 tasks with Nx=200 data-points and =3=3 inducing variables in each GP module. The scenario for 100k is similar but divided into K tasks with 40- 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 the predictive performance. We also tested the data with the inference setup Ga et al. (2014) obtaining NLPD of 2.58++0.11 with 250 nodes for 100K data. It slightly better than and the baseline methods. Note that it only distributes the computation of matrix products and inversions, which is equivalent the standard model in Hensman et 2013)
combination of solutions, if is bad-fitted, it has direct effect
(iii) Modular Meta-GPs: For a large synthetic dataset (N=10%), tested the modular GP framework with K K== 10 tasks as shown in Table However, if join large amounts of modular GPs, e.g. K 10 it is often problem- atic for baseline methods, due partitions must be revisited for building final predictions. So data revisiting is totally avoided. This is not constraint for us, as we work with the ex- pectation of modules. Table 2 (1 shows that we are closer to the exact inference solutions, even large-scale cases. Additionally, we repeated the experiment in pyramidal way see appendix. That is, building meta-GPs from modular meta-GPs, inspired in Deisen- roth and Ng (2015). Our method obtained [NIPD=4.15, MMFF-2.71. MAE-2.22
Congut,
Stere
Congett
Tupul,
Figure 4: Modular GPs with synthetic data. Three of five tasks are plotted in the upper row.
4.2 Classification. We adapted the modular GP framework to accept non- Gausssn likelihoods, and in particular, binary classification with Bernoulli distributions. We used the sigmoid mapping 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
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 78. x05x100+ 14 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 {ff 1 and we are interested in the difference

---

Banana GP Eusemble
London household data
Lattuude (ra)
Langitude (ea)
Figure 5: Modular GP classifier for banana dataset.
Figure 6: Modular MOGP f2(-)
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= one obtain meta-MOGP linking me-module per output. The results
shown in Table Bl where
see that the relative error difference (A) 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 pargitude-lattude coordinates and we use an heterogeneous model with an heteroscedastic Gaussian and the two Bernoulli likelihood distributions. We trained GP modules 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 meta- meta-MOGP from the two modules and new binary data register (type of house). The predictive curves of the meta- MOGP for the binary module are provided in Figure 6 We obtained NLPD=4.1811 10t....88000666 in regression and NLPD=4. 7..-
the classification module. The meta-MOGP obtained NLPD,, 0.17 NLPD 4.PD,,,464400.0000 0.06 and NLPD,, f, =4.51=00 =0.11 11 in the Gaussian and two Bernoulli tasks respectively. The performance has slight decrease on the regression task, while improves in the binary output w.r.t. the module. The dataset size is 2-20 We used Q-3 and M- 16 inducing inputs.
5 Conclusion
We introduced new framework for building meta-modek, 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 otate-of-the-arr works. Experimental results show remarkable evidence that the method 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 af.) (017) for large-scale image processing and functional regularisation (Titsias et al. | Moreno- Muñoz 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, do with data to the preserve privacy of individuals.
Acknowledgements
The authors want thank Daniel Hernández- 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 BES- 2016-077626 and ERC funding under the EU's Horizon 2020 research and innovation programme (grant agreement n 757360). AAR acknowledges the grant SESCC007--2552-BBB by Ministerio de Ciencia, Innovación y Universidades and the grants TEC2017-86921-C22
(aMBITION)
(CAIMAN) and RTT201 8-099655--B ----0 (CLARA) jointly with the European Comission (ERDF). He also acknowledges the grant Y2018/TCS-477 (PRACTICO-CM) by the Comunidad de Madrid. MAA has been financed by the EPSRC Research Projects EP/T00343X/2 and EP/V029045/L.

---

References
M. Alvarez and N. D. Lawrence. Sparse convolved Gaussian processes for multi-output regression. In NIPS 21, pages 57-64, 2009.
M. Álvarez, D. Luengo, M. Titsias, and N. Lawrence. Efficient multioutput Gaussian processes through variational inducing kernels. In AISTATS, pages 25- 32, 2010.
M. A. Alvarez. L. Rosasco, N. D. Lawrence, et al. Kernels for vector-vaie valued functions: A review. Foundations and Trends@ in Machine Learning, (33):1952666) 2012.
T. D. Bui and R. E. Turner. Tree-structured Gaussian process approximations. In Advances in Neural Information Processing Systems (NIPS), pages 2213-2221, 2014.
T. D. Bui, C. Nguyen, and R. E. Turner. Streaming sparse Gaussian process approximations. In Advances in Neural Information Processing Systems (NIPS), pages 5299-3307 2017.
T. D. Bui, C. V. Nguyen, S. Swaroop, and R. E. Turner. Partitioned variational inference: A unified framework encompassing federated and continual learning. arXiv preprint arXiv:1811 11200 2018.
D. R. Burt, C. E. Rasmussen, and M. Van der Wilk. Rates of convergence for sparse variational Gaussian process regression. International Conference Machine Learning (ICML), pages 862-871, 2019.
Y. Cao and D. J. Fleet. Generalized product of experts for automatic and principled fusion of Gaussian process predictions. arXiv preprint urXiv:1410.782: 2014.
M. Deisenroth and J. Ng. Distributed Gaussian processes. International Conference on Machine Learning (ICML), pages 1481-1490, 2015.
D. Duvenaud, Lloyd, Grosse, Tenenbaum, and G. Zoubin. Structure discovery in nonparametric regression through compositional kernel search. In International Conference on Machine Learning (ICMI.), pages 1166-1113 174, PMLR, 2013.
Y. Gal, M. Van Der Wilk, and C. E. Rasmussen. Distributed variational inference in sparse Gaussian process regression and latent variable models. In Advances in Neural Information Processing Systems (NIS), pages 3257-3265 2014.
I. J. Goodfellow, M. Mirza, D. Xiao, A. Courville, and Y. Bengio. An empirical investigation of catastrophic forgetting in gradient-based neural networks. In International Conference Learning Representations (ICLR), 2014.
J. Hensman, N. Fusi, and N. D. Lawrence. Gaussian processes for big data. In Uncertainty in Artificial Intelligence (UAI), pages 282-290,, 2013.
Hensman, Matthews, and Z. Ghahramani. Scalable variational Gaussian process classification. Artificial Intelligence and Statistics (AISTATS), pages 351-360, 2015.
Hinton. Training products experts by minimizing contrastive divergence. Neural Computation, 14(8): 1771-1800, 2002.
M. D. Hoffman, D. M. Blei, Wang, and J. Paisley. Stochastic variational inference. Journal of Machine Learning Research (JMLR), 14(1):1303-1347, 2013.
R. A. Jacobs, M. I. Jordan, S. J. Nowlan, and G. E. Hinton. Adaptive mixtures of local experts. Neural Computation, 3(1):79-87, 1991.
A. G. Journel and C. J. Huijbregts. Mining Geostatistics. Academic Press, London, 1978.
D. P. Kingma and Ba. Adam: A method for stochastic optimization. 2015.
Lawrence. Probabilistic non-linear principal component analysis with Gaussian process latent variable models. Journal of Machine Learning Research, 6(11), 2005.
A. Mallasto and A. Feragen. Learning from uncertain curves: The 2- Wasssersteinn metric for Gaussian processes. In Advances in Neural Information Processing Systems (NIPS), pages 5665-5674, 2017.
A. G. d. G. Matthews, J. Hensman, R. Turner, and Z. Ghahramani. On sparse variational methods and the Kullback-Leii Leibler divergence between stochastic processes. In Artificial Intelligence and Statistics (AISTATS), pages 231-239, 2016.
10

---

P. Moreno-Muñooz, A. Ariés-Rodriguez, and M. A. Álvarez. Heterogeneous multi-output Gaussian process prediction. In Advances in Neural Information Processing Systems (NeurIPS), pages 671 1711-0720, 2018.
Moreno-Mee Moreno-Muñoz, A. Artés-Rodrígee and M. Álvarez, Continual multi-tast task Gaussian processes. arXiv preprint 0rXiv:1911.00000 2019.
M. Neumann, K. Kersting, Z. Xu, and Schulz. Stacked Gaussian process leamning. IEEE International Conference on Data Mining, pages 387-396. IEEE, 2009.
Ng M. P. Deisenroth. Hierarchical mixture-of expertss model large-scale Gaussian process regression. arXiv preprint (rXi::/412.3078, 2014.
V. Nguyen, Li, T. Bui, and E. Tumer. Variational continual learning. In International Conference on Learning Representations (ICLR), 2018.
D. Peterson, P. Kanani, and V. J. Marathe. Private federated learning with domain adaptation. Workshop on Federated Learning for Data Privacy and Confidentiality at NeurIPS, 2019.
C. Rasmussen information processing systems (NIPS), pages 881-888 1888, 2002.
Ghahramani. Infinite mixtures Gaussian process experts. Advances neural
C. E. Rasmussen and C. Williams. Gaussian processes for machine learning, volume MIT Press, 2006.
F. J. Ruiz, M. K. Titsias, A. B. Dieng, and D. M. Blei. Augment and reduce: Stochastic inference for large categorical distributions. International Conference on Machine Learning (ICML), 2018.
A. D. Saul, J. Hensman, A. Vehtari, and N. D. Lawrence. Chained Gaussian processes. In Artificial Intelligence and Statistics (AISTATS), pages 1431 1440, 2016.
M. Seeger. PAC-Bayesian generalisation emror bounds for Gaussian process classification. Journal Machine Learning Research (JMLR). 3(0e)):233-269, 2002.
V. Smith, C.-K. Chiang. M. Sanjabi, and A. S. Talwalkar. Federated multi-task learning. In Advances in Neural Information Processing Systems (NIPS), pages 4424 4434, 2017.
M. Titsias. Variational learning of inducing variables in sparse Gaussian processes. In Artificial Intelligence and Statistics (AISTATS), pages 567-574, 2009a.
M. K. Titsias. Variational model selection for sparse Gaussian process regression. Technical Report, University of Manchester, 2009b.
M. K. Titsias, J. Schwarz, A. G. de G. Matthews, R. Pascanu, and Y. W. Teh. Functional regularisation for continual learning with Gaussian processes. International Conference on Learning Representations (ICLR), 2020.
Tresp. A Bayesian committee machine. Neural Computation, 12(11):719-274: 2000.
M. Van der Wilk, C. E. Rasmussen, and Hensman. Convolutional Gaussian processes. In Advances in Neural Information Processing Systems (NIPS), pages 2849-2858, 2017.
D. Velychko, B. Knopp, and Endres. Making the coupled Gaussian process dynamical model modular and scalable with variational approximations. Entrupy, 20(10)722 2018.


Reviews:
[{"review_id": "lcZ7fw123dN", "reviewer": "Reviewer_ntYs", "summary": "The paper proposes an ensemble GP method (where each model is termed a 'module') that can fuse a set of pre-existing GPs to train a new \"meta-GP\".  The authors also present a closed form ensemble lower bound that considers the parametrisations of the 'modules'.", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I rather enjoyed the motivating dialogue and examples. \n\n\"Each module is considered to be a sparse variational GP, containing only variational parameters and hyperparameters.\" Is it a necessary condition that each module be specifically a sparse variational GP? Could this approach be extended/directly applied to spectral approximations like in [R1]? I ask because given the generality of the proposal \"We present a framework based on modular Gaussian processes\" it is not so clear that the method applies only to a particular _class_ of GP derivatives (of which sparse vi methods are just one), or a more general class of GP approximations (of which, e.g., [R1] is also one).\n\nClarity:\nThe introduction starts with a rather involved sentence: \"We develop a module-based method...revisiting any data\". The sentence is rather convoluted and I would suggest condensing or splitting it up.\n\n\"n. However, GP models are not immune to settings where we need to adapt to irregular ways of observing data, e.g. asynchronous samples or missing inputs.\" I'm not so sure this is true (or at least the way it is phrased here). GPs are typically derived on continuous domains; i.e. ones that do not need data to be 'gridded' (although approximate methods have been proposed to accommodate structured data such as KISS-GP).\n\nRelated work:\nI think [R0] should be included in the related work since it is addressing a similar problem but from the optimal transport perspective and applies in a potentially more general sense than KL divergence between stochastic processes.\nAlso, the approach of a \"meta\" GP seems directly related to the classic concept of model stacking. This has been attempted in the past (e.g. [R2]) and has also been termed a meta-gp. Could the authors comment on the relation to such prior work?\n\nCompression: It seems to me that using a sparse VI gp, there is a consequent loss of information. If it was a full GP, there would be no information loss in the sense of the full covariance (of the 'past') data. \n\n[R0]  Mallasto, Anton, and Aasa Feragen. \"Learning from uncertain curves: The 2-Wasserstein metric for Gaussian processes.\" Advances in Neural Information Processing Systems. 2017.\n[R1] Lázaro-Gredilla, Miguel, et al. \"Sparse spectrum Gaussian process regression.\" The Journal of Machine Learning Research 11 (2010)\n[R2] Neumann, Marion, et al. \"Stacked Gaussian process learning.\" 2009 Ninth IEEE International Conference on Data Mining. IEEE, 2009.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "NA", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "4"}, "extra_text": {}}, {"review_id": "f2VVxs9ONEJ", "reviewer": "Reviewer_abDr", "summary": "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. ", "questions": "", "limitations": "", "rating": 9, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "To be able to learn first on a large dataset and then quickly transfer to other tasks of smaller datasets is critical for modern machine learning tasks. Thus, I think the problem addressed by the authors is of high importance. The method proposed is novel and interesting (with small caveats not fully exposed), combining inspirations from previous work to address the transfer learning problem. In particular, the transfer does not require storing the full datasets of previously encountered tasks, which is appealing to many practical problems. Further, this method can be extended to multi-output and multi-dimensional GP latent, and can be stacked on top of each other according to tasks splitting to obtain a hierarchical meta-GP model. \n\nThe experiments are thorough, covering small and large datasets. although I am very curious about the results in other well-controlled settings. Despite this small imperfection and other minor issues I list below, I really enjoyed reading this paper. The paper is also written fairly clearly with nice figures explaining the intuitions, but some descriptions are not intuitive to readers like me, so please provide more details. \n\nI am not completely up-to-date with the literature in this area, but I have not seen any methods that are trivially similar to this proposed here. The benchmark methods compared in the experiments are also quite different in methodology (though distributed in different senses), some are significantly less well-developed than this work. This argues for the novelty of this work, and if it is true that this method is that novel, the authors should further highlight this new contribution. \n\nCould the authors provide details to the questions below? If addressed well, I am willing to further raise my score. \n1. How does the approximation in (4) impact the bound? This is probably the only loose end in this method, and I hope to see some analysis in this respect. Are there cases when this approximation is particularly bad so that the objective is no longer a bound on the likelihood? Are there any tricks used to ensure the quality of this approximation?\n2. Although we do not need to retain all the data for each task, we do retain the inducing points to compute $q_k$ for (6) and (7). How does the relative size of the inducing set to the task data set affect the performance of this method?\n3. In the toy concatenation experiment, I think it would be great to analyse all the moving parts of this method. What are the effects of \n   1) inducing points and its approximate posterior,\n   2) approximation made in (4) \n   3) modular training of hyperparameters? Importantly, I'd love to see a comparison to the best solution possible without modular training, for example, a full posterior computed exactly without inducing points. \nIf it's not quite possible to clearly dissect the contribution of these factors, please explain why, as it is not immediately clear to me. \n4. In the description of the banana dataset experiment, why is it not possible to control the number of inducing points? How would the two methods compare if this number is controlled?\n4. All the figures need proper captions. I know the space is limited, but the figures make less sense and the readers have to guess. \n5. The bold and normal Z's refer to unrelated objects, please consider using other symbols?\n6. In line 229, what does $i$ index? Is this related to LMC?\n7. In line 313, what is the exact inference solution? Exact in what sense? \n8. In line 303, what is the \"combination of solutions\" referring to? \n9. Figure 5 is not referred to. \n10. There should at least be some discussions in the main text on training and optimization. \n\n===== post rebuttle ======\n\nI thank the authors for their reply. I also do not see serious technical concerns from other reviewers, except for the presentation and minor details (I now recall some initial difficulty in reading which is expressed in other reviews). This is a great paper and I raise my score and confidence. The value of this paper may be more obvious in the near future. \n\nI would strongly encourage the authors to revise the paper to include the content of ALL their responses provided here.", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "The authors mentioned data retention its impact on privacy. The authors could discuss how inducing points encapsulate or hides the information about each individual data, and if the dataset in each task reflects properties of a user group, does the inducing point reveal specific information about this group?", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "3"}, "extra_text": {}}, {"review_id": "WV1la5pBdVy", "reviewer": "Reviewer_2qsL", "summary": "The authors propose a framework of reusing previously trained sparse variational GP models, without the need for revisiting the original data.\nEach variational sparse GP model is called a \"module\".\nThe 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)).\nThe trained modules are subsequently combined into a meta-GP, which is defined as a global variational approximation of a large-dimensional GP.\nThis is achieved by deriving a global variational bound, where the global likelihood is approximated by utilising the local variational distribution in eq. (4).\nEventually, 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$.\nThe framework is also extended to a multi-output context following MOGP.\nThe experimental campaign demonstrates that the methodology proposed is highly competitive on a number of large-scale datasets. ", "questions": "", "limitations": "", "rating": 7, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The variational framework that is proposed to combine pre-trained GP modules into a single approximating model is well-motivated and novel, to the best of my knowledge.\nThe paper is well-written and the related work is adequately discussed.\nThe methodology proposed is demonstrated on a number of examples, while the numerical experiments confirm that the method is competitive to the state-of-the-art on a number of datasets, some of which involve millions of instances.\nI think this work is of interest to the community, as it may have significant applications not only due to its scalability properties, but also due to the possibility of reusing trained models based on a rigorous framework.\n\n\nThe are however a few points that need to be clarified.\n\nI found the discussion of the VEM algorithm (Algorithm 1) in the appendix somewhat confusing.\nNo EM algorithm has been discussed or motivated in the main paper.\nBefore reading the appendix, I had been under the impression that the variational parameters ($\\mu$ and $L$), the hyperparameters $\\psi$ and the inducing locations $Z$ and were being jointly optimised with respect to the ELBO in Equation (6).\nIn the joint optimisation case, $\\psi$ and $Z$ are treated in terms of empirical Bayes (due to the ELBO).\nHowever, according to the VEM algorithm, these three classes of parameters are optimised independently in three nested optimisation loops.\nIn particular, the optimisation of $\\mu$ and $L$ is called as the E-step, followed by 2 M-steps for $\\psi$ and $Z$, which are treated as latent variables.\nI do not understand the motivation behind this choice.\n\nIn line 176, it is stated that the sum of expectations in (5) works as an average contrastive indicator.\nI am not sure that I understand that \"contrastive\" means in this context. \nCould you elaborate on this?", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "In terms of societal impact, the authors comment on the need for protection of fitted models, as these may contain sensitive information.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "8"}, "extra_text": {}}, {"review_id": "QAPidSxrMv5", "reviewer": "Reviewer_FB5F", "summary": "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. ", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Weakness:\n\nThe presentation made the paper hard to follow, without a clear story line to seize. For example, why is it necessary to introduce the inducing variable $z_k$, and how are $z_k$ computed from D_k given x and y? Things get even worse when the subscript notation * is introduced, making the notations harder to interpret. Many notations are not defined，such as $L_k^*$, $f_{+\\neq u_*}, V_{+ \\neq *}$. Most importantly, the paper lacks discussions that why decompostion on the data is advantageous for GP learning. I listed some of my confusions below. \n\n79-81 “consider that N_k for all k ∈ {1, 2, . . . , K} is sufficiently large for not accepting exact GP inference due to temporal and computational demands.” What does this sentence mean? \n\n89-90 ”…being $u_k$ values of the same function $f(\\cdot)$” (grammarly wrong).\n\n144: notation $f_{+\\neq u_*}$ is not defined here.   \n\n153: why $p(y|f_+)=\\sum_k \\log p(y_k|f_+)$? Are the marginal of $ p(y) $and $p(y_k)$ expected to be same (I think so, because k seems to be indexing a subsets of the entire data., i.e., samples of size $N_k$ of the entire data)? If so, then there should be an additional multiplicative 1/K constant on the r.h.s term.\n\n178-179: “Without the need of revisiting any data, the GP predictive $q_C(u_k)$ is playing a different role in contrast with the usual one.” This sentence is seemingly trying to support the novelty claim, but failed to deliver the message. To be honest, the whole discussion paragraph that claims the advantage of such modular method (175-182) is difficult to interpret even if one with proper GP background is 100% focused. I spent hours on this paragraph but still cannot get the insights behind. \n\n189: $L_k^*$ is not defined. \n\n240: the reasoning behind such factorization is unclear, partly because the notation V_{+ \\neq *} is undefined. \n\nIt leaves one wonder how the GP helps the transfer learning task in the empirical section at all. Please clarify this point. It seems to me the tasks are conventional Regression/Classification GP tasks, without the need of transfer learning. Please clarify how the advantage of the proposed method is reflected in the experiments section.\n\nThe usage of word \"accept\" has been abused everywhere. ", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "Yes, the authors adequately addressed the limitations and potential negative societal impact of their work.", "review_text": "", "extra_scores": {"needs_ethics_review": "No", "ethics_review_area": "", "time_spent_reviewing": "6 hours"}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper proposes a framework for reusing previously learned sparse variational GP models without revisiting any data and a further extension to the multi-output context. All four knowledgeable reviewers recommend acceptance and I agree with them. As a clarification to the authors, yes, you are allowed and encouraged to incorporate amendments/additions to address the reviewers feedback in the final version of the paper. 

Decision: Accept (Poster)