Sparse Gaussian Process Variational Autoencoders
Authors: Matthew Ashman
Jonathan So
William Tebbutt
Vincent Fortuin
Michael Arthur Leopold Pearce
Richard E Turner
Venue:  (2021)
Keywords: Gaussian process
variational inference
variational autoencoders
Bayesian inference

Abstract: Large, multi-dimensional spatio-temporal datasets are omnipresent in modern science and engineering. An effective framework for handling such data are Gaussian process deep generative models (GP-DGMs), which employ GP priors over the latent variables of DGMs. Existing approaches for performing inference in GP-DGMs do not support sparse GP approximations based on inducing points, which are essential for the computational efficiency of GPs, nor do they handle missing data -- a natural occurrence in many spatio-temporal datasets -- in a principled manner. We address these shortcomings with the development of the sparse Gaussian process variational autoencoder (SGP-VAE), characterised by the use of partial inference networks for parameterising sparse GP approximations. Leveraging the benefits of amortised variational inference, the SGP-VAE enables inference in multi-output sparse GPs on previously unobserved data with no additional training. The SGP-VAE is evaluated in a variety of experiments where it outperforms alternative approaches including multi-output GPs and structured VAEs.

Under review as a conference paper at ICLR 2021.
SPARSE GAUSSIAN PROCESS VARIATIONAL AUTOEN- CODERS
Anonymous authors Paper under double-blind review
ABSTRACT
Large, multi-dimensional spatio-temporal datasets are omnipresent modern sci- ence and engineering. An effective framework for handling such data are Gaus- sian process deep generative models (GP-DGMs), which employ GP priors over the latent variables of DGMs. Existing approaches for performing inference in GP-DGMs do not support sparse GP approximations based on inducing points, which are essential for the computational efficiency of GPs, nor do they handle missing data - natural occurrence in many spatio-temporal datasets in a prin- cipled manner. We address these shortcomings with the development of the sparse Gaussian process variational autoencoder (SGP-VAE). characterised by the use of partial inference networks for parameterising sparse GP approximations. Lever- aging the benefits of amortised variational inference, the SGP-VAE enables infer- ence in multi-output sparse GPs on previously unobserved data with no additional training. The SGP-VAE is evaluated in variety of experiments where it outper- forms alternative approaches including multi-output GPs and structured VAEs.
INTRODUCTION
Increasing amounts of large, multi-dimensional datasets that exhibit strong spatio-temporal depen- dencies are arising from wealth of domains, including earth, social and environmental sciences (Atluri et al., 2018). For example, consider modelling daily atmospheric measurements taken by weather stations situated across the globe. Such data are (I) large in number; (2) subject to strong spatio-temporal dependencies; (3) multi-dimensionall and (4) non-Gaussian with complex depen- dencies across outputs. There exist two venerable approaches for handling these characteristics: Gaussian process (GP) regression and deep generative models (DGMs). GPs provide a frame- work for encoding high-level assumptions about latent processes, such as smoothness or periodicity, making them effective in handling spatio-temporal dependencies. Yet, existing approaches do not support the use of flexible likelihoods necessary for modelling complex multi-dimensional outputs. In contrast, DGMs support the use of flexible likelihoods; however, they do not provide natural route through which spatio-temporal dependencies can be encoded. The amalgamation of GPs and DGMs, GP-DGMs, use latent functions drawn independently from GPs, which are then passed through a DGM at each input location. GP-DGMs combine the complementary strengths both approaches, making them naturally suited for modelling spatio-temporal datasets.
Intrinsic to the application of many spatio-temporal datasets is the notion of tasks. For instance: medicine has individual patients; each trial in a scientific experiment produces individual dataset; and, in the case of a single large dataset, it is often convenient to split it into separate tasks to im- prove computational efficiency. GP-DGMs support the presence of multiple tasks in a memory efficient way through the use of amortisation, giving rise to the Gaussian process variational au- toencoder (GP- VAE), model that has recently gained considerable attention from the research community (Pearce, 2020; Fortuin et al., 2020; Casale et al., 2018; Campbell & Liò, 2020; Ram- chandran et al., 2020). However, previous work does not support sparse GP approximations based on inducing points, necessity for modelling even moderately sized datasets. Furthermore, many datasets contain an abundance of missing data: weather measurements are often absent due to sensor failure, and in medicine only single measurements are taken at any instance. Handling partial observations in a principled manner essential for modelling spatio-tempont data, but is yet to be considered.

---

