Revisiting Active Sets for Gaussian Process Decoders
Authors: Pablo Moreno-Muñoz
Cilie W. Feldager
Søren Hauberg
Venue: NeurIPS 2022 Accept (2022)
Keywords: Gaussian Processes

Abstract: Decoders built on Gaussian processes (GPs) are enticing due to the marginalisation over the non-linear function space. Such models (also known as GP-LVMs) are often expensive and notoriously difficult to train in practice, but can be scaled using variational inference and inducing points. In this paper, we revisit active set approximations. We develop a new stochastic estimate of the log-marginal likelihood based on recently discovered links to cross-validation, and we propose a computationally efficient approximation thereof. We demonstrate that the resulting stochastic active sets (SAS) approximation significantly improves the robustness of GP decoder training, while reducing computational cost. The SAS-GP obtains more structure in the latent space, scales to many datapoints, and learns better representations than variational autoencoders, which is rarely the case for GP decoders.

Revisiting Active Sets for Gaussian Process Decoders
Pablo Moreno-Muñoz"
Cilie W. Feldager*
Søren Hauberg
Section for Cognitive Systems Technical University of Denmark (DTU) {pabmo, cife, sohau)adtuuuuuuauu dk
Abstract
Decoders built on Gaussian processes (GPs) enticing due to the marginalisation over the non-linear function space. Such models (also known as GP-LVMs) are often expensive and notoriously difficult to train in practice, but can be scaled using variational inference and inducing points. In this paper, we revisit active set approxi- mations. We develop new stochastic estimate of the log- -agrgial likelihood based on recently discovered links to cross- volss-vali,,ionn, and we propose a computationally efficient approximation thereof. We demonstrate that the resulting stochastic active sets (SAS) approximation significantly improves the robustness of GP decoder training, while reducing computational cost. The SAS-GP obtains more structure in the latent space, scales to many datapoints, and learns better representations than variational autoencoders, which is rarely the case for GP decoders.
Introduction
Generative models can be viewed as regression models from unknown inputs. That is. we assume = ==),,, where f' is an unknown mapping from latent variables z to observations æ. Given the inherent difficulty of this task, it is perhaps sensible to marginalize the unknown mapping f to avoid the brittleness of point estimates. This is the driving idea in the Gaussian process latent variable (GP-LVM) (Lawrence, 2005), which places Gaussian process (GP) prior over the unknown mapping and marginalizes accordingly. This contrasts contemporary generative models that predominantly operate with point-estimates of f (Kingma and Welling, 2013; Sohl-Dickstein et al., 2015). While conceptually elegant, the GP-LVM is, however, notoriously difficult to train and the conceptual benefits often not realized in practice.
Exact inference involves computing the marginal likelihood, but (like other GP methods) its cubic complexity in the number of observations O(N3) limits the scalability of the GP-LVM. However, the idea of marginalizing the decoder is sufficiently attractive to motivate the development of scalable and reliable training techniques: Its Bayesian formulation (Titsias and Lawrence, 2010) variationally integrates out the latent variables z to obtain an evidence lower bound. Using auxiliary inducing variables, Snelson and Ghahramani (2006) expanded GP regression, which is also applicable in the unsupervised learning setting (i.e. the GP-LVM).
However, considering inducing variables here involves dangers. First, the convergence of inducing points is well-studied the supervised GP scenario, where inputs are fixed, but it differs from the unsupervised case where the inputs are estimated. Bauer et al. (2016) notes that even in the supervised setting, inducing points are "not completely trivial to optimise, and often tricks [...] are required". and we hypothesize that this is further complicated in the unsupervised setting where we optimize both latent and inducing variables while they interact.
In this paper, we revisit active sets for scaling GP decoders, a sparse approximation predominantly used before the seminal work of Snelson and Ghahramani (2006). From practical viewpoint, active
Equal contribution.
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

sets are fixed inducing variables that belong to the training dataset. We make links between such active sets and voos-validation, allowing us lean on recent result from Fong and Holmes (2020) which, in turn, links cross validation and the log-marginal likelihood. We show how these links allow for a stochastic estimate of the log-marginal likelihood, and that active sets can be seen as a computationally efficient approximation of this. Practically, this amounts to repeatedly and randomly sampling active sets rather than trying to find the optimal active set. We denote this framework as stochastic active sets (SAS). We demonstrate that SAS consistently results in significantly better-fitted GP decoders over models trained using inducing points.
Historical remarks. Methods based on subsets of data diminish the computational demand and were first introduced in a GP context in the foundational work on sparse approximations by Smola and Bartlett (2001). Back then, QuifooerroCCandel and Rasmussen (2005) had already pointed out the main difficulties behind the optimal selection of subsets:
"Traditionally, sparse models have very often been buill upon carefully chosen subset of the training inputs. [...] In sparse Gaussian processes it has also been sug- gested to select the inducing inputs from among the training inputs. Since this involves a prohibitive combinatorial optimization, greedy optimization approaches have been suggested [...]. Recently, Snelson and Ghahramani (2006) have proposed to relax the constraint that the inducing variables must be subset of training/test cases, turning the discrete selection problem into one of continuous optimization."
This explains how inducing variables reshaped the Gaussian process community, effectively banishing other subset-based methods. Our work builds the advances made in stochastic optimization in the time since active sets were left behind. We show a third way over those considered by Quiñonero- Candela and Rasmussen (2005): instead of optimizing the active set, average with respect to it. This simplifies matters notably and makes them more robust.
To justify approach, we establish a link between active sets and cross validation (CV). The latter has a long history for model selection in GPs, dating at least to the seminal work of Wahba (1990). For probabilistic models, Rasmussen and Williams (2006) point to the utility of CV variants within negative log- probabilities.. Building on results from Fong and Holmes (2020) linking CV and log- marginal likelihoods, we argue that, for GP-LVMs, active sets combine more gracefully with stochastic optimization. The remainder of this paper elaborates on this viewpoint and demonstrates empirically.
2 Gaussian Process Decoders
The Gaussian process latent variable model (GP-LVM) (Lawrence, 2005) defines a decoder which is non-linear mapping² x f(z) from the latent space Z E Rọ to observation space X' E R.. The prior on this map is a Gaussian process (GP) so it is drawn like f ~ CP ko(,,)) where ke is the covariance function kernel and KN denotes the evaluated kernel function so the i, jthh element of KNN equals ko to(2;, (z,, 2z). For clarity, we omit the dependence on covariance function parameters, O.
The original version of the GP-LVM starts from one-dimensional observations x {<r,, JN and latent variables z = } and factorises the joint distribution of the model as p(æ, flz) = p())),,))p(fl Here the conditional distributions correspond to the likelihood model and the prior
N
p(zlf,2)
II ([,
p(flz) = N(f=))0 KNN).
(1)
When the data dimensionality is D> the model factorises across dimensions, and we have different mappings ff per dt feature. One of the principal assumptions of the GP-LVM is that the prior p(f)) regularises the smoothness of all mappings equally, so we only consider one kernel. This assumption can be relaxed if needed, but at increased computational cost and with more learnable hyperparameters.
We use (æ, } to denote observations and latent variables respectively, since we do not consider inducing variables. Notice that Lawrence (2005) use the notation (y. x).

---

Mapping marginalization. A GP prior over the non- non-lineat linear decoder f alloos for marginalisation of obtain a closed-form expression of the marginal likelihood of the GP-LVM
p( 2)
pla=|S, =))([)))d)))))) = =a NNN + d²T).
On a log-scale, this gives the following objective function (Lawrence, 2004), which can optimized w.r.t. both hyperparameters 0 and latent variables
DN log 2T
D
log KNN + + - u + + xx
(2) The difficulties of training the GP-LVM using this objective function are evident above, as the evaluation cost grows cubically with the number of observations N. Furthermore, notice that observations æ are no longer independent (in contrast with Eq. 1) once f is integrated out. The log-marginal likelihood will be the starting point for our approach in Sec. 3.
C = 2
2
Bayesian extension. In the seminal works of Lawrence (2004, 2005), the GP-LVM is first derived as a non-linear extension of probabilistic principal component amalysis (PPCA) (Tipping and Bishop, 1999). Considering the isotropic prior on the latent variables z, such that p(=,) (=n, [0, 1) Vz, E z. the general idea is to optimize them rather than introducing marginalization. From a full probabilistic perspective, one could also interested in the posterior distribution over 2, which leads to using Bayesian inference for the GP-LVM approach. This is the driving idea of Titsias and Lawrence (2010), where variational methods are introduced. In particular, latent variables are not easy to marginalize, mainly due to their presence in the kernel mappings, so lower-bound on the log-marginal likelihood log p(x) = log p(æ) 2z 1))((2 (=)ddz of the model is derived.
So far, the Bayesian GP-LVM model (Titsias and Lawrence, 2010) has been considered as the standard methodology to apply GPs to large unsupervised datasets with 100-10" observations, e.g. in regression (Bui and Turner, 2015), classification (Gal et al., 2015) and representation learning (Märtens et al., 2019) tasks.
3 Stochastic Active Sets
Our key objective is a computationally efficient estimate of the log- marginal likelihood in Eq. 2, as this is known be good measure of generalization performance (Rasmussen and Ghahramani, 2000; Germain et al., 2016). Another popular measure of generalization performance is cross validation (CV) (Geisser and Eddy, 1979; Vehtari and Lampinen, 2002), which is arguably mostly used outside the realm of Bayesian models. Recently, Fong and Holmes (2020) linked these two measures, effeetiiely showing that the marginal likelihood is equivalent to the average over exhaustive leave-Rt-out CV scores. In particular, the average is w.r.t. the size of the hold-out set. More precisely, let
C 1 Scv(a|R) = Σ
Σ
Log (2,, E. EER,
log p(z, Ap 2)
(3)
p=1 BER,
denote the leave-R-oud CV using prg-iiicietive scoring functions Dog A, 2). Here Ap denotes the active set indices of the training data, such that Ap C (1, N) and R, (1,2..... NJ Ap are the remaining hold-out samples. The subscript E denotes the permutation, and we average over all
(N) possible hold-out sets. We use use R to indicate the size of the hold-out set RR (R = Rp and let = |Ap = N R. we average Sev(x)R) over all possible sizes of the hold-out set, then Fong and Holmes (2020) has shown that recover the exact -oramnagi likelihood,
R
N
logp(æ//z) Σ Scv( Secv(x))B) + Spev(x R).
(aalogp(æ=lzz... and Sccv R) = Σ Scv Scv (x|r) the cumulative CV score, which reduces to sum of expectations over the predictive factors. 3 Further details Spcv and Secv can be found in the Appendix. Fong and Holmes (2020) used Eq. 4 to argue in favor of using the marginal likelihood over cross-validation for model selection.
Here, (pcv(
We drop the permutation subscript p in A and R to avoid cluttered notation.
3

---

D2
D2 EAA2
D2 E
21
31
CiC.Ag
24 E
B4 E-A2
2g E A
2i EA
A1
23 E.A2
33
Figure I: Schematic graphical model of the correlation structure given different permutations of the active set A2....... e for five observations (.),.... 2t Thick red lines indicate that we build full covariance densities between observations included in 2.,4- Only-blue variables are considered conditionally independent among them w.r.t. to the red colored ones. Dashed lines indicate the conditional probability factors
z) that we can easily compute.
Stochastic approximation. We take a slightly different view than Fong and Holmes (2020) and argue that Eq. can be the grounds for an efficient stochastic gradient (Robbins and Monro, 1951) of the log- -arrginal likelihood suitable for training. In the context of GPs, we have that conditional probabilities for Spcv (ar) R) and Secv (ar RR) Eq. are
p( (A) = Al0 KAa + II)
p(zn| =) = mm CalA);
(5)
where we used Eq. 2.22 from Rasmussen and Williams (2006) to obtain
mnIA Knd(KAn + -
CalA = Knn + to KnA KAA + TT) K nA:
and K KAA AA RRAAA has entries k(zi, zz) with Zi, Zy E ΣĄ. The computational cost of Spcv R) is ((4)), while log p(æn #) can reuse the matrix inversion K AAA to only have additional linear cost. Clearly, we can obtain an unbiased stochastic estimate of the log-margina likelihood by first
uniformly sampling R, second making a random split permutation into train and hold-out data and finally by evaluating
log p(xp)))
Σ nER
Logp(zz,/au x) + log p(z))
(6)
where we remark that the summation can be mini-batched. This approach is equivalent to the decomposition log p(z) 2)) =)= = p(o x) ( + log p(æAalzA) and also assumes conditional inde- pendence among observations 2, for n N.R. This is similar to the standard active set approximation (Seeger et al., 2003) previously discussed, and we may think of A as stochastic active set (SAS).
However, the estimate of log p(z]2) still has the same computational complexity O(N3) as the usual deterministic approach, since may sample A= N in the worst case. Instead, we propose to make a stochastic approximation where we choose the size of the active set deterministically through a user-specified parameter, such that the computational cost reduces to ((43), in sparse approximations based on inducing points. This does not need to be unbiased; empirically we find that most often the approximation behaves as lower bound the true marginal likelihood, and that, in all instances, it is rather close approximation. We include a longer discussion on this point in the Appendix and the training methodology using SAS is in Alg, 1.
3.1 Extension for Bayesian GP decoders
We next seek to extend the previous SAS approach to the Bayesian GP-LVM, where we aim to obtain the posterior distribution p(z|2z). From this perspective, we are interested in marginalising latent variables z to obtain the marginal likelihood p(æ) of the model. However, this integration is not possible, as latent variables appear non- linearly the kernel function (Titsias and Lawrence, 2010). Alternatively, we consider the variational inference scheme, where an auxiliary distribution q(z) is
*Notice that probabilistic objective function changes between standard and Bayesian versions of the GP-LVM. In the former case, we usually look for p(z) 22) as the marginal likelihood w.r.t. the mapping f. This the one usually considered in supervised GP tasks. In the latter, we refer to p(œ) as the marginal likelihood of the model, since also integrated out.