Under review as a conference paper at ICLR 2021.
Our key technical contributions as follows:
i) We develop the sparse GP- VAE (SGP- (SGP-VAE), which uses inference networks to parameterise multi-output sparse GP approximations.
ii) We employ a suite partial inference networks for handling missing data the SGP-VAE.
iii) We conduct rigorous evaluation of the SGP-VAE in variety of experiments, demonstrat- ing excellent performance relative to existing multi-output GPs and structured VAEs.
A FAMILY OF SPATIO-TEMPORAL VARIATIONAL AUTOENCODERS
Consider the multi-task regression problem in which we wish to model T datasets {D(1)) J1=1+ each of which comprises input/output pairs ( fxx (f)
In=1, Xn E RR and (1) E Rp. Fur- ther, let any possible permutation of observed values be potentially missing, such that each observa- tion yn y U Uyn" contains a set of observed values yn and unobserved values (!),), with ol denoting the index set of observed values. For each task, model the distribution of each ob- servation yn conditioned on a corresponding latent variable { E RK as a fully-factorisee Gaus- sian distribution parameterised by passing f through decoder deep neural network (DNN) with parameters 02. The elements of f(() correspond to the evaluation of K K--dimensional latent function at input x") That is, f(()) = (()) (xx Each latent function f() is modelled as being drawn from K independent GP priors with giving rise to the complete probabilistic model:
f(() = (f(t) ((t) ...., f(t
01 = 10 1K
K
p(
II gP (0, Kor, (x,x'))
y0
II (12 (P( diag 0g ((r(o)
(
Poy (y2 1800 xn (o(
((F(") are the outputs the decoder indexed by O We shall refer to the set = (0,, }2} as the model parameters, which are shared across tasks. The probabilistic model in equation explicitly accounts for dependencies between latent variables through the GP prior. The motive of the latent structure is twofold: to discover a simpler representation each observation, and to capture the dependencies between observations at different input locations.
where peg, (f(!) and og 2
2.1 MOTIVATION FOR SPARSE APPROXIMATIONS AND AMORTISED INFERENCE
The use of amortised inference in DGMs and sparse approximations in GPs enables inference in these respective models scale large quantities of data. To ensure the same for the GP-DGM described in equation 1, we require the use of both techniques. In particular, amortised inference is NNt)..
necessary to prevent the number of variational parameters scaling with ( Further,
the inference network can be used condition previously unobserved data without needing to learn new variational parameters. Similarly, sparse approximations are necessary to prevent the computational complexity increasing cubically with the size of each task O (E (EL,NO)))
Unfortunately, it is far from straightforward to combine sparse approximations and amortised infer- ence computationally efficient way. To see this, consider the standard form for the sparse GP approximate posterior, q(f) pa, (f (ulluuq(j) where qu) N (u; m, S) with m. S and Z, the inducing point locations, being the variational parameters. (u)) does not decompose into a product over (e) factors and is therefore not amendable to per-datapoint amortisation. That is, m and S must be optimised as free-form variational parameters. naive approach to achieving per-datapoint amortisation is to decompose q(u) into the prior Po, (u) multiplied by the product of approximate likelihoods, one for each inducing point. Each approximate likelihood is itself equal to the product of per-datapoint approximate likelihoods, which depend on both the observation y", and the distance of the input x,, to that of the inducing point. An inference network which takes these two values of inputs be used to obtain the parameters of the approximate likelihood factors. Whilst we found that this approach worked, is somewhat umprincipled. Moreover, it requires passing each data- point/inducing point pair through an inference network, which scales very poorly. In the following
2

---

Under review as a conference paper at ICLR 2021.
section, we introduce a theoretically principled decomposition of q(f) term the sparse structured approximate posterior which will enable efficient amortization.
2.2 THE SPARSE STRUCTURED APPROXIMATE POSTERIOR
By simultaneously leveraging amortised inference and sparse GP approximations, we can perform efficient and scalable approximate inference. We specify the sparse structured approximate poste- rior, (()()),, which approximates the intractable true posterior for task t:
N II
po(s x(0) = (100))
Zp por (5(0)
Pod (y
(1
15(0 x(0 1000)
(2)
N
Lo yu 2 Po. (
II
xx 2) - q(5lo)
Analogous to its presence in the true posterior, the approximate posterior retains the GP prior, yet replaces each non-conjugate likelihood factor with approximate likelihood, Lov yu (1) over set of KM 'inducing points', Zmk- For tractability, we restrict the approximate likelihoods to be Gaussians factorised across each latent dimension, parameterised by passing each observation through partial inference network:
Xn (1) Z),
- U Wmak 'inducing locations' Z = Lk U
L0 (us yn
X (1) Zx) =N (Mon.LYY
kydt pu K ug, o (y
where Φi denotes the weights and biases of the partial inference network, whose outputs are shown in red. This form is motivated by the work of Bui et al. (2017), who demonstrate the optimality of approximate likelihoods of the form N (gn k 520 K. uk pn): result we prove in Ap- pendix A.I. Whilst, in general, the optimal free-form values gn. and vn, depend on all of the data points, we make the simplifying assumption that they depend only (t), For GP regression with Gaussian noise, this assumption holds true as In = yv and U, o (Bui et al., 2017).
The resulting approximate posterior can be interpreted the exact posterior induced by a surrogate regression problem, in which "pseudo-observations' In are produced from a linear transformation of inducing points with additive 'pseudo-noise" Un. In = ,,tt u. uguu Uy Vone. The inference network learns to construct this surrogate regression problem such that results in posterior that is close to our target posterior.
By sharing variational parameters o {or, Z} across tasks, inference amortised across both datapoints and tasks. The approximate posterior for single task corresponds to the product of K independent GPs, with mean and covariance functions
(x) = k K us
D
(()) (t) K Kergy bikk Mork
(4)
2.00) (x, x') p(
R2 - k & us
K-1
(1) uu uuue kunf + k k us k kun f0
where $ (1) Kujui
K
(()) K Kepun [
(0)
(y and de.k ij
p(
8yj A (0). ())(()) See Appendix A.2 for a complete derivation. The computational complexity asso- ciated with evaluating the mean and covariance functions is O (KM)N()) significant improve- ment over the O (paN( cost associated with exact multi-output GPs for KM² < p0) We refer to the combination of the aforementioned probabilistic model and sparse structured ap- proximate posterior as the SGP-VAE.
The SGP- VAE addresses three major shortcomings of existing sparse GP frameworks. First, the inference network can be used to condition on previously unobserved data without needing to learn new variational parameters. Suppose we use the standard sparse GP variational approximation q(f) = pe, (f,,, (,)(()((( where q(u) A (u; m, S). If more data are observed, m and S have to be re-optimised. When an inference network is used to parameterise quu), the approximate posterior
3

---

Under review as a conference paper at ICLR 2021.
is "automatically" updated by mapping from the new observations to their corresponding approxi- mate likelihood terms. Second, the complexity of the approximate posterior can be modified as desired with no changes to the inference network, or additional training, necessary: any changes in the morphology of inducing points corresponds to a deterministic transformation of the inference network outputs. Third, if the inducing point locations parameters does not depend on the size of the dataset, even as more inducing points are added. This contrasts with the standard approach, in which new variational parameters are appended m and S.
fixed, then the number of variational
2.3 TRAINING THE SGP- VAE
Learning and inference in the SGP-VAE are concerned with determining the model parameters 0 and variational parameters Φb These objectives can be attained simultaneously by maximising the evidence lower bound (ELBO), given by LELBO Ei=
ELBO where
((t)
po(y(0,,100) (((())
(1) ELBO =Bggrin)
E (erttt
KL (q) ll peo pee (u))
(5)
and ( x Po, (u) II Le (u: ( X Z). Fortunately, since both (u) and pa, (u) are multivariate Gaussians, the final term, and its gradients, has an analytic solution. The first term amounts to propagating a Gaussian through non-linear DNN, so must be approximated using a Monte Carlo estimate. We employ the reparameterisation trick (Kingma & Welling, 2014) account for the dependency of the sampling procedure on both and @ when estimating its gradients.
(t) We mini-batch over tasks, such that only single ELBO is computed per update. Importantly in combination with the inference network, this means that avoid having to retain the O (TM*) terms associated with T Cholesky factors if we were to use free-form q(u) for each task. Instead, the memory requirement is dominated by the O OKM2² + KNN + (rl) terms associated with (1) storing Ku,u. K p(t) and di, as instantiating & and (t) u.
b involves only O (KN) terms.' This corresponds to considerable reduction in memory. See Appendix C for thorough comparison of memory requirements.
2.4 PARTIAL INFERENCE NETWORKS
Partially observed data is regularly encountered in spatio-temporal datasets, making it necessary to handle it a principled manner. Missing data is naturally handled by Bayesian inference. However, for models using inference networks. it necessitates special treatment. One approach to handling partially observed data is to impute missing values with zeros (Nazabal et al., 2020; Fortuin et al., 2020). Whilst simple to implement, zero imputation is theoretically unappealing as the inference network can no longer distinguish between a missing value and true zero.
Instead, we turn towards the ideas of Deep Sets (Zaheer et al., 2017). By coupling the observed value with dimension index, we may reinterpret each partial observation as a permutation invariant set. We define family of permutation invariant partial inference networks' as
(Hg(Y)), log (Dr)) = P4z
Σ ho, (smp) PEO..
(6)
where her 122 R and Ps. R 2P are DNN mappings with parameters di and Φ2. respectively. Sap denotes the couples of observed value Yup and corresponding dimension index p. The formulation in equation 6 is identical to the partial variational autoencoder (VAE) framework established by Ma et al. (2019). There are number of partial inference networks which conform to this general framework, three of which include:
This assumes input locations are shared across tasks, which true for all experiments we considered. Whilst the formulation in equation 6 can describe any permutation invariant set function, there is caveat: both hag and Po. be infinitely complex, hence linear complexity is guarantced.
4

---

Under review as a conference paper at ICLR 2021.
PointNet Inspired by the PointNet approach of Qi et al. (2017) and later developed by Ma et al. (2019) for use in partial VAEs, the PointNet specification uses the concatenation of dimen- sion index with observed value: Sup (p. Bnp). This specification treats the dimension indices as continuous variables. Thus, an implicit assumption of PointNet is the assumption of smoothness between values of neighhouring dimensions. Although valid in computer vision application, it is ill-suited for tasks in which the indexing of dimensions is arbitrary.
IndexNet Alternatively, one may use the dimension index to select the first DNN mapping: hes (Smp) = heg. (Ump). Whereas PointNet treats dimension indices as points in space, this specification retains their role as indices. We refer to it as the IndexNet specification.
FactorNet A special case of IndexNet, first proposed by Vedantam et al. (2017), uses sepa- rate inference network for each observation dimension. The approximate likelihood is fac- torised into product of Gaussians, one for each output dimension: L,y (ux; yū: Xn, Zs)
Ipeo, N/ (Hab. (Hap); Kfann K Kugng Ug, &orpA (UmP)) We term this approach FactorNet.
See Appendix G for corresponding computational graphs. Note that FactorNet is equivalent to IndexNet with Pe, defined by the deterministic transformations of natural parameters of Gaussian distributions. Since IndexNet allows this transformation to be learnt, we might expect it to always produce better partial inference network for the task at hand. However, it is important to consider the ability of inference networks to generalise. Although more complex inference networks will better approximate the optimal non-amortised approximate posterior on training data, they may produce poor approximations
it on the held-out data.² In particular, FactorNet is constrained to consider the individual contribution of each observation dimension, whereas the others are not. Doing so is necessary for generalising to different quantities and patterns of missingness, hence anticipate FactorNet to perform better in such settings.
3 RELATED WORK
We focus our comparison approximate inference techniques; however, Appendix D presents unifying view of GP-DGMs.
Structured Variational Autoencoder Only recently has the use of structured latent variable priors in VAEs been considered. In their seminal work, Johnson et al. (2016) investigate the combination of probabilistic graphical models with neural networks to learn structured latent variable represen- tations. The authors consider a two stage iterative procedure, whereby the optimum of a surrogate objective function containing approximate likelihoods in place of true likelihoods - is found and substituted into the original BLBBO. The resultant structured VAE (S VAE) objective is then op- timised. In the case of fixed model parameters 0, the SVAE objective is equivalent to optimising the ELBO using the structured approximate posterior over latent variables (z) x po ( Accordingly, the SGP-VAE can be viewed as an instance of the SVAE. Lin al. (2018) build upon the SVAE, proposing a structured approximate posterior of the form g(z) α qa 1o (z)y). The authors refer to the approximate posterior as the structured inference network (SIN). Rather than using the latent prior Pa (z), SIN incorporates the model's 's latent structure through (z). The core advantage of SIN is its extension to more complex latent priors containing non-conjugate factors qs (z) can replace them with their nearest conjugate approximations whilst retaining a similar latent structure. Although the frameworks proposed by Johnson et al. and Lin et al. are more general than ours, the authors only consider Gaussian mixture model and linear dynamical system (LDS) latent priors
Gaussian Process Variational Autoencoders The earliest example of combining VAEs with GPs is the GP prior VAE (GPPVAB) (Casale et al., 2018). There are significant differences between our work and the GPPVAE, most notably in the GPPVAE's use of a fully-factorised approximate posterior an approximation that is known to perform poorly in time-series and spatial settings (Turner & Sahani, 2011). Closely related the GPPVAE Ramchandran al.'s 's (2020) longitu- dinal VAF, which also adopts factoriat approximate posterior, yet uses additive covariance functions for heterogeneous input data. Fortuin et al. (2020) consider the use of Gaussian ap-
This kind of 'overfitting' is different to the classical notion of overfitting model parameters. It refers to how well optimal non-amortised approximate posteriors are approximated the training versus test data.
5

---

Under review as a conference paper at ICLR 2021.
proximate posterior with tridiagonal precision matrix A, q(f) N (f: m, A-),, where and A are parameterised by an inference network. Whilst this permits computational efficiency, the parameterisation is only appropriate for regularly spaced temporal data and neglects rigorous treat- ment of long term dependencíes. Campbell & Liò (2020) employ an equivalent sparsely structured variational posterior that used by Fortuin et al., extending the framework to handle more general ppatio-temporal data. Their method is similarly restricted to regularly spaced spatio-temporal data. fundamental difference between our framework and that of Fortuin et al. and Campbell & Liò is the inclusion of the GP prior in the approximate posterior. As shown by Opper & Archambeau (2009), the structured approximate posterior identical in form the optimum Gaussian approxi- mation to the true posterior. Most similar to ours is the approach of Pearce (2020), who considers the structured approximate posterior q(f) 2, I[= (0 yn). We refer to this as the GP- PAE. Pearce's approach special case of the SGP-VAE for u f and no missing data. Moreover, Pearce only considers the application to modelling pixel dynamics and the comparison to the standard VAF. See Appendix B for further details.
4 EXPERIMENTS
We investigate the performance of the SGP- illustrative bouncing ball experiments, followed by experiments in the small and large data regimes. The first bouncing ball experiment provides visualisation of the mechanics of the -GP-VAE, and a quantitative comparison to other structured VAEs. The proceeding scall-scale experiments demonstrate the utility of the GP-VAE and show that amortisation, especially in the presence of partially observed data, is not at the expense of predictive performance. In the final two experiments, showease the efficacy of the SGP-VAE on large, multi-output apatio-tempors datasets for which the use of amortisation is necessary. Full experimental details are provided Appendix E.
4.1 SYNTHETIC BOUNCING BALL EXPERIMENT
The bouncing ball experiment first introduced by Johnson et al. (2016) for evaluating the SVAE and later considered by Lin et al. (2018) for evaluating SIN considers a sequence of one- dimensional images of height 10 representing ball bouncing under linear dynamics, (xh, R 1! (O). The VAE is able to significantly outperform both the SVAE and SIN in the origi- nal experiment, shown in Figure 1a. To showcase the versatility of the SGP-VAE complexity of the original experiment to consider a sequence of images of height 100, y 10000 representing two bouncing balls: one under linear dynamics and another under gravity. Furthermore, the images are corrupted by removing 25% of the pixels at random. The dataset consists of 80 noisy image sequences, each of length N 500, with the goal being to predict the trajectory of the ball given prefix of a longer sequence.
extend the
LDS
SVAE
SIN
IGP
GP-VAE
Mean prsdiction STmur
1.0
SGRAAE
06
R(51 3g
0.
Postrroe bnctt GP
<00
30
10
100
303
40
E-siqp alocal
Tim:
(a)
(b)
Figure 1: a) Comparing the GPE's predictive performance to that of the SVAE, SIN, an LDS and independent GPs (IGP). b) Top: sequence of images representing two bouncing balls. Middle: mean of the SGP-VAEL predictive distribution conditioned on partial observations up to the red line. Bottom: latent approximate GPs posterior, alongside the inducing point locations (crosses).
Using two-dimensional latent space with periodic kernels, Figure 1b compares the posterior latent GPs and the mean predictive distribution with the ground truth for a single image sequence. Observe that the SGP-VAE has 'disentangled" the dynamics of each ball, using single latent dimension

---

Under review as a conference paper at ICLR 2021.
model each. The SGP- -AE reproduces the image sequences with impressive precision, owing in equal measure to (1) the ability of the GPs prior to model the latent dynamics and (2) the flexibility of the likelihood function map to the high-ding high-dimennsioal observations.
4.2 SMALL-SCALE EXPERIMENTS
EEG Adopting the experimental procedure laid out by Requeima et al. (2019), we consider an EEG dataset consisting of N = 256 measurements taken over a one second period. Each measurement comprises voltage readings taken by seven electrodes. FZ and F1-F6, positioned on the patient's scalp (x, E R yn E R ). The goal is to predict the final 100 samples for electrodes FZ, F1 and F2 having observed the first 156 samples, as well as all 256 samples for electrodes F3-F6
Jura The Jura dataset is a geospatial dataset comprised of N = 359 measurements of the topsoil concentrations of three heavy metals Cadmium Nickel and Zinc- - collected from 14.5km² region of the Swiss Jura (x, E RR2 Nn E IRS)) (Goovaerts, 1997). Adopting the experimental procedure laid out by others (Goovaerts, 1997; Álvarez & Lawrence, 2011; 1; Requeima et al., 2019), the dataset is divided into a training set consisting of Nickel and Zinc measurements for all 359 locations and Cadmium measurements for just 259 locations. Conditioned on the observed training set, the goal is to predict the Cadmium measurements at the remaining 100 locations.
Table 1: A comparison between multi-output GP models on the EEG and Jura experiments.
GP-VAE
Metric
IGP!
GPAR!
ZI
PointNet
IndexNet
FactorNet
FF
SMSE
1.75
0.26
0.27 (0.03)
0.60 (0.09)
0.24 (0.02)
0.28 (0.04)
0.25 (0.02)
EEG
NLL
2.60
1.63
2.24 (0.37)
3.03 (1.34)
2.01 (0.28)
2.23 (0.21)
2.13 (0.28)
MAE NLL
0.57
0.41
0.42 (0.01)
0.45 (0.02)
0.44 (0.02)
0.40 (0.01)
0.41 (0.01)
Jura
1.13 (0.09)
1.13 (0.08)
1.12 (0.08)
1.00 (0.06)
1.04 (0.06)
Results taken directly from Requeima et al. (2019).
Table compares the performance of the VAE using the three partial inference networks pre- sented in Section 2.4, as well as zero imputation (ZI), with independent GPs (IGP) and the GP autoregressive regression model (GPAR), which, our knowledge, has the strongest published per- formance on these datasets. We also give the results for the best performing GP- VAE using a non-amortised, or 'free-form' (FF), approximate posterior, with model parameters 9 kept fixed to the optimum found by the amortised GP- VAE and variational parameters initialised to the output of the optimised inference network. All GP-VAE models use a two- and three-dimensional latent space for EEG and Jura, respectively, with squared exponential (SE) kernels. The results highlight the poor performance of independent GPs relative to multi-output GPs, demonstrating the impor- tance of modelling output dependencies. The GP-VAE achieves impressive SMSE and MAE on the EEG and Jura datasets using all partial inference networks except for PointNet. In Appendix F demonstrate superior performance of the GP-VAE relative the GPPVAE, which can be at- tributed to the use of the structured approximate posterior over the mean-field approximate posterior used by the GPPVAE. Importantly, the negligible difference between the results using free-form and amortised approximate posteriors indicates that amortisation is not at the expense of predictive performance.
Whilst GPAR performs as strongly as the GP- EAP in the small-scale experiments above, it has two key limitations which severely limit the types of applications where it can be used. First, it can only be used with specific patterns of missing data and not when the pattern of missingness is arbitrary. Second, it is not scalable and would require further development to handle the large datasets con- sidered in this paper. In contrast, the SGP- VAE is far more flexible: it handles arbitrary patterns of missingness, and scales to large number of datapoints and tasks. A distinct advantage of the SGP- VAE is that it models P outputs with just K latent GPs. This differs from GPAR, which uses P
'i. using Net for EEG and FactorNet for Jura.
'The two different performance metrics adopted enable comparison the results of Requeima al..
7

---

Under review as a conference paper at ICLR 2021.
GPs. Whilst this is not issue for the small-scale experiments, it quickly becomes computation- ally burdensome when P becomes large. The true efficacy of the SGP- VAE is demonstrated in the following three experiments, where the number of datapoints and tasks is large, and the patterns of missingness are random.
4.3 LARGE-SCALE EEG EXPERIMENT
We consider an alternative setting to the original small-scale EEG experiment, in which the datasets formed from T 60 recordings of length N 256, each with 64 observed voltage readings (yn E IR)4) For each recording, we simulated electrode "blackouts' by removing consecutive sam- ples random. We consider two experiments: in the first, we remove equal 50% of data from both the training and test datasets; in the second, we remove 10% of data from the training dataset and 50% from the test dataset. Both experiments require the partial inference network generalise to different patterns of missingness, with the latter also requiring generalisation to different quantities of missingness. Each model trained on 30 recordings, with the predictive performance assessed on the remaining 30 recordings. Figure 2 compares the performance of the SGP- VAE with that of independent GPs as the number of inducing points varies, with = 256 representing use of 10-dimensional latent space with SE kernels. The SGP-VA VAE using PointNet results in substantially worse performance than the other partial inference networks, achieving average SMSE and NLL of 1.30 and 4.05 on the first experiment for M 256. Sim- ilarly, using a standard VAE results in poor performance, achieving an average SMSE and NLL of 1.62 and 3.48. These results excluded from Figure for the sake of readability.
the GP-VAE. In each case,
+ Z1
IndexNet
FactorNet
NLL
Training time
IIP
Experiment
1000
3.00
0.50
2.75
JGP
Expeiment
1.00
1500
3.25
0.75
3,00
ICP
500
256
Figure 2: Variation in performance of the SGP-VAE on the large-scale EEG experiment as the number of inducing points, M, varies.
For all partial inference networks, the SGP- VAE achieves a significantly better SMSE than inde- pendent GPs in both experiments, owing to its ability to model both input and output dependencies. For the first experiment, the performance using FactorNet is noticeably better than using either In- dexNet or zero imputation; however, this comes at the cost of a greater computational complexity associated with learning inference network for each output dimension. Whereas the performance for the SGP- VAE using IndexNet and zero imputation significantly worsens on the second exper- iment, the performance using FactorNet is comparable the only partial inference network that is able to accurately quantify the contribution of each output dimension to the latent posterior, enabling it generalise different quantities of missing data.
the first experiment. This suggests it is
The advantages using a sparse approximation are clear - using M = 128 inducing points results in slightly worse average SMSE and NLL., yet significantly less computational cost.
4.4 JAPANESE WRATHER EXPERIMENT
Finally, consider dataset comprised of 731 daily climate reports from 156 Japanese weather stations throughout 1980 and 1981, total of 114,036 multi-dimensional observations. Weather reports consist of date and location, including elevation, alongside the day's maximum, minimum and average temperature, precipitation and snow depth (xf E R2, yn E IR²), any number of which is potentially missing. We treat each week single task, resulting in T 105 tasks with
8

---

Under review as a conference paper at ICLR 2021.
Day
Day
Day
Day
Day
Day
Day
Tan
Trass
Tag
Figure 3: An illustration of the Japanese weather experiment. The dotted red lines highlight the missing data, with the CGP-VAE's predictive mean shown below.
N 1092 data points each. The goal is to predict the average temperature for all stations on the middle five days, illustrated in Figure 3. Each model is trained on all the data available from 1980. For evaluation, we use data from both 1980 and 1981 with additional artificial missingness - the average temperature for the middle five days and a random 25% of minimum and maximum temperature measurements". Similar to the second large-scale EEG experiment, the test datasets have more missing data than the training datasets. Table 2 compares the performance of the SGP- VAE using 100 inducing points to that of a standard VAE using FactorNet and a baseline of mean imputation. All models three-dimensional latent space with SE kernels.
Table 2: A comparison between model performance on the Japanese weather experiment.
CGP-VAE
Metric
MI
IGP
VAE
Z1
PointNet
IndexNet
FactorNet
RMSE
9.16
2.49 (0.00)
2.31 (0.35)
3.21 (0.35)
3.53 (0.13)
2.59 (0.11)
1.60 (0.08)
1980
NLL
2.84 (0.00)
2.68 (0.15)
2.68 (0.15)
2.80 (0.05)
2.44 (0.04)
2.18 (0.02)
RMSE
9.27
2.41 (0.00)
3.20 (0.37)
3.20 (0.37)
3.61 (0.14)
2.55 (0.11)
1.68 (0.09)
1981
NLL
2.73 (0.00)
2.67 (0.15)
2.67 (0.15)
2.83 (0.06)
2.43 (0.04)
2.20 (0.02)
All models significantly outperform the mean imputation baseline (MI) and are able to generalise inference to the unseen 1981 dataset without any loss in predictive performance. The SGP-VAE achieves better predictive performance than both the standard VAE and independent GPs, showcas- ing its effectiveness in modelling large spatio-temporal datasets. The SGP-VAE using FactorNet achieves the best predictive performance on both datasets. The results indicate that FactorNet is the only partial inference network capable of generalising to different quantities and patterns of missingness, supporting the hypothesis made in Section 2.4.
CONCLUSION
The SGP-VAI is scalable approach to training GP-DGMs which combines sparse inducing point methods for GPs and amortisation for DGMs. The approach is idcally suited to spatio-temporal data with missing observations, where outperforms VAEs and multi-output GPs. Future research directions include generalising the framework to leverage state-space GP formulations for additional scalability and applications streaming multi-output data.
REFERENCES
Mauricio A Álvarez and Neil D Lawrence. Computationally efficient convolved multiple output Gaussian processes. The Journal of Machine Learning Research, 12: 12:0409-1500, 2011.
'This prevents the model from using solely output dependencies impute the missing temperature readings.

---

Under review as a conference paper at ICLR 2021.
Gowtham Atluri, Anuj Karpatne, and Vipin Kumar. Spatio-temporal data mining: A survey of problems and methods. ACM Computing Surveys (CSUR), 51(4):1-41, 2018.
Edwin V Bonilla, Kian M Chai, and Christopher Williams. Multi-task Gaussian process prediction. In Advances Neural Information Processing Systems, pp. 153-160, 2008.
Thang D Bui, Josiah Yan, and Richard Turner. unifying framework for Gaussian process pseudo-point approximations using power expectation propagation. The Journal of Machine Learning Research, 184433364, 2017.
Alex Camphell and Pietro Liò. tvGP- VAE: Tensor-variate Gaussian process prior variational au- toencoder. arXiv preprint arX1v:2006.04788, 2020.
Francesco Casale, Adrian Dalca, Luca Saglietti, Jennifer Listgarten, and Nicolo Fusi. Gaussian process prior variational autoencoders. In Advances in Neural Information Processing Systems, pp. 10369-10380, 2018.
Andreas Damianou and Neil Lawrence. Deep Gaussian processes. In Artificial Intelligence and Statistics, pp. 207-215, 2013.
Vincent Fortuin, Dmitry Baranchuk, Gunnar Rätsch, and Stephan Mandt. GP-VAE: Deep probabilis- tic time series imputation. In International Conference on Artificial Intelligence and Statistics, pp. 1651-1661, 2020.
Pierre Goovaerts. Geostatistics for natural resources evaluation. Oxford University Press on De- mand, 1997.
Matthew Johnson, David K Duvenaud, Alex Wiltschko, Ryan Adams, and Sandeep R Datta. Composing graphical models with neural networks for structured representations and fast infer- ence. In Advances in Neural Information Processing Systems, pp. 2946-2954, 2016.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. International Conference on Learning Representations, 12 2014.
Diederik P Kingma and Max Welling. Auto-encoding variational Bayes. In 2nd International Conference on Learning Representations, 2014.
Neil D Lawrence. Gaussian process latent variable models for visualisation of high dimensional data. In Advances in Neural Information Processing Systems, pp. 329-336, 2004.
Neil D Lawrence and Andrew Moore. Hierarchical Gaussian process latent variable models. In Proceedings of the 24th International Conference Machine Learning, pp. 481-488, New York, NY, USA, 2007. Association for Computing Machinery.
Wu Lin, Nicolas Hubacher, and Mohammad Emtiyaz Khan. Variational message passing with struc- tured inference networks. In International Conference on Learning Representations, 2018.
Chao Ma, Sebastian Tschiatschek, Konstantina Palla, José Miguel Hernández-Lobato. Sebastian Nowozin, and Cheng Zhang. EDDI: efficient dynamic discovery of high-value information with partial VAE. In Proceedings of the 36th International Conference on Machine Learning, vol- ume 97, pp. 4234 4243, 2019.
Alfredo Nazabal, Pablo M Olmos, Zoubin Ghahramani, and Isabel Valera. Handling incomplete heterogeneous data using VAEs. Pattern Recognition, pp. 107501, 2020.
Manfred Opper and Cédric Archambeau. The variational Gaussian approximation revisited. Neural computation, 21(3):786-792, 2009.
Michael Pearce. The Gaussian process prior VAE for interpretable latent dynamies from pixels. In Symposium on Advances in Approximate Bayesian Inference, pp. 1-12, 2020.
Charles R Qi, Hao Su, Kaichun Mo, and Leonidas J Guibas. PointNet: Deep learning on point sets for 3D classification and segmentation. Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 652-660, 2017.
10

---

Under review as a conference paper at ICLR 2021.
Siddharth Ramchandran, Gleb Tikhonov, Miik Koskinen, and Harri Lähdesmäki. Longitudinal variational autoencoder. arXiv preprint arXiv:2006.0976. 09763, 2020.
James Requeima, William Tebbutt, Wessel Bruinsma, and Richard E Turner. The Gaussian process autoregressive regression model (GPAR). In The 22nd International Conference on Artificial Intelligence and Statistics, pp. 1860-1869, 2019.
Yee Whye Teh, Matthias W Seeger, and Michael Jordan. Semiparametric latent factor models. In Proceedings of the Tenth International Workshop on Artificial Intelligence and Statistics, 2005,
Michael E Tipping and Christopher M Bishop. Probabilistic principal component analysis. Journal of the Royal Statistical Sociery: Series B (Statistical Methodology), 61(3):611 1--622 1999.
Michalis Titsias. Variational learning inducing variables in sparse Gaussian processes. In Artificial Intelligence and Statistics, pp. 567-574, 2009.
Rich Turner and Maneesh Sahani. Two problems with variational expectation maximisation for time-series models. In D. Barber, T. Cemgil, and S. Chiappa (eds.), Bayesian Time series models, chapter 5, pp. 109-130. Cambridge University Press, 2011.
Ramakrishna Vedantam, Ian Fischer, Jonathan Huang, and Kevin Murphy, Generative models of visually grounded imagination. arXiv preprint arXiv:1705.10780 2017.
Byron M Yu, John P Cunningham, Gopal Santhanam, Stephen Ryu, Krishna V Shenoy, and Ma- Gaussian-process
neesh Sahani. Gaussian process factor analysis for diw-dimensional single-trial analysis of neural population activity. In Advances in Neural Information Processing Systems, pp. 1881 1881-1888, 2009.
Manzil Zaheer, Satwik Kottur, Siamak Ravanbakhsh, Barnabas Poczos, Russ R Salakhutdinov, and Alexander J Smola. Deep sets. In Advances in Neural Information Processing Systems, pp- 3391-3401,2017
A MATHEMATICAL DERIVATIONS
A.1 OPTIMALITY OF APPROXIMATE LIKELIHOODS
To simplify notation, shall consider the case 1 and K 1. Separately, Opper & Archam- beau (2009) considered the problem of performing variational inference in a GP for non-Gan non-Gaussian likelihoods. They consider multivariate Gaussian approximate posterior, demonstrating that the optimal approximate posterior takes the form
((f) - P(r) II N ( Ins un),
requiring a total of 2NN variational parameters ((gn, vn 1.
In this section, derive result that generalises this inducing point approximations, showing that for fixed M the optimal approximate posterior can be represented by max [M (M +1)/2+ 2N). Following Titsias (2009), we consider an approximate posterior the form
q(f) q(u)p(f(,lu))
where qu)) = N mu. Kau) is constrained to be a multivariate Gaussian with mean m, and covariance kmn. The ELBO is given by
LELBO
Ey(rr [logp(y[f)] KL (q(u) || p(u)
Eg(i)) [Ep(f(n) [logp(y]f]]] KL (q(u) puu)
(9)
N
Σ Bg(u) EN(Ja: A.utan K,,1a) [logp(y,(fn]] KL (g(u) (u))
n=l

---

Under review as a conference paper at ICLR 2021.
where
A, = K fr K mn
(10)
On m. K f.u K mu.
(11)
Recall that for twice-differentiable scalar function h
Vy p. 2) [h(w]] EN( p. 2) [u)
(12)
where Ha (u) is the Hessian of h at u. Thus, the gradient the ELBO with respect to K.u can be rewritten as
N
I 1 Kun
Ruu
LELBO Σ
N(u tha Ruo)
Hnn (u)]
Kuu +
(13)
2
2
where (u) = E W Anutan, k)lu)) plogply,lfll
To determine an expression for Hnn, we first consider the gradients of hn. Let
an (Bn) (N &n. Ky,,,.) [ly.p(y,,..n))
(14)
Bn (u) = A, + an.
(15)
The partial derivative of hn with respect to the jth element of u can expressed
Ōhn (u) Oa, (3, (u)
(u).
(16)
Ouj 03,
Du,
Taking derivatives with respect to the th element u gives
hn (u)
an
0Bn (u) (u) + 03n
3gn
dan
0
(8n (u))
(3, (u))
(u).
(17)
Ou 032
Ou, Ou:
OuyJu,
Thus, the Hessian is given by
pan
dan (Bn (u) Ho. (u) 03,
H (u) 052
(pn (u))
VB (u) (u)
(18)
Nx1
1xN
Since (u)
A, + ans we have (B, = An and H3. (u) = O. This allows us to write
V K... LELRO as
Pan (Bn(w)) A.A. 1 Kin 032
1 k
V Ru, CELBO
Σ
Ex(u: aa. Rmn)
(19)
2
The optimal covariance therefore satisfies
[oan (Sn(w)) A,A 083
K
K- 2 Σ
Ex(m: ia. Rmn)
(20)
Similarly, the gradient of the ELBO with respect to i, can be written as
LELBO
V iun E N(u: m.. R.n) nn, (u)] K K (m. mu)
m.
(21)
N
Σ
E E ma. Kun) [u)( K-1 (m, m.)
where we have used the fact that for differentiable scalar function
V Bu Σ) [g(u)] EN Σ) [Vg(u)]
(22)
12

---

Under review as a conference paper at ICLR 2021.
Using equation 16 and (,(u)= An + an get
dan (Bn(u)))) (n)) 03,
Vhn(u)
(23)
giving
N
Jan (Sn(u)) 05,
V Via., LELBO =
Σ
B N(u mu. Ruu)
K mu).
(24)
The optimal mean is therefore
N Bui(u: ria.. Run)
0an (Bn(u)) 03,
m mu
KunAn.
(25)
Equation 20 and equation 25 show that each 1h observation contributes only rank-1 1 term to the optimal approximate posterior precision matrix, corresponding to an optimum approximate posterior of the form
(f) x p(f) II N (K fa uuuuu In. Un) ==1
(26)
where
Oan (Bn( 03n
In
E (w: mar Ruu)
kn Kun + A m,
(27)
O Jaan an (3n, (u))
1/0, -2E 22 N(u
(28)
mun Run)
0,3
For general likelihoods, these expressions cannot be solved exactly so Ja and en are freely optimised as variational parameters. When N M, the inducing points are located at the observations and A, A is zero everywhere except for the 1hth element of its diagonal we recover the result of Opper & Archambeau (2009). Note the key role the linearity each 8n in this result without Ha, would not necessarily zero everywhere and the contribution of each n term could have arbitrary rank.
A.2 POSTERIOR GAUSSIAN PROCESS
For the sake of notational convenience, we shall assume 1. First, the mean and covariance of q(u) = = mun Kmn) о Po, (u) I Lo (u: Xn, Z) are given by
m, = Kun &k d.
(29)
K uu Kuu PK
where ---- -KK Eu - K The approximate posterior over some latent function value S. obtained by marginalisation of the joint distribution:
a(f+) / po, (fx (f=l))((u)du
/ N (fr: k ky.uKuu ky.s. - ks. kuf. (u; mu, Buu) du
(30)
-N N (Io: ky Kndi ky.f. Ky.nkulku + kf K.lk KinkKuukkkk
Kutkins.)
Substituting in equation results mean and covariance function of the form
nx(x) = ky K K- 1k --1 K k(x, x') = k,p - k K kuy + kfu Pkuf'
(31)
13

---

Under review as a conference paper at ICLR 2021.
B THE GP-VAE
As discuss in Section 3, the GP-VAE is described by the structured approximate posterior
q() Zg () Pol (1) II lo (fni
(32)
where Lo (fniyn) = -I N (fni us (yr), diag o² (rr)) and corresponding ELBO
Pa, (()pp,(y))) (y(r)
= Q())) log Pa, (y(f)
LELBO = =By(S)
log
+ log (0, Φ).
(33)
P,(0,o7Po (f) (fiy)
10, (f:y)
B.1 TRAINING THE GP-VAE
The final term in equation 33 has the closed-form expression
K K
(0
= -I
log
N (Horki O, Keef. + E..k)
(34)
R= A=
log Zoe (0,6)
which can derived by noting that each Zq (0, Φ) corresponds to the convolution between two multivariate Gaussians:
zq (0, N O, K (Mock- -f: O, Deck) df
(35)
Similarly, closed-form expression for Ex(r) (o. (f:y)] exists:
(E))) [log Lo (f:
Byttan) [log bor (frki
N
(Sak Mob (y
I Feg(fak)
log 2 2 (y
20 $1 2 (yn)
2
N
(Aknnn Mon
log 12 o. (yr))
20 (y))
n=1
K
N
Σ 20 2. (y
log M (pu.ss (o) 2
nn
m=1
K log N (fegi Porks Eogi)
N
Σ$
(36)
20 01,
(y
n=1
where Eg = Rr (x,X') and ih mm (X), with
my (x) = kauu (Kuuu + Eer xx(x) = kig., kyyuy (Kn,y, Eg, kunfe-
Mark
(37)
Eg(F) [log Po, [y] is intractable, hence must be approximated by a Monte Carlo estimate. Together with the closed-fome form expressions for the other two terms form unbiased estimate of the ELBO, the gradients of which can be estimated using the reparameterisation trick (Kingma & Welling, 2014).
14

---

Under review as a conference paper at ICLR 2021.
B.2 AN ALTERNATIVE SPARSE APPROXIMATION
An alternative approach to introducing a sparse GP approximation is directly parameterise the struc- tured approximate posterior at inducing points
N
(f) (0,0, po, (f)
II Loi (u; Xn.
Xn, Z))
(38)
1=
where Lor (u; yh Xn: 2), the approximate likelihood, is fully-factorised Gaussian distribution pa- rameterised by a partial inference network:
M
Lo: (u:y nn Xx, Z)
II Jo m=1
II N V (umki (o, (y)))
(39)
In general, each factor Lo (mikk: Y Zrnkk Xn)) conditions on data at locations different to that of the inducing point. The strength of the dependence between these values is determined by the two input locations themselves. To account for this, we introduce the use of an inference network that, for each observation/inducing point pair (umbe: Y yn), maps from (Zmk: X y%) to parameters of the approximate likelihood factor.
Whilst this approach has the same first order computational complexity as that used by the SGP- VAE, having to making forward and backward passes through the inference network K. N M renders significantly more computationally expensive for even moderately sized datasets. Whereas the approach adopted by the SGP-VAL employs deterministic transformation of the outputs of the inference network based on the covariance function, this approach can be interpreted as learning an appropriate dependency between input locations. In practice, found the use of this approach to result in worse predictive performance.
C MEMORY REQUIREMENTS
Assuming input locations and inducing point locations shared across tasks, require storing {K u + Ku ] and K p in memory, which is O (KMN + KM² + +?)).
(t) For the SGP-VAE, we also require storing & and instantiating {u bi,k bo which O (1 + KMD + 2KN).
Collectively, this results in the memory requirement O (KNM + KM² + 2 @i + KMD+ 2K
If were to employ the same sparse structured approximate posterior, but replace the output of the inference network with free-form variational parameters, the memory require- ment is O OKKNM- + KM N + KMD + 2TN Alternatively, if we were to let q(u) to be parameterised by free-form Cholesky factors and means, the memory requirement is O KKNM KM2 N2 + KMD TK ((M + 1)/2 TKM). Table 3 compares the first der approximations. Importantly, the use of amortisation across tasks stops the memory scaling with the number of tasks.
Table 3: A comparison between the memory requirements of approximate posteriors.
((u)
Amortised?
Memory requirement
p(u) II, (u) p(u) II, (u) ((u)
Yes
O (KNM + KM² + NN2 + + O (KNN + KM TTKN)
No
No
O (KNM (KMM²
D MULTI-OUTPUT GAUSSIAN PROCESSES
Through consideration of the interchange of input dependencies and likelihood functions, we can shed light on the relationship between the probabilistic model employed by the SGP- VAE and other multi-output GP models. These relationships are summarised in Figure 4.
Note we only require evaluating single K (0)(()) at each update.
15

---

Under review as a conference paper at ICLR 2021.
SGP-VAE
fP ~ GP(O, x(x,x)) y,lfa ((((f)), E(f..)
DGP
GP-FA
fx ~ CP(O, x(x.x') 5 CP(O, k(I,f))
fa GP(O, x(x,x)) y,Ifn N(WL,, E)
GP likelihood
Linear likelihood
Remove input dependencies
GP-LVM
Factor Analysis
f. No0, I bolf GP(O, k(f,
fn OO y, N(Wf,, E)
VAE
fn (O V, ~ (f((f,) E(f.)
Figure 4: A unifying perspective on multi-output GPs.
Linear Multi-Output Gaussian Processes Replacing the likelihood with a linear likelihood func- tion characterises family of linear multi-output GPs, defined by linear transformation of K inde- pendent latent GPs:
II GP (0, Ko, (x,
k=1
(40)
N
ylʃ II Wy,s Wf,, E).
n=1
The family includes Teh et al.'s (2005) semiparametric latent factor model, Yu et al.'s (2009) GP factor analysis (GP-FA) and Bonilla et al.'s (2008) class of multi-task GPs. Notably, removing input dependencies by choosing Ao,,.. (x,x') = 6(x, x') recovers factor analysis, or equivalently, probabilistic principal component analysis (Tipping & Bishop, 1999) when Σ = =_,, Akin to the relationship between factor analysis and linear multi-output GPs, the probabilistic model employed by standard VAEs can be viewed as special, instantaneous case of the SGP- VAE's.
Deep Gaussian Processes Single hidden layer deep GPs (DGPs) (Damianou & Lawrence, 2013) are characterised by the use of a GP likelihood function, giving rise to the probabilistic model
Il gP (0, ko, (x,x'))
km
(41)
P
vff I ko (0) (f(x)f(x)))
where y, >xn). The GP latent variable model (GP-LVM) (Lawrence & Moore, 2007) is the special, instantaneous case single layered DGPs. Multi-layered DGPs are recovered using hierarchical latent space with conditional GP priors between each layer.
E EXPERIMENTAL DETAILS
Whilst the theory outlined in Section 2 describes general decoder parameterising both the mean and variance the likelihood, we experienced difficulty training SGP-VAES using learnt variance, especially for nigh-dimensional observations. Thus, for the experiments detailed in this paper use shared variance across all observations. We use the Adam optimiser (Kingma & Ba, 2014) with constant learning rate of 0.001. Unless stated otherwise, estimate the gradients of the
16

---

Under review as a conference paper at ICLR 2021.
ELBO using single sample and the ELBO itself using 100 samples. The predictive distributions are approximated as Gaussian with means and variances estimated by propagating samples from q(f) through the decoder. For each experiment, normalise the observations using the means and standard deviations of the data in the training set.
The computational complexity of performing variational inference (VT) in the full GP-VAE. per update, is dominated by the O NK)* cost associated with inverting the set of K N x N matrices, (Krate Eo, I This can quickly become burdensome for even moderately sized datasets. A pragmatic workaround is to use biased estimate of the ELBO using N < N data points:
N N
Po, (f/f)
CN ELBO
E BO0
log
+ log (0,0,
(42)
la yf
y and denote the mini-batch of Ñ observations and their corresponding latent variables, re- spectively. The bias introduced due to the normalisation constant, which does satisfy N log [ (). = E [log (0, o)]. Nevertheless, the mini-batch estimator will be a reason- able approximation to the full estimator provided the lengthscale of the GP prior is not too large." Mini-batching cannot be used to reduce the O (KN3) cost of performing inference at test time, hence sparse approximations necessary for large datasets.
E. SMALL-SCALE EEG
For all GP-VAE models, we use a three-dimensione latent space, each using squared exponential (SB) kernels with lengthscales and scales initialised to 0.1 and 1, respectively. All DNNs, except for those in PointNet and IndexNet, use two hidden layers of 20 units and ReI.U activation functions. PointNet and IndexNet employ DNNs with a single hidden layer of 20 units and a 20-dimensional intermediate representation. Each model is trained for 3000 epochs using batch size of 100, with the procedure repeated 15 times. Following (Requeima et al., 2019), the performance of each model is evaluated using the standardised mean squared error (SMSE) and negative oielliilldood (NLL). The mean + standard deviation of the performance metrics for the 10 iterations with the highest ELBO is reported."
E.2 JURA
We use two-dimensional latent space for all GP-VAE models with SE kernels with lengthscales and scales initialised to 1. This permits a fair comparison with other multi-output GP methods which also use two latent dimensions with SE kernels. For all DNNs except for those in IndexNet, two hidden layers of 20 units and ReLU activation functions. IndexNet uses DNNs with a single hidden layer of 20 units and 20-dimensional intermediate representation. Following Goovaerts (1997) and Lawrence (2004), the performance of each model is evaluated using the mean absolute error (MAE) averaged across 10 different initialisations. The 10 different initialisations are identified from a body of 15 as those with the highest training set ELBO. For each initialisation the GP- VAE models are trained for 3000 epochs using a batch size of 100.
E.3 LARGE- EEG
In both experiments, for each trial in the test we simulate simultaneous electrode "blackouts' by removing any sample period at random with 25% probability. Additionally, we simulate individual electrode "blackouts' by removing any 16 sample period from random with 50% probability from the training set. For the first experiment, we also remove any 16 sample period at random with 50% probability from the test set. For the second experiment, we remove any sample period at random with 10% probability. All models are trained for 100 epochs, with the procedure repeated five times, and use 10-dimensional latent space with SE kernels and lengthscales initialised
'In which case the off-diagonal terms in the covariance matrix will be large making the approximation po; (f) = ro, (f) extremely cnide.
We found that the VAE oecasionally got stuck in very poor local optima. Since the BL is calculated
the training set alone, the experimental procedure is still valid.
17

---

Under review as a conference paper at ICLR 2021.
and 0.1, respectively. All DNNs, except for those in PointNet and IndexNet, use four hidden layers of 50 units and ReLU activation functions. PointNet and IndexNet employ DNNs with two hidden layers of 50 units and 50-dimensional intermediate representation.
E.4 BOUNCING BALL
To ensure a fair comparison with the SVAE and SIN, we adopt an identical architecture for the inference network and decoder the original experiment. particular, we use DNNs with two hidden layers of 50 units and hyperbolic tangent activation functions. Whilst both Johnson et al. and Lin et al. use eight-dimensional latent spaces, we consider a GP-VAE with a one-dimensional latent space and periodic GP kernel. For the more complex experiment, we use SGP- VAE with fixed inducing points placed every 50 samples. We also increase the number hidden units in each layer of the DNNs to 256 and two-dimensional latent space one for each ball.
E.5 WEATHER STATION
The spatial location of each weather station is determined by its latitude, longitude and elevation above sea level. The rates of missingness in the dataset vary, with 6.3%, 14.0%, 18.9%, 47.3% and 93.2% of values missing for each of the five weather variables, respectively. Alongside the average temperature for middle five days, we simulate additional missingness from the test datasets by removing 25% of the minimum and maximum temperature values. Each model is trained on the data from 1980 using a single group per update for 50 epochs, with the performance evaluated on the data from both 1980 and 1981 using the root mean squared error (RMSE) and NLL averaged across five runs. We use three-dimensional latent space with SE kernels and lengthscales initialised to 1. All DNNs, except for those in PointNet and IndexNet, ReLU activation functions. PointNet and IndexNet employ DNNs with two hidden layers of 20 units and 20-dimensional intermediate representation. Inducing point locations are initialised using k- means clustering, and are shared across latent dimensions and groups The VAE uses FactorNet. We consider independent GPs modelling the seven point time series for each variable and each station, with model parameters shared across groups. No comparison to other sparse GP approaches made and there no existing framework for performing approximate inference in sparse GP models conditioned on previously unobserved data.
four hidden layers of 20 units and
F FURTHER EXPERIMENTATION
F.1 SMALL SCALE EXPERIMENTS
Table 4: A comparison between multi-output GP models on the EEG and Jura experiments.
Metric
GP-VAE
GPPVAE
SMSE
0.24 (0.02)
0.524 (0.11)
EEG
NLL
2.01 (0.28)
1.92 (0.05)
MAE
0.40 (0.01)
0.49 (0.02)
Jura
NLL
1.00 (0.06)
0.99 (0.04)
Table 4 compares the performance of the GP- VAE to that of the GPPVAE, In all cases, FactorNet is used to handle missing data. We emphasise that the GP-VAE and GPPVAE employ identical probabilistic models, with the only difference being the form of the approximate posterior. The perior predictive performance of the GP- VAE can therefore be accredited to the use of the structured approximate posterior as opposed to the mean-fiel approximate posterior used by the GPPVAE.
F.2 SYNTHETIC BOUNCING BALL EXPERIMENT
The original dataset consists of 80 12-dimensional image sequences each of length 50, with the task being to predict the trajectory of the ball given a prefix longer sequence. The image sequences generated at random by uniformly sampling the starting position of the ball whilst keeping the
18

---

Under review as a conference paper at ICLR 2021.
E
Postector latent GP
100
150
200
Time
Figure A comparison between the mean of the GP-VsEEs posterior predictive distribution (mid- dle) and the ground truth (top) conditioned noisy observations up to the red line. The latent approximate GP posterior is also shown (bottom).
bouncing frequency fixed. Figure 5 compares the posterior latent GP and of the posterior predictive distribution with the ground truth for single image sequence using just a single latent dimension. As demonstrated in the more more complex experiment, the GP-VAE able recover the ground truth with almost exact precision.
Following Lin al. (2018), Figure Ia evaluates the r-steps ahead predictive performance of the GP-VAE using the mean absolute error, defined as
1 (T r)d IIy r
( [n,,4rlll
(43)
where Niett is the number of test image sequences with T time steps and y denotes the noiseless observation time step + T.
G PARTIAL INFERENCE NETWORK COMPUTATIONAL GRAPHS
19

---

Under review as a conference paper at ICLR 2021.
(1,31)
shared
Σ
(P.yp)
rp
(a) PointNet
11
h1
:
P
n
yP
hr
Fp
(b) IndexNet
31
hn
m
I
"
hp
1P
(c) FactorNet
Figure 6: An illustration of the three different partial inference network specifications discuss in Section 2.4. r7 denotes the vector of natural parameters of the multi-variate Gaussian being parame- terised.
20


Reviews:
[{"review_id": "wvQC_xEi0z8", "reviewer": "AnonReviewer2", "summary": "", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "This paper attempts to enhance the inference in multi-output GPs on previously unobserved data using amortised variational inference. This model is proposed to overcome the shortcomings in current GP-DGMs. The major concerns regarding this paper are summarized as below.\n\nOne of the two main contributions in this paper is introducing sparse GP via amortized variational inference in order to replace the time-consuming GP, the novelty of which however seems to be incremental.\n\nThe modeling of task dependency should be highlighted. In my understanding, this paper employs a MLP to mix K latent functions in a nonlinear manner through equation 1 for each output. The NN parameters are shared across tasks, which however might be a too strong assumption and thus could deteriorate the performance of multi-task learning.\n\nThe comparative results in Table 1 cannot showcase the superiority of GP-VAE in comparison to the state-of-the-art multi-output GP. It performs worse than GPAR in terms of NLL on the EEG dataset, and is comparable to GPAR on the two datasets in terms of SMSE.\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "mI3riFNkufu", "reviewer": "AnonReviewer4", "summary": "", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "The model proposed by the paper is a multi-task, multi-output GP, with likelihoods driven by a decoder network. While this model is not new --- same as that of Pearce 2020 as acknowledged by the authors the related work section --- this paper includes the sparse approximation and also handles missing data.\n\nThe paper can be made clearer:\n1. The paper introduces \"partial inference network\" in section 2.1. Since this plays a central part in the contribution of the paper, I hope the authors can expand on this. For example, in what way is it \"partial\" and how is it different from normal inference networks, and what are the challenges. From a GP and Bayesian inference perspective, missing data is no problem, as acknowledge by the authors in section 2.3. Hence, why is it that the normal mechanism fails here? Although section 2.3 is section dedicated to this, it is insufficient. \n\n2. N and P in the last para of section 2.1 are not defined.\n\n3. The first and third shortcomings of sparse GP frameworks claimed by the authors in the last three sentences of section 2.1 do not really exist in existing sparse approximations to GP. Some clarification is needed here.\n\n4. Some figures accompanying the different approaches in section 2.3 will be very helpful.\n\n5. The section title for section 2 is \"spatial-temporal\". It is good to lay out what is the spatial and what is the temporal aspect for each of the data sets in section 4. For example, I do not really find the \"temporal\" dimension in Jura. One way around this is to have a more general section title.\n\nThe paper is contributes to the community and the work is correct, so *accept*. It is *not a stronger accept* because I find the \"sparse\" contribution minimal given the vast literature of sparse approximations to GP; and I find the \"missing-data\" contribution confusing.\n\n\n\n\n\n\n", "extra_scores": {}, "extra_text": {}}, {"review_id": "_8gd1mt-sTr", "reviewer": "AnonReviewer1", "summary": "", "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": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "In this work generative models using a GP as prior and a deep network as likelihood (GP-DGMs) are considered. In the VAE formalism for inference, the novelty of this paper is located in the encoder: It is sparse and the posterior can be computed even when part of the observations are missing. Sparsity is obtained using inducing inputs and the missing observations are handled through the use of deep sets, i.e. the observations aren't given as a vector, but as a permutation-invariant set of (index, value) pairs.\n\nThe idea is  interesting and well executed, and the experimental results are certainly good when compared with the provided baselines. However, there are many unclear details in the experiments. I was left with quite a few unanswered questions:\n- In the cases in which no imputation is needed, how does the proposed SGP-VAE compare with a standard encoder (instead of one based on deep sets). I.e., are we paying a price for the ability to handle missing data, or are we getting even better results than with a traditional encoder?\n- How does the model compare with prior similar work, such as the GPPVAE? In principle, it should be easy to improve on this baseline.\n- How good is the model itself? Here we see results based on amortized inference, but what if we weren't using it? Even though this case would be slow and maybe not able to handle the entire dataset, this can be done in a subset. In general, it is unclear how much is being lost by a) amortized inference; b) sparsity; c) use of deep sets (in this latter case, it could that something is won instead).\n- In the Japanese Weather experiment, when the standard VAE is used (with no mean imputation), how are the unobserved values being dealt with?\n- In the Japanese Weather experiment, the baselines don't seem very strong. A standard VAE and independent GPs don't seem to be particularly suited to this spatiotemporal weather prediction.\n- In general, I'm not sure where the advantage seen in these experiments is coming from. For instance, when compared with published results on the same data (e.g., GPAR on EEG), the SGP-VAE doesn't seem to offer a big advantage. In the other experiments, no strong baselines for that specific  data seem to be provided.\n- I would be curious to know in which of these datasets the model itself is superior (if we were to skip the encoder and use, say, MCMC for inference) and in which ones it isn't.\n\nThanks for your response, clarifying.\n", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: The paper proposes a method for inference in models with GP priors and neural network likelihoods for multi-output modelling, dealing with the problem of scalability and missing data. The paper builds upon previous work on inducing variables for scalability on GP models and inference networks for amortization (reducing the number of parameters to estimate) and dealing with missing data.

There are several concerns about the paper in terms of generality/flexibility of the approach, as the proposed model shares the NN parameters across tasks and the results on the small datasets do not show improvements wrt baseline such as GPAR. The authors’ comments provide somewhat satisfactory replies to these issues. Nonetheless, the major drawback of this paper is its novelty as the ideas on the paper have been explored extensively in the GP literature. Although the authors do make a case for scalability when using inference networks, there are other previous works that perhaps the authors are unaware of, for example, https://arxiv.org/abs/1905.10969   and even more sophisticated inference algorithms than can serve as truly state-of-the-art competing approaches (for example based on stochastic gradient Hamiltonian Monte Carlo, https://arxiv.org/abs/1806.05490). 
Decision: Reject