---

introduced into the formulation. Therefore, we are able to build the evidence lower bound (ELBO) of the model using Jensen's inequality as
log)(x) 2
g(2) p())pp((()))d
By(=) -
which equivalent to the one obtained by Titsias and Lawrence (2010, Eq. 8). At this point, the computational cost is ((N²), since the ELBO requires evaluating log p(z] where we invert KN The expectation in Eq. can also benefit from stochastic SAS approximation, just as with inducing points (Hensman et al., 2013). Thus, the lower bound can be approximated as
LELBO
Σ Bete.) 1og (g((a) [log((0A|zz)) nER
Σ KLLg(ze))ll(za)) (8)
where we consider mean-field VI to factorize the distribution g(z). The Bayesian GP-LVM shares high-level similarities with other generative models that marginalize the latent variable according to simple prior (Rezende and Mohamed, 2015). The proposed SAS approximation (8) scales similarly to mini-batched inducing point approximations, but will later see that SAS behaves notably better in practice. Algorithmically, the approach is simple, and the summary code is provided in Alg. 2.
Algorithm 1 SAS for GP decoders
Algorithm 2 SAS for Bayesian GP decoders
1: Input: Observed data æ
1: Input: Observed data
2: Parameters: Initialize 0, db
// 110, NA
2: Parameters: Initialize 0, Φ // 1,, NA
3: for in epochs do
for e in epochs do
4: for in batches do
4: for b in batches do
5: Sample Pbatch ~
Sample Fbacch x
6: FR. A < random_split(annn
6: DR; AAA +- random_nnspppp (Xaach))
7: if amortized then
if amortized then
8:
Get (2R, 2A] +- TA
8: Get μ/z 9p (R Ab
9:
end if
Get 0= 4 9o (TR, (A
10:
Compute K
1/ via Cholesky
10:
end if
11:
Evaluate log pl A (.AAlz... (
11:
Sample {#R, ~ alkxz 0z) /1 RT
12:
Evaluate log pl R. A 2), Van E
12:
Compute AAA
/1 via Cholesky
13:
Evaluate Eq. 6
13: Evaluate L in Eq. 8
14:
do Adam(O. 0 Φ) step for L
14: do Adam(0, +) step for LFLBO
15:
end for
15:
end for
16:
end for
16:
end for
NA: Non-amortized.
NA: Non-amortized, RT: Reparametrization trick.
3.2
The Role of Amortization
Early after the emergence of the seminal GP-LVM (Lawrence, 2005), the lengthy optimization required obtaining a result in which all latent representations 2, became noticeable. An additional consideration is that, while most approaches non-linear wwwddimensoonaltty reduction focus on pre- serving similarities, the GPLVM does the opposite. This property was initially discussed by Lawrence and Quiñonero- Candella (2006), since some sense the GP-LVM is dissimilarity preserving, such that different observations will generally be represented far away from each other. In practice, we often more interested in embeddings that reflect the true distance between the observed objects in their representations, particularly those that are close together. This observation inspired buck constraints for locality preservation (Lawrence and Quiñonero-Camma Candela, 2006), which enforces latent variables g(x((b parameterized by Φ. This is similar to the
to be an explicit function of observations stochastic encoder applied variational autoencoders (Kingma and Welling, 2013; Rezende and Mohamed, 2015). This idea also extended to the VI framework in GP-LVMs by Bui and Turner (2015) using a recognition model, e.g. ((z,) = (2 (xnl()), go (z,, Φ)) and more recently, to accelerate hyperparameter learning in GPs with hierarchical attention networks (Liu et al., 2020).
In our context, we assume the mappings to be neural networks (NNs) like Bui and Turner (2015), referring to them as amortization networks. We find such networks accelerate learning very nicely when used in conjunction with SAS. It is also worth noting that amortization has been empirically shown to improve generalization performance (Shu al., 2018).
5

---

4 Related Work
Marginal likelihood approximations were used in GPs (Smola and Bartlett, 2001; Csató and Opper. 2002) before the apparition of pseudo-inputs (Snelson and Ghahramani, 2006) and the associated variational inference framework (Titsias, 2009). In the former case, stochastic approximations ELBO were first presented by Hensman et al. (2013, 2015), in line with the Bayesian counterpart of SAS. In terms of active sets, Seeger et al. (2003) empirically observed that the approximation was generally stable for optimisation, even if the size of A was small fraction of the training size only. However, they also observed that random selection of active sets led non-smooth fluctuations, making it difficult to converge through exact gradient ascent. Particularly, the issue with re- re-selecting of A motivated Snelson and Ghahramani (2006) as a smoother optimization alternative, and we find that SAS also circumvents this problem via stochastic estimates shown in Sec. 5.
the
The connection between cross-validation and GPs was already described in Rasmussen and Williams (2006) as an alternative for model selection. However, the equivalence between exhaustive CV and the log- -ng-magginal likelihood provided in Fong and Holmes (2020) provides novel perspective that we exploit. Additionally, the notion of back constraints has recently been considered in Lalchand al. (2022) for GP-LVMs with inducing points, where doubly stochastic formulation is used. More recently, amortization networks have been used to drastically reduce the number of inducing points in supervised GPs.
Experiments
In this section, we evaluate the performance of the SAS approach for stochastic learning of GP decoders, both the deterministic GP-LVM (Lawrence, 2004) and Bayesian counterpart (Titsias and Lawrence, 2010). For this purpose, we consider three different datasets: MNIST (LeCun et al., 1998), FASHION MNIST (Xiao et al., 2017) and CIFAR-10 (Krizhevsky, 2009). For a fair comparison of our model with baseline methods, use the same amortization across models, namely neural network (three linear layers ReLU activation functions) for all models in all experiments and all GPs use an quadratic exponentiated kernel. In all experiments, the learning rates are set in the range [10- 10-0) 21. the maximum number of epochs considered is 300 and we use the ADAM optimizer (Kingma and Ba, 2015). For SAS experiments, we only consider batch sizes greater than the active set size, as this is a requirement for SAS.
Performance metrics of the SAS-GP decoders are given in terms of the negative log-predictive density (NLPD), root mean-square error (RMSF) and mean absolute error (MAE). all cases, we optimize w.r.t. approximation to the log- log-marginal likelihood log p(x|2) in the deterministic scenario or w.r.t. a lower bound on the logp(x)) of the model in the Bayesian cases. We also provide PYTORCH code that allows for reproducing all experiments and models." We monitor the run-time of convergence as we suspect that rotating active sets (see Fig. I) across the dataset is a fast way to capture the correlation structure of the regression problem.
5.1 Representation performance
First, we analysed our SAS approach and the approximate optimization procedure on an unsupervised version of MNIST, FMNIST and CIFAR- 10, where we took the full training corpus for learning two- dimensional latent representations of images. The approximation curves are included in Fig. 2, where we observe convergence in less that 2h runtime in CPU for most cases. For all experiments, we can observe that for larger active set sizes A, the SAS approaches take longer time to complete the 300 epochs, as the computational cost of inversion is higher.
Evaluation with SOTA methods. We also tested the performance of representation learning in the state-of sfate-of-the-eee oh-ta--rttttttto methods with and without GPs for unsupervised scenarios. Results are shown in Fig. 4. We include a short description of the models considered:
SAS-GP DECODER - uses stochastic active sets to approximate the log-marging likeli- hood log =)) The training methodology is described in Alg. 1.
*Please see the supplementary material for more details. 'The code is publicly available in the repository: https://githut commppoorreoo/SAASPP including baselines.

---

SAS A MINIST
X30
SAN VS A. PANIST
SAS VR. A CIBAR
-0.5
-0.6
-0.5
1.0
00.5
-122
-1.0
0.5
1.0
Fim
Run Time [ouus
Rum Time [horas]
X107
Bayes MNIST
x107
Boges
PMNIST
X300
Bages A CHPAR
0.3
00.5
0.00
1.3
1.5
2.35
0.3
0.5
1.0
1.5
0.0
0.5
3.0
1.5
Run Time [houn]
Run Time (horza]
Rim Tioe (beu))]
Figure 2: Approximate log-marginal likelihood (upper row) for SAS and ELBO curves (lower row) for Bayesian SAS. We fix the batch size in Alg. to different active set sizes A. All values in the curves displayed are per-epoch.
1024 and study the convergence for
SAS Latent representation
Bayesian SAS
MNIST
PASHION MNIST
FASHION MNIST
1.2
-10
Figure 3: Latent representation a 2-dimensional space for the 10 MNIST and FMNIST classes learnt with SAS and Bayesian SAS. Notice that the likelihood model of the GP decoder is controlled by a vanilla RBF kernel. The examples have been obtained using an active set size A= 800.
SAS-GP
Bayesian SAS-GP
GP-LVM
Bayesian GP-LVM
VAE
MNIST
FashionMNIST
000
Figure 4: Illustration of latent space mappings Zn E on test data for different SAS models and baselines (rows) and different datasets (columns). The models have been trained on full MNIST and FMNIST.

---

BAYESIAN SAS-GP DECODER It uses stochastic active sets to approximate the ELBO in Eq. (8) on log p(z). See Alg.
BAYESIAN GP-LVM - is based on the model Titsias and Lawrence (2010). Means and variance parameters are generated by amortization NN in Bui and Turner (2015). The model is trained using stochastic variational inference (Hensman et al., 2013).
GP-LVM - We used the model proposed by Lawrence (2005) enhanced with an NN encoding to the latent space. The model is trained using maximum likelihood (ML).
VARIATIONAL AUTOENCODER (VAB) (Kingma and Welling, 2013) The NN encoder has the same architecture as the amortization function used in the SAS-GP models.
For the SAS-GP decoder and GP-LVM, a neural network encodes the latent variables. Likewise, for the Bayesian SAS-GT decoder and the Bayesian GP-LVM, two neural networks each encode the latent means and latent variances. We refer to this encoding amortization. All models use Gaussian likelihood'.
Table I: Comparative metrics for SAS and Bayesian SAS on MNIST, FMNIST and CIFAR-10.
MODEL
SAS
RAYESIAN SAS
ACTIVE SET SIZE
A = 100
A = 200
A = 400
A = 100
A 200
A 400
MNIST// RMSE
2.53 + 0.98
2.47 + 0.98
2.41 0.93
2.16 + 0.02
2.08 + 0.02
1.99 0.02
MNIST/ / MAE
1.61 0.97
1.55 + 0.99
1.51 = 0.96
1.11 + 0.02
1.04 # 0.02
0.96 + 0.01
MNIST NLPD
2.99 + 1.41
2.92 + 1.38
2.84 1.31
2.33 - 0.03
0.03
2.26 + 0.02
2.17 0.02
EMNIST / RMSE
2.37 + 0.95
2.31 + 0.94
2.25 = 0.90
1.99 + 0.17
1.88 + 0.20
1.85 + 0.13
FMNIST MAB 1
1.48
0.91
1.42 0.91
1.39 = 0.89
1.11 - 0.02
1.02 + 0.03
0.98
0.02
2.76
1.33
1.31 2.71
2.65
2.16 0.18
2.07 +
2.04 0.12
CIPARIO RMSE
2.66 + 1.08
2.55 + 1.06
2.55 = 1.03
2.74 + 1.07
2.64 1.08
2.57 + 1.02
CIPARIO/
1.77
1.06
1.69
1.06
1.84 1.03
1.76 +
1.71 1.03
CIFARTO / NLPD
3.20 + 1.55
3.07 + 1.44
3.32 1.89
3.24 + 1.53
3.14 1.53
3.06 + 1.45
All metrics (x10-)).
5.2
Evaluation metrics
In this section, we are interested in the evaluation of the GP decoder with the standard error metrics used for GPs. In Tab. 1 we provide RMSE, MAE and NLPD for the three datasets considered in this experiment. Interestingly, the performance usually improves with larger active set sizes A, as the SAS model captures captures the underlying correlation of datapoints better and this enhances the approximation. This happens for both deterministic and Bayesian cases. Additionally, we observe that the NL.D metric is generally better for smaller active set sizes A. This also be noticed in Fig. 2, where loss curves have a better convergence for the lowest A.
Classification accuracy. We interested in evaluating the structure of the representation. For this purpose, we trained (one) nearest neighbour classifier on the encoded, two-dimensional latent variables and tested the accuracy using encoded test data. Tab. shows the mean and standard deviations of test the accuracy.
Table 2: Classification accuracy (t) on 2-dim. latent space Z.
MODEL
MNIST
FMNIST
BAYESIAN SAS-GP DRC. (ours)
0.63 0.022
0.63 + 0.020
BAYESIAN GP-IVM
0.18 + 0.033
0.24 + 0.043
VAR
0.54 + 0.026
0.58 + 0.008
Runtime and convergence of SAS approximation. The loss curves obtained for the (sparse) Bayesian GP-LVM model and the Bayesian SAS-GP decoder are shown in Fig. 5. We observe that the two methods scale similarly, but that SAS is faster by notable constant within stochastic optimization. Additionally, we observe the performance of the SAS approximation to closely fit the exact log-marginal likelihood log p(æ) (x) 22) For the plot in Fig. 5, we computed the exact probability for subset of MNIST 4040, 000 samples.
PPytorch (Paszke et al., 2019) and Pyro (Bingham et al., 2019) available public repositories.

---

X10
Bayesion SAS-GP
<101
Bayrssan GP-IVM
X107
SAS log
10.5
11.00
-0.25
0.0
EL.00
Suttitt
Baaiaato buckslad poiza
0.55
-1.00
-1.25
-1.0
201
Ram Time [Lound]
Fun Time bours]
Run Time [houss)
Figure 5: (Left-center) ELBO curves for Bayesian SAS-GP and Bayesian GPLVM for different batch-sizes and active points. (Right) SAS loss function compared with the exact log-marginal likelihood computed. In all curves, 4040, 000 samples of MNIST were considered and five different initializations per batch and A setup.
Conclusion
State-of-the-ae representation learning is generally based on neural networks, as this allows for scaling to large datasets. However, often we want reliable uncertainty estimates from the model and we can achieve these with Gaussian process decoders if we can scale them sufficiently. We have reviewed the main difficulties to obtain decent performances with GP-LVM approaches when applied to large-scale learning, even with inducing variables. Revisiting active set approximations, we considered stochastic viewpoint to approximate the marginal likelihood while simultaneously keeping the model marginalized. We formulated our stochastic active sets (SAS) approach for both deterministic and Bayesian versions of GP decoders. We found that our approach works well with amortization, such that a neural network encoder approximately inverts the GP decoder. While amortization also helps when using inducing points, we found the combination with SAS to be particularly efficient and robust.
Empirically, we illustrated the advantage of our method first on image-based observations, where our approach learns better representations using fewer computational resources compared to inducing point methods. We further demonstrated that our approach easily scales to nearly 10t observations. In this experiment, we found that the learnt representations are qualitatively on par with those attained by comparable autoencoder. This is an important finding as, beyond small datasets, GP decoders generally recovers less useful representations compared with models based on neural networks. From this result, we speculate that improvements in training might be enough to get state-of-the-art representations with GP decoders.
Additional benefits. Besides the empirical benefits demonstrated by SAS in the previous section, we have also observed other practical benefits worth reporting. First, we have observed that SAS easily runs in 32-bit numerical precision, unlike inducing point methods that generally require 64- bits of precision (when reporting running times we consistently used 64 bits). Similarly, the jitter usually added the Cholesky factorization is of less importance in SAS, Second, we note that our implementation is surprisingly free of additional tricks and no numerical heuristics were needed to realize reliable implementation.
Limitations and future work. Stochastic active sets rely the Gaussian likelihood, and this is perhaps the strongest limitation. This works well for continuous data, but many data sources are inherently discrete and this requires suitable likelihood, e,g. the discretized mixture of logistics (Salimans et al., 2017). Having more powerful likelihoods would surely improve the GP decoders, but this requires realization of further developments using SAS.
Future work will focus on applying the SAS approach the supervised setting as well, and building SAS-like methods for discrete data. Other possible directions include extending the decoder with deep kernels (Wilson et al., 2016) to capture more features the data and applying convolutional GPs (Van der Wilk et al., 2017) which are more suited to high-dimensional images.

---

Acknowledgements
The authors want to thank Simon Bartels for the fruitful discussions during the early stages of our investigation. This work was supported by research grants (15334, 42062) from VILLUM FONDEN. This project has also received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement 757360). This work was funded part by the Novo Nordisk Foundation through the Center for Basic Machine Learning Research in Life Science (NNF200C0062606) This work was further supported by the Pioneer Centre for AI, DNRF grant number Pl.
References
M. Bauer, M. van der Wilk, and C. E. Rasmussen. Understanding probabilistic sparse Gaussian process approximations. In Advances neural information processing systems (NIPS), pages 1533- 1531,541,20 2016.
E. Bingham, Chen, M. Jankowiak. F. Obermeyer. N. Pradhan, T. Karaletsos, R. Singh, P. A. Szerlip. Horsfall, and Goodman. Pyro: Deep universal probabilistic programming. Journal of Machine Learning Research (JMLR), 20:28:1-28:6, 2019.
D. Bui and R. Turner Stochastic variational inference for Gaussian process latent variable models using back constraints. In Black Box Learning and Inference NIPS Workshop, 2015.
L. Csató and M. Opper. Sparse on-line Gaussian processes. Neural Computation, 14(3):641-66 2002.
E. Fong and C. C. Holmes, On the marginal likelihood and cross-validation. Biometrika, 107(():489-44 496, 2020.
Y. Gal, Y. Chen, and Z. Ghahramani. Latent Gaussian processes for distribution estimation of multivariate categorical data. In International Conference Machine Learning (ICML), pages PMLR, 2015.
Geisser and W. F. Eddy. predictive approach to model selection. Journal the American Statistical Association, 74(365):153-16,, 1979.
Germain, Bach, Lacoste, and Lacoste-Julis Julien. PAC-Bayesian theory meets Bayesian inference. Advances Neural Information Processing Systems (NeurIPS), 29, 2016.
J. Hensman, N. Fusi, and N. D. Lawrence. Gaussian processes for big data. In Uncertainty in Artificial Intelligence (UAI), pages 282-290, 2013.
J. Hensman, A. Matthews, and Z. Ghahramani. Scalable variational Gaussian process classification. In Artificial Intelligence and Statistics (AISTATS), pages 351-360, 2015.
D. P. Kingma and J. Ba. Adam: method for stochastic optimization. International Conference Learning Representations (ICLR), 2015.
P. Kingma and Welling, Auto-encoding variational Bayes, arXiv preprint arXiv:1312.6114, 2013.
Krizhevsky. Learning multiple layers features from tiny images. Technical report, University of Toronto, 2009.
Lalchand, A. Ravuri, and N. D. Lawrenee. Generalised GPLVM with stochastic variational inference. International Conference on Artificial Intelligence and Statisties (AISTATS), pages 7841-7864 2022.
N. Lawrence. Probabilistic non-linear principal component analysis with Gaussian process latent variable models. Journal Machine Learning Rexearch (JMLR). 6(11), 2005.
N. D. Lawrence. Gaussian process latent variable models for visualisation of high dimensional data. In Advances Neural Information Processing Systems (NIPS), pages 329-33 336, 2004.
D. Lawrence and Quiftonero- Ciiiooner--Canela. Local distance preservation the GP-LVM through back constraints. In International Conference on Machine Learning (ICML), pages 513-520, 2006.
LeCun, L. Bottou, Y. Bengio, and Haffner. Gradient-basse learning applied to document recognition. Proceedings of the IEEE, 86(1 1):278 2324, 1998.
S. Liu, X. Sun, P. J. Ramadge, and R. P. Adams. Task-agnosti amortized inference of Gaussian process hyperparameters. Advances in Neural Information Processing Systems (NeurIPS),
2020.
10

---

K. Märtens, K. Campbell, and C. Yau. Decomposing feature-level variation with covariate Gaussian process latent variable models. In International Conference Machine Learning (ICML), pages 4372-4381. PMLR. 2019.
A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, et al. Pytorch: imperative style, high-performance deep learning library. Advances in neural information processing systems (NeurIPS), 32, 2019.
Cichoonrr--Candell and C. F. Rasmussen. unifying view of sparse approximate Gaussian process regression. Journal of Machine Learning Research (JMIR), 19999)99939-1999, 2005.
C. Rasmussen and Z. Ghahramani. Occam's razor. Advances in Neural Information Processing Systems (NIPS), 13, 2000.
C. E. Rasmussen and C. K. Williams. Gaussian processes for machine learning, volume 2. MIT Press, 2006.
D. Rezende and S. Mohamed. Variational inference with normalizing flows. In International Conference on Machine Learning (ICML), pages 1530-1538. PMLR, 2015.
H. Robbins and S. Monro. A stochastic approximation method. The Annals of Mathematical Statistics, pages 400-407, 1951.
Salimans, A. Karpathy, X. Chen, and D. Kingma. PixcICNN++: Improving the PixeICNN with Dis- cretized Logistic Mixture Likelihood and Other Modifications. In International Conference on Learning Representations (ICLR), Toulon, France, 2017.
M. W. Seeger, C. K. Williams, and N. D. Lawrence. Fast forward selection to speed up sparse Gaussian process regression. In International Workshop Artificial Intelligence and Statistics (AISTATS), pages 254-261. PMLR, 2003.
R. Shu, H. H. Bui, S. Zhao, M. Kochenderfer, and S. Ermon. Amortized inference regularization. Advances in Neural Information Processing Systems (NeurIPS), 31, 2018.
A. J. Smola and Bartlett. Sparse greedy Gaussian process regression. In Advances in Neural Information Processing Systems (NIPS), pages (19-625. 2001.
E. Snelson and Z. Ghahramani. Sparse Gaussian processes using pscudo-inpure In Advances in Neural Information Processing Systems (NIPS), pages 1257-1264 2006.
Sohl-Dickstein, Weiss, N. Maheswaranathan, and S. Ganguli. Deep unsupervised learning using nonequilib- rium thermodynamics. In International Conference on Machine Learning (ICML), pages 2256 2256-2265. 2265.2265222 PMLR, 2015.
M. F. Tipping and C. M. Bishop. Probabilistic principal component analysis. Journal of the Royal Statistical Society: Series (Statistical Methodology). 61()):611-62 1999.
M. Titsias. Variational lcarning of inducing variables in sparse Gaussian processes. In Artificial Intelligence and Statistics (AISTATS), pages 567-574, 2009.
M. Titsias and N. D. Lawrence. Bayesian Gaussian process latent variable model. Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics (AISTATS), pages 444-851,20 2010.
M. Van der Wilk, C. F. Rasmussen, and J. Hensman. Convolutional Gaussian processes. In Advances in Neural 2017.
Information Processing Systems (NIPS), pages
A. Vehtari and Lampinen. Bayesian model assessment and comparison using cross-validation predictive densitics. Neural Computation, 14(10):449-2468 2002.
G. Wahba. Spline models for observational data. SIAM, 1990.
A. Wilson, Z. Hu, R. Salakhutdinov,
Xing. Deep kernel learning. Artificial Intelligence and Statistics (AISTATS), pages 370-378. PMLR, 2016.
H. Xiao, K. Rasul, and R. Vollgraf. Fashion-MNIST: A novel image dataset for benchmarking machine learning algorithms. arXiv preprint |rXiv: 2017.

---

Checklist
1. For all authors:
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contribu- tions and scope? [Yes] We put particular attention on active sets, scalability and the connection between marginal likelihood of GPs and CV.
(b) Did you describe the limitations your work? [Yes] See Section 6.
(c) Did you discuss any potential negative societal impacts your work? [N/A]
(d) Have you read the cthics review guidelines and ensured that your paper conforms them? [Yes]
2. If you are including theoretical results:
(a) Did you state the full set of assumptions of all thcoretical results? [Yes] In Section 3 and the
Appendix.
(b) Did you include complete proofs all theoretical results? [Yes] See Appendix.
3. If you experiments:
(a) Did you include code, data, and instructions needed to reproduce main experimental results (cither in the supplcmental material or as a URL)? [Yes] We provide an URL for the Pytorch code, which will be available in public repository. Additionally, experiments files will be included and details the data can be found in the Appendix.
(b) Did you specify all the training details (e.g., data splits, hyperparameters. how they were chosea)? [Yes] Set up of GP decoder variables (A. initial 0, etc) used in the experiments, further details on initialization, algorithms and optimization are provided in the Appendix.
(c) Did you report error bars (e-g., with respect the random seed after rumning experiments multiple times)? [Yes] our experiments consider multiple initialization and seeds randomness. which are reflected on the metrics and figures.
(d) Did you include internal cluster, or cloud provider)? [Yes] A deseription of computational resources used is provided in the Appendix.
total amount of compute and the type resources used (e.g.. type GPUs,
If you are using existing assets (e-g., code, data, models) or curatingfreleasing new assets:
(a) If your work uses existing assets, did you cite the creators? [Yes] In the experiments and in the scripts.
(b) Did you mention the license of the assets? Yes) In the appendix and scripts.
(c) Did you include any new assets either in the supplemental material or as URL? [N/A]
(d) Did you discuss whether and how consent was obtained from people whose data you're 're us- ing/curating? [N/A]
(c) Did you discuss whether the data you are using/curating contains personally identifiable informa- tion offensive content? [N/A]
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if applicable?
[N/A]
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, if applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent participant compensation? [N/A]
12


Reviews:
[{"review_id": "YR8CmX0DX1M", "reviewer": "Reviewer_38P8", "summary": "The authors investigate a novel stochastic gradient ascent method for fitting Gaussian Process Latent Variable Models. In practice, if I understood correctly, the method samples a \"minibatch\" from the full data to build a GP, and measures posterior likelihood on the remaining hold out data. This is justified by recent mathematical results that show the full log marginal likelihood ($O(n^3)$ cost) is equivalent to the marginal likelihood over subsets ($O(b^3)$ with batchsize $b$) of all subset sizes and ( $O(n-b)$, linear cost) summations over terms excluded from the subset. This alternative marginal likelihood evaluation is a a sum and an expectation and is amenable to stochastic Monte-Carlo approximation, and therefore unbiased stochastic gradient ascent.\n\nThe authors apply this trick to both standard and Bayesian GP-LVMs implemented both with and without amortization over a range of datasets and demonstrate the benefits.", "questions": "- in the Standard GP-LVM case, can the authors confirm that the only practical novelty (over naïve subsampling data) is the random batch size and additional term in Equation (6) and (8)? (If so, this can be portrayed as a strength)\n- Can the authors include such a naive subsampling data baseline forming an ablation study?", "limitations": "The authors clearly discuss the Gaussian likelihood assumption. No more unaccounted limitations come to my mind.", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "**Strengths**\n- I felt the paper was very well written, providing a clear story of GP-LVMs and their milestones over the years and with clear interpretation\n- the result of Fong and Holmes 2020 appears to show the marginal likelihood can be (stochastically) approximated using just partitions of data, which is very impactful for GPs that require cubic cost effectively reducing cost to smaller cubes and linear costs.\n- I felt the application to GP-LVMs and Bayesian prior GP-LVMs with and without amortization covered a wide extensive range of empirical use cases and provides value to the community.\n\n**Weaknesses**\n- I feel that evaluating marginal likelihood on a stochastically subsampled data to reduce the cost of GP training is a fairly common trick, but obviously biased. Admittedly I cannot seem to find papers using it exactly, a google search for \"Gaussian Process minibatch\" returned some results [1, 2].\n- Hence the only _practical_ novelty (over naive subsampling) in this work is \n    - (1) the first term of Equations 6 and Equation 8 and\n    - (2) the batch size is random (but it doesn't need to be? See nitpick)\n- The rest of the paper seems to be technical background on GP-LVMS and Fong and Holmes (2020).\n- Can the authors include a naive random subsampling baseline? Simply SAS methods without the First terms of equation 6? This would serve as an ablation to show the impact of the theoretical result.\n\n**Nit Picks**\n- Equation 4: $R$ is said to be averaged over $R$ yet it is a constant on the right term but not in the centre or left terms, this is explained in the Appendix \"Additionally, Eq. (3) holds for every size of the hold-out data R ∈ [1, 2, · · · , N]\" and would make Eqn 4 easier less confusing to me.\n- Equation 3: for me, seeing the bold $x_n$ automatically made me think it was a set of correlated points requiring whose probability would have led to $n^3$ matrix inversion.\n- Appendix Equation 4: I was very curious to know the intermediate steps to deriving $S_{PCV}(x, R)$ and I couldn't find the answer during a (very) quick skim through Fong and Holmes (2020).\n\n\n[1] https://arxiv.org/abs/2111.10461\n\n[2] https://nipunbatra.github.io/blog/ml/2021/09/03/param-learning-sgd.html", "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": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "PuA92F0Q1D", "reviewer": "Reviewer_QibS", "summary": "The paper proposes a computationally efficient way to train the GP-LVM based on a connection between marginal likelihood and cross-validation in Fong and Holmes (2020). To make the estimation practical, stochastic estimation is used for sampling the size of the hold-out set and also the size of the active set. Amortization is also used for approximating $p(z|x)$. In the experiments, the performance of the method is compared with GP-LVM and VAE. ", "questions": "1. How is R picked? Also, will the uniform sampling of r results in high variance of the estimate?\n2. Why does Bayesian GP-LVM perform much worse than the proposed method or VAE? Some explanation and ablation studies would be helpful.", "limitations": "N.A.", "rating": 5, "confidence": 4, "soundness": 2, "presentation": 1, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n- Applying the idea of estimating MLL with cross-validation scores to GP-LVM is interesting.\n- Several ablation studies on different active set sizes.\n\nWeaknesses:\n\n- It is not clear how different hold-out set sizes (R) will affect the proposed estimation. Also, I have difficulty understanding what “by first uniformly sampling R” means. Does it mean first pick R and then uniformly sample $r \\in [1:R]$? I could not infer this from the pseudo-algorithm or experiments as well. If the sampling is by uniformly sampling $r \\in [1:R]$, it seems necessary to have more discussion about whether the sampling is of high variance and if there are ways to reduce the variance.\n- The presentation of the experimental results can be improved. Some important experiment details are deferred to the appendix, such as what kernel, NN architecture is used. They serve as very important information to understand the results. Also, not sure why RMSE/MAE is used as a metric as the MNIST or CIFAR10 are classification tasks. For evaluating the quality of the representation learned, it is not clear why 2 is selected as the number of dimensions. It might be better to consider trying a range of different dimensions.\n- There are two major components in the proposed method: the SAS approximation and amortization. It would be better if there is some ablation study on both SAS and amortization. For example, by comparing GP-LVM-SAS-Amortization with GP-LVM-InducingPoint-Amortization and GP-LVM-SAS with GP-LVM-InudcingPoint without amortization.\n- It seems that the baseline methods are only compared with the proposed method in terms of representation quality. There is no comparison on NLPD/RMSE/MAE. And it would be better to include some regression tasks as well.", "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": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "Jff8P1P0Vu", "reviewer": "Reviewer_YAoZ", "summary": "This paper presents a new approach to train Gaussian process latent variable models (GP-LVMs). The idea relies on the (previously shown) connection between cross-validation (CV) and marginal likelihood (ML). More specifically, the authors interpret the leave-R-out CV objective from an optimization perspective, which allows for stochastic and scalable optimization of the ML objective of vanilla and Bayesian GP-LVM formulations. The experiments demonstrate that the framework leads to more compact latent representations compared to previous formulations. On MNIST, the deterministic variant of the framework is shown to match the exact ML with a certain hyperparameter choice (batch and active set size).", "questions": "Presentation/clarity/notation:\n- What is $K_{NN}$ in (1)?\n- Shouldn't $S_{CV}(x|R)$ in (3) be also conditioned on $z$?\n- How about including the formal definitions of $S_{PCV}$ and $S_{CCV}$ in the main paper instead of supplements?\n- $r$ appears on the middle term in (4) but only $R$ appears on the right hand side. \n- The first equation in line 114 involves $R$ on the left but not on the right. If I understood correctly, $|\\mathcal{A}|=N-R$ (hence $R$ implicitly appears on the right. In that case, I would suggest using modifying the notation a bit to reflect this.\n- Are $\\mathcal{A}$ and $\\mathcal{A}_p$ used interchangeably?\n\nExperiments:\n- The loss curves in Figure 2 oscillate a lot.\n- The MNIST latent space with Bayesian SAS is missing.\n- Only $2D$ latent spaces are considered, which significantly restricts the generalizability of the findings.\n- NLPD consistently increases with active set size (for which I would use the shorthand $|\\mathcal{A}|$ instead of $A$). Is this just the opposite of what we would expect?\n\nMisc:\n- Although the proposed framework seems suitable for GP models in general, the authors particularly focus on GP-LVMs. I wonder why authors restrict themselves to GP-LVM. It would also be nice to see how the framework performs on a standard GP regression task with simple noisy inputs and outputs. ", "limitations": "Everything sounds good.", "rating": 5, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Overall, I liked the paper as it nicely brings theoretical and practical aspects together. My major comments concern the derivation/notation and experiments, which strangely both escalate and diminish my score.\n\nStrengths:\n- The framework has solid theoretical foundations. I find the connection with CV important as it distinguishes the proposed framework from heuristics-based optimization objectives.\n- The results are impressive. The improvement over existing deterministic and Bayesian GP-LVMs is very clear. Looking at the latent nearest neighbor classification experiment, one can even argue that the learned latent space is more structured than that of VAE.\n\nWeaknesses:\n- The presentation can be improved. I list several points in the subsequent \"Question\" section. Most importantly, I do not understand\n  - how to obtain (5) from (4), and\n  - how $S_{PCV}$ in (5) has $x_A$ on the left hand side of the conditioning while (3) has $x_A$ on the right.\n- The claims in lines 29-34 are not well-supported. It would be much better to list the problems of training with inducing points (with references) instead of quotations.\n- The experiments are slightly worrying (please see the next \"Question\" section again)\n\nI would be glad to increase my score if my concerns are addressed (even verbally).\n\n**UPDATE:** After the author response, I increase my score to 5.\n\n", "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": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept