Input Dependent Sparse Gaussian Processes 
Authors: Bahram Jafrasteh
Carlos Villacampa-Calvo
Daniel Hernández-Lobato
Venue: ICLR 2022 Submitted (2022)
Keywords: gaussian processes
variational inference
neural networks
sparse approximations

Abstract: Gaussian Processes (GPs) are Bayesian models that provide uncertainty estimates associated to the predictions made. They are also very flexible due to their non-parametric nature. Nevertheless, GPs suffer from poor scalability as the number of training instances $N$ increases. More precisely, they have a cubic cost with respect to $N$. To overcome this problem, sparse GP approximations are often used, where a set of $M \ll N$ inducing points is introduced during training. The location of the inducing points is learned by considering them as parameters of an approximate posterior distribution $q$. Sparse GPs, combined with variational inference for inferring $q$, reduce the training cost of GPs to $\mathcal{O}(M^3)$. Critically, the inducing points determine the flexibility of the model and they are often located in regions of the input space where the latent function changes. A limitation is, however, that for some learning tasks a large number of inducing points may be required to obtain a good prediction performance. To address this limitation, we propose here to amortize the computation of the inducing points locations, as well as the parameters of the variational posterior approximation $q$. For this, we use a neural network that receives the observed data as an input and outputs the inducing points locations and the parameters of $q$. We evaluate our method in several experiments, showing that it performs similar or better than other state-of-the-art sparse variational GP approaches. However, with our method the number of inducing points is reduced drastically due to their dependency on the input data. This makes our method scale to larger datasets and have faster training and prediction times. 

Under review as a conference paper at ICLR 2022
INPUT DEPENDENT SPARSE GAUSSIAN PROCESSES
Anonymous authors Paper under double-blind review
ABSTRACT
Gaussian Processes (GPs) are Bayesian models that provide uncertainty estimates associated to the predictions made. They are also very flexible due to their non- parametric nature. Nevertheless, GPs suffer from poor scalability the number of training instances N increases since their cost is cubic in N. To overcome this problem, sparse GP approximations are often used, where set of M < N inducing points is introduced during training. The location of the inducing points is learned by considering them as parameters of approximate posterior distribution q. Sparse GPs, combined with variational inference for inferring q, reduce the cost of GPs per iteration (MM)) Critically, the inducing points determine the flexibility of the model and they are often located in regions of the input space where the latent function changes. A limitation is, however, that for some learning tasks a large number of inducing points may be required to obtain good prediction performance. To address this limitation, we propose here to amortize the computation of the inducing points locations, the variational posterior approximation q. For this, we use a neural network that receives the observed data and the parameters of q- We evaluate our method in several experiments, showing that it performs similar or better than other state-of-the-art sparse variational GP approaches. However, with our method the number of inducing points is reduced drastically due their dependency on the input data. This makes our method scale to larger datasets and have faster training and prediction times.
well as the parameters of
input and outputs the inducing points locations
INTRODUCTION
Gaussian Processes (GPs) are non-parametric models that can be used to address regression and clas- sification machine learning problems (Rasmussen & Williams 2006). GPs become more expressive as the number of training instances N grows and, since they are Bayesian models, they provide pre- dictive distribution that estimates the uncertainty associated to the predictions made. This uncertainty estimation or ability to know what is not known is critical in many practical applications (Gal) 2016). Nevertheless, GPs suffer from poor scalability as their training cost is O(N3) per iteration due to the need of computing the inverse of N x N covariance matrix. Another limitation is that approximate inference is needed with Non-Gaussian likelihoods (Rasmussen & Williams 2006).
Sparse approximations can improve the cost of GPs (Rasmussen & Williams, 2006). The most popular ones introduce set of M < N inducing points (Snelson & Ghahramani] 2006 Titsias 2009). The inducing points and their associated posterior values completely specify the posterior process at test points. In Snelson & Ghahramani (2006), the computational gain is obtained by assuming independence among the process values at the training points given the inducing points and their values. This can also be seen as using an approximate GP prior (Quiñonero-Canee Rasmussen 2005). By contrast, in Titsias (2009) the computational gain obtained by combining variational inference (VI) with posterior approximation q that has a fixed part and tunable part. In both methods the cost is reduced to (N M²) per iteration and the inducing points, considered as model's byper-parameters.
&
learned by maximizing estimate of the marginal likelihood.
Importantly, the VI approach of Titsias 2009) maximizes lower bound on the log- marginal likelihood as indirect way of minimizing the KL-divergence between approximate posterior distribution for the process values at the inducing points and the corresponding exact posterior. The advantage is that the objective is expressed as sum over the training instances, allowing for

---

Under review as a conference paper at ICLR 2022
mini batch training and stochastic optimization techniques to be applied the objective (Hensman et al. 2015b). This reduces the cost to ( per iteration, making GPs scalable to large datasets.
In sparse approximations one often observes in practice that after the optimization process the inducing points are located in regions of the imput space in which the latent function changes (Sncison Hensman al.] 2015a Bauer et al., 2016). Therefore, the expressive power of the model depends on the number of inducing points M and their correct location
& Ghahramani 2006 Titsias 2009
on the input space. Some problems may require large number of inducing points, in the order of
thousands, to get good prediction results (Hensman et al. 20155 Shi et al. 2020; Tran et al. 2020)
This makes training inducing point based sparse GPs diffcult in those problems.
There have been some attempts to improve the training cost of sparse approximations, including using different sets of inducing points for the computation of the posterior mean and variance (Cheng & Boots, 2017). Other approaches use an orthogonal decomposition of the GP that allows to introduce an extra set of inducing points with less cost (Shi et al., 2020). Finally, other methods consider a large set of inducing points, but restrict the computations for a particular data point to the nearest neighbors to that point from the set of inducing points (Tran et al. 2020).
In this work we are inspired by Tran et al. 2020) and propose novel method to improve the training cost of sparse GPs. Our method also tries to produce a set of inducing points (and associated variational approximation q) that are specific of each input data point. For that, note that some works in the literature have observed that one can learn the mappings from inputs to proposal distributions instead of directly optimizing their parameters (Kingma & Welling 2014; [hhu et al. 2018). This approach, known as amortized variational inference, is a key contribution of variational auto-encoders (VAE) (Kingma & Welling, 2014), and has also been explored in the context of GP to solve other types of problems such as multi-class classification with input noise (Villacampa-Calvo ct al. 2021.. Amortized inference has also been empirically shown to lead useful regularization properties that improve the generalization performance (Shu et al. 2018).
Specifically, here we combine sparse GPs with neural network architecture that computes, for each potential data point, the associated inducing points to be used for prediction. We also employ a neural network to carry out amortized VI to compute the parameters of the variational distribution q approximating the posterior of the process values for the outputted inducing points. While the number of parameters that need optimization may increase with the use of neural network, this approach allows for big reduction in the total number of inducing points without losing expressive power. In particular, it enables different sets of inducing points associated to each input location. The inducing points are simply given by a mapping from the inputs provided by neural network. We show on several experiments that the proposed method is able to perform similar or better than standard sparse GPs and competitive methods for improving the cost of sparse GPs (Tran et al. 2020 Shi et al. 2020). However, the training and prediction times of our method are much better.
GAUSSIAN PROCESSES
Gaussian Process (GP) is a stochastic process for which any finite set of variables has a Gaussian distribution (Rasmussen & Williams 2006). In learning task, use GP prior over latent function. Then, Bayes' rule is used to get a posterior for that function given the observed data. Consider a dataset D = (xx, ) where each scalar y; is assumed to be obtained as y = f(x)+ + e, with f(-) latent function and Gaussian noise with variance o², i.e., E, NO, o²). We specify a prior distribution for f in the form of GP, which is described by mean function m(x) (often set to zero) and covariance function k(x,xx such that k(xx x') = E[[f (x) - m mfx))(((xx) m(x)))] Covariance functions typically have some parameters O. Given D. the predictive distribution for f at new test point x* is Gaussian with mean and variance given by
[(x*) k(x ) + o21)-ly,
o²
k(xx
k(x " + 021) kx("),
(1)
where u(x*) and 0² (x*) ) are the prediction mean and variance, respectively. k(x*) is vector with the covariances between f(x") and each fxx Similarly, K has the covariances between f(x)) and f(x,) for for = 1, .,....J.. .N. Finally, I stands for the identity matrix. Popular covariances functions k(.,)) are the squared exponential and the Matérn (Rasmussen & Williams 2006). Their parameters, 0, and o² can be found by maximizing p(y) (Rasmussen & Williams 2006). The computational
2

---

Under review as a conference paper at ICLR 2022
complexity of this approach is ((N2) since it needs the inversion of K, N N matrix. This makes GPs unsuitable for large data sets.
2.1 SPARSE VARIATIONAL GAUSSIAN PROCESSES
Sparse approximations improve the cost of GPs. The most popular methods introduce, in the same input space as the original data, new set of M & N points called the inducing points, denoted by Z (Zi..... ZM) )d Snelson & Ghahramani 2006 Titsias 2009). Let the corresponding latent function values be observed data and their location can be learned during training. A GP prior is placed on u. Namely, p(u) ~ (O0 Kz), where Kz is a matrix with the covariances associated to each pair of points from Z. The idea that the posterior for f can be approximated in terms of the posterior for u.
(f(z),..... ((z)) )) The inducing points are not restricted to be part of the
In this work we focus on a widely used variational inference (VT) approach to approximate the posterior for f [Titsias] 2009). f= ((()), [ (xx))) VI, the goal is to find an approximate posterior for f and g(f, [ u). that resembles
much as possible the true posterior p(f, uly). Critically, q is constrained to be q(f, u) = p(f))))((u) with p(flu) fixed and q(u) a tunable multi- variate Gaussian. To find q()) lower bound of the marginal likelihood is maximized. The evidence lower bound (or ELBO) is obtained via Jensen's 's inequality, leading to (after some simplifications):
C= C= (logg p(us f(] (L(g(u))p(u)),
(2)
where f(y) is the model's likelihood for the i-th point and KL[-I-] the Kullback-Leibler di- vergence between probability distributions. In Titsias (2009). they optimize q(u) in foored-form. The resulting expression is then maximized to estimate Z, 0 and d². This leads to a complexity of O(NM²). However, if the variational posterior g(u) is optimized alongside with Z, 0 and 0², as pro- posed in Hensman et al. (2013), the ELBO can be expressed as sum over training instances, which allows for mini-batch training and stochastic optimization techniques. Using stochastic variational inference (SVI) reduces the traiming cost to (M33) per iteration (Hensman et al. 2013). Importantly, the first term in (2) is expectation that has closed-form solution in the case of Gaussian likelihoods. It needs to be approximated for other cases, egg., binary classification, either by quadrature or MCMC methods (Hensman et al.) 2015b). The second term is the KL--iiergeence between the variational posterior and the prior, which can be computed analytically since they are both Gaussians.
The expressive power of the sparse GP critically depends on the number of inducing points M and on their correct placement in the input space via optimizing (2) (Titsias) [2099] Hensman et al. 501 Bauer al. 2016). In some problems several thousands of inducing points may be required get good results (Hensman et al. [2015b] Shi et al.) [2200] TTran et al.) 2020). This makes difficult and expensive using sparse GPs in those problems. In the next section we describe how alleviate this.
INPUT DEPENDENT SPARSE GPs
We develop new formulation of sparse GPs which for every given input computes the corresponding inducing points to be used for prediction, and also the associated parameters of the approximate distribution q. To achieve this, we consider a meta-point x that used to determine the inducing points Z and the corresponding u. Namely, now u depends on X, i.e., u p(u(x). In particular, we set p(u(x) = (O, Kz(z)) where the inducing points Z depend non-linearly, e.g., via a deep neural network, x. The joint distribution of and x is then given by p(u, x) p(u pxx())p(x) for some prior distribution p( ()) Following (Tran et al. 2020), we can consider an implicit distribution p(x).
That is, its analytical form is unknown, but we can draw samples from it. Later on, we specify p(x).
Note that the marginalized prior p(u) is no longer Gaussian. However, we can show that this formulation does not impact on the prior over f. For an arbitrary selected meta-point we have that
u(x) =N (0.
K
Kx,z(8) K
1)
(3)
K K
where Kx,z(3) are the cross-covariances between f and u. Therefore, if us marginalized out in (3) the prior for f is the standard GP prior and does not depend on x. Hence, p(f[x) = p(f). Thus, p(f, u) = p(f, ul))p((X)d) dx is a mixture of Gaussian densities, where the marginal over is the same for every component of the mixture. In the standard sparse GP, the inducing points also have an
3

---

Under review as a conference paper at ICLR 2022
impact the variational approximation via the fixed conditional distribution uf (Titsias 2009). Therefore, we also incorporate the input dependence on x in q. This is done in the next section.
3.1 LOWER BOUND ON THE LOG-MARGINAL LIKELIHOOD
We follow Tran et al. 2020) derive a lower bound the log- marginal likelihood of the ex- tended model described above. Consider a posterior approximation of the form g(f, u, x) = p( ([[]))q(u(X)) (p((()/())(()))))))))))))))))( where only q(u|&) can be adjusted and the other factors are fixed. Using this posterior's factorization and Jensen's 's inequality we obtain the lower bound after some simplifications:
L = ii
E J ((x) (p(falu)g(
log p(ue fa)dfdu N (L (t
[p([[))] dx.
(4)
Now, assuming that p(x) is an implicit distribution, we can draw samples from it and approximate the expectation w.r.t p(x). Thus, for meta-point sample x, from p(x). (4) is approximated as
R Ei [Eppf// [log (())
We can evaluate 5 and its gradients to maximize the original objective in (4) using stochastic optimization techniques. This is valid for any implicit distribution p(x). Consider now that we use mini-batch-bret based training for optimization, and we set x, x,. In this case, the value of x remains random, as it depends on the points (x, yi) that are selected in the random mini-batch. This results in method that computes different inducing points for each input location. In practice, we the same sample to approximate the expectation w.r.t. p(*) and the sum across the data in (5) This could introduce a bias the objective. However, such reusing of the samples done in Tran with good empirical results. Moreover, our experiments in Section also validate this approximation.
(2020)
3.2 AMORTIZED VARIATIONAL INFERENCE AND DEEP NEURAL NETWORKS
Maximizing the lower bound finds the optimal approximate distribution q. problem, how- ever, is that we have potential large number of parameters to fix, corresponding to each x.). In particular, if we set q(u[x,) to be Gaussian, we will have to infer different means and covariance matrices for each differ- ent x,. This is expected to be memory inefficient and to make optimization more difficult. To re-
duce the number of parameters of our method use amortized variational inference and spec- ify a function that can generate these parame- ters for each x, (Shu et al. 2018). More pre- cisely, we set the mean and covariance matrix of g(u|x,) xx) to be mxx;) and S(x;), for some non- non-linear functions.
Deep neural networks (DNN) are flexible mod- els that can describe complicated functions. In these models, the inputs go through several lay- ers of non-linear transformations. We use these models to compute the con-liieasss that gen- erate from x; the inducing points, Z(x;), and the means and covariances of q( ((|xx; xx). ), L..., m(x;) and S(x,). The architecture employed is displayed in Figure [ At the output of the DNN we obtain Z, mean vector and the Cholesky factor of the covariance matrix S LL. The maximization of the lower bound in (4) when using DNNs for the non-linearities is shown in Algorithm [] The required expectations are computed in closed-fore in regression. In binary classification, we use I-dimensional quadrature. as in Hensman et al. 2015a).
Figure 1:
The network's inputs is x. The outputs are the inducing points Z, the mean vector m and the Cholesky factor, L, of g(u|x;).
3.3 PREDICTIONS AND TRAINING COST
At test the instances are not randomly chosen. In that case, we simply set p(x) be a deterministic distribution placed on the candidate point x*. The DNN is used to obtain the associated information. Namely, Z, and the parameters of q(u(x*), m and S. The predictive distribution for I f((x*) is:
f(xx) ~ N (Kx. .Kz m, 1(xx _oxxxx'> x") + Kx 2 2Kz (S Kz) K KT z)
(6)
4

---

Under review as a conference paper at ICLR 2022
Algorithm 1 Training input dependent sparse GPs
Input: D, .M. neural network NNet with L hidden layers and hidden units
Output: Optimal parameters of the model
initialize neural network's weights and kernel's parameters 0 while stopping criteria is False do
Logix = 0, KL =
gather mini-batch Mb of size from D
for (x) in Mb do
(Zx; My, ]x,] NNet(xi): Logi: (B( wl ploggllg/lll p(x]][]] :KLaiv KLlg(u/x))llullss
ELBO + N x Logn _ x KLdiv
Update parameters of the model using the gradient of ELBO
Given this distribution for fxx",, the probability distribution for y" can be computed in closed form in the case of regression problems and with 1-dimensional quadrature in the case of binary classification. Note that (6) is only suitable for predictions at individual test points, as in Tran et al. (020). This can be a limitation in applications needing covariances. As solution, one could consider for prediction the union of all input dependent inducing points for test points. This would be inconsistent with the proposed training method. Nonetheless, such approach can also be modified to consider the union of input dependent inducing points for the training points within mini-batch, in (Tran et al. (020).
The cost of our method is smaller than the cost of a standard sparse GP if a smaller number of inducing points M is used. The cost of DNN with L layers, P hidden units, d, dimension of the input data, and output dimension d, is Onn, P + n n2P2i. + n nPd, + n(T. + 1)). The cost of the sparse GPs is (n.M), with n the mini-batch size. Therefore, the cost of our method per iteration is O(nd P + rnP2 L + nPd, + nlL + 1) + ") Since in our method the inducing points are input dependent, we expect to obtain good prediction results even for M values that are fairly small.
RELATED WORK
Early works on sparse GPs simply chose a subset of the training data for inference based on an information criterion (Csató & Opper 2002. Lawrence et al.] 2003) Seeger et al.) [003 Henao & Winther 2012). This approach limited practice and more advanced methods in which the inducing points need not be equal to the training points are believed to be superior. In the literature there are several works analyzing and studying sparse GP approximations based on inducing points.
Some of these works include Duinonero-Cannee Candela & Rasmussen (2005): [Snelson & Ghahramani & Hernández-Lobato/ (2016). We focus here on variational approach (Titsias 2009) which allows for stochastic optimization and mini-batch training (Hensman et al., 2013) 2015a). This enables learning in very large datasets with cost of ((3)) per iteration, with .7 the number of inducing points.
(2006); (oi06) Slaish-Guzman & Holden y; Titsias (2009 Bauer et al. (2016); Hernáneez-
In some problems, however, several thousands of inducing points may be needed to get good prediction results (Hensman et al. 2015b; Shi al. 2020 Tran al. 2020). There is hence need to improve the cost of sparse GPs, without losing expressive power. One work addressing this task is that Cheng Boots 2017). In that work is proposed to decouple the process of inferring the posterior mean and variance, allowing to consider a different number of inducing points for each one. Importantly, the computation of the mean has linear complexity, which allows to have more expressive posterior means at a lower cost. A disadvantage is that such an approach suffers from optimization difficulties. An alternative decoupled parameterization adopts an orthogonal basis in the mean Salimbeni et al.] (2018a). Such method can be considered as specific case of Shi et al. (2020). There, the authors introduce a new interpretation of sparse variational approximations for GP using inducing points. For this, the GP is decomposed leads to tighter lower bounds on the marginal likelihood and new inference algorithms considering two different sets of inducing points. This enables using more inducing points at linear cost.
of two independent processes. This
Our work is closer to Tran et al. (2020). There, mechanism is also described to consider input dependent inducing points in the context of sparse GP. However, the difference is significant. In particular, in Tran et al. (2020) a very large set of inducing points M is considered initially. Then, for each input point, subset of these inducing points is considered. This subset is obtained by finding
5

---

Under review as a conference paper at ICLR 2022
the K < M nearest inducing points to the current data instance x,, This approach significantly reduces the cost of the standard sparse GP described in Titsias (2009). However, it suffers from the difficulty of having to find the K nearest neighbors for each point in a expensive. Therefore, the final cost is higher than what would be thought initially. Our method is expected to be better because of the extra flexibility by the non-linear relation between x, and Z given by the DNN. Furthermore, the DNN is expected to make a better use of GPU acceleration.
which is very
Another method to improve the training cost of GP is described in Wilson & Nickisch (015); [Evans & Nair (2018); Gardner et al. (2018). It consists in placing the inducing points on a grid. This allows to perform fast computation exploiting the inducing points structure. One can easily consider values for M that are even larger than N. However, to get such benefits the inducing points need to be fixed due to the structure constraints. This may be detrimental in high-dimensional problems.
Instead of using inducing points, there are some works that scale GPs by approximating the posterior GP process using an inference network (Shi et al., 2019) [unn et al. 2019). An inference network receives some random noise and outputs function values for each input. Particular examples include among others Bayesian DNNs. Inference networks are expected to lead to flexible stochastic processes. However, it is difficult to enforce that the approximate posterior process looks similar to the prior GP in regions where there is no data. For this, approximate inference is carried out on a finite subset of points chosen at random from the input space. This is expected to lead to poor results in high-dimensional spaces. Moreover, another problem of using an inference network is that tuning the prior GP hyper-parameters, is challenging and has often to be done in separate step.
Amortized variational inference (Shu et al. 2018) has also been explored in the context of GPs in Villacampa-Calvo et al. (202T). There, input noise is considered in multi-class learning problem and the performance of the final model is improved by amortizing the variational parameters of the posterior approximation for the noiseless inputs, using a DNN that receives both x; and yi-
Other sparse GPs in the literature do not fully rely on inducing points, e.g., (Trsp) 2000; Snelson 2007 Gramacy & Apley] 2015). These techniques. however, cannot use, in general, stochastic optimization and do not scale to very large problems. Finally, sparse GPs, and our method, can benefit from natural- -raaiinnt (Salimbeni et al., 2018b). They could result in an orthogonal improvement.
EXPERIMENTS
We evaluate the performance of the proposed method, to which we refer to as Input Dependent Sparse GP (IDSGP). We consider both regression and binary classification with probit likelihood. In this later case, we approximate the expectation in the lower bound using I-dimensional quadrature, as in
Hensman al. (2015a). The code of the proposed method in Tensorflow 2.0 (Abadi et al.] 2015) is provided in the supplementary material. In the experiments we compare results with the standard variational sparse GP (Titsias 2009). We refer such a method as VSGP. We also compare results with two of the methods described in Section 4 Namely, the sparse within sparse GP (SWSGP) described in Tran al. 2020), and the sparse GP based on an orthogonal decomposition that allows to consider two different sets of inducing points (Shi et al., (020). We refer to this last method as SOLVE. All methods use Matérn 3/2 covariance function (Rasmussen & Williams 2006). The DNN architecture employed in IDSGP is described in detail in Appendix B
5.1 TOY PROBLEMS
We show the posterior mean and standard deviation of each method on a 1-dimensional regression problem (Snelson & Ghahramani 2006). We compare results with full GP. Figure shows the results obtained, including the learned locations of the inducing points. In the case of IDSGP we show the locations of the inducing points for the point represented with a star at x = 3.9. The number of inducing points, for each method, are indicated in the figure's 's caption. We consider a small number of inducing points to study the benefits of having input dependent inducing points. IDSGP uses smaller number of inducing points than the other methods. The figure shows that, in regions with observed data, IDSGP's predictions look closer to those of the full GP. Appendix [.] has results for an increasing number of inducing points M. They show that as M increases IDSGP becomes more similar to the full GP. Figure 3 shows the decision boundary of each method on the banana classification dataset (Hensman et al. 2013). IDSGP produces the most accurate boundaries.

---

Under review as a conference paper at ICLR 2022
VSGP
IDSGP
menn precktion
mean presktion
menn SP prudiction
mean podditinn
Standart dediation
Intual InSuuny perto
indtl incucing porto Inducing polnes
Incucing poins paint
SWSGP
SOLVE
mean prediction
Suastard decacie
Stansard deciitionn
Standard deckesion oP Thandin data
Thaiing dats
11 ITSLCIn poRtos
Toy data set with N 200 points. Initial and final locations for the inducing points are shown on the top and bottom of each figure. In IDSGP, the inducing points correspond to the point drawn with star. The posterior mean and standard deviation of full GP are shown with blue and brown dashed lines, respectively. VSGP method with M 4. IDSGP with M 2 and a neural network with 2 layers with 50 units. SWSGP with M = 4 and 2 neighbors. SOLVE with M1 = M2 = 2.
Figure 2:
Training data lobserved values)
Text data (observed valueal Inducing point
Induring points
SWSGP
sols
(observed valuesi
Test data (observed values) Induring points Inducing points v
Figure 3: Banana classification data set with N = 5300 points. The final location of inducing points are shown inside the figures. For IDSGP, we show the location of inducing points related to the green colored point. VSGP with M = 4. IDSGP with M = 2 and neural network with 2 hidden layers each contains 50 hidden nodes. SWSGP with M = 4 and 2 neighbors. SOLVE with M1 = M2 = 2.
5.2 EXPERIMENTS ON UCI DATASETS
We consider several regression and binary classification datasets extracted from the UCI repository Dua & Graff 2017). The number of inducing points of IDSGP is set to M 15. In SOLVE we use M, 1024 and M2 = 1024 inducing points. In VSGP we set M = 1024. In SWSGP we set M- 1024 and K 50 neighbors. All the methods are trained using ADAM (Kingma & Ba 2015)
7

---

Under review as a conference paper at ICLR 2022
with mini- batch size of 100 and learning rate of 0.01. In the classification setting setup, but the number of inducing points of IDSGP is set to M = 3. All methods are trained on a Tesla P100 GPU with 16GB of memory. On each dataset use 80% of the data for training and the rest for testing. We report results across 5 splits of the data since the datasets are already quite big.
the same
The average negative test log- likelihood of each method on each dataset is displayed in Table I for the regression datasets, and in Table for the classification datasets, respectively. The average rank of each method is also displayed at the last row of each table. The RMSE and prediction accuracy results are similar to those displayed here. They can be found in Appendix D.3 and Each table also shows the number of instances N and dimensions d of each dataset. We observe that in the regression datasets, the proposed method, IDSGP, obtains best results in 6 out of the 8 datasets. IDSGP also obtains the best average rank (closer to always performing best on each train test data split). This is remarkable given that IDSGP a much smaller number of inducing points (e.g.. M = 15 for IDSGP vs. M = 1024 for VSGP). In classification, however, all the methods seem to perform similar to each other and the differences between them are smaller. Again IDSGP uses here smaller number M = 3 inducing points. Increasing M in IDSGP does not improve the results.
Table I: Avg. neg. test log- ligellil values for the UCI regression datasets. The numbers in parentheses are standard errors. Best mean values are highlighted in bold face.
N
d
VSGP
SOLVE
SWSGP
IDSGP
Kin40k 32,000
-0.047 (0.003)
-0.415 (0.006)
-0.110 (0.007)
-1.461 (0.019)
Protein 36,584
9
2.848 (0.002)
2.818 (0.003)
2,835 (0.002)
(0.007)
KeggDirected
42,730
19
-1.955 (0.013)
-1.756 (0.073)
-22256 (0.012)
-2.410 (0.012)
KEGGU
51,686
26
-2.344 (0.012)
-2.531 (0.015)
-2.396 (0.006) -2.908 (0.042)
3dRoad
347,899
3.691 (0.006)
(0.000)
3.879 (0.026)
3.399 (0.009)
Song
412,276
90
3.613 (0.003)
3.608 (0.002)
3.618 (0.004)
3.637 (0.002)
Buzz
466,600
77
6.272 (0.012) 6.297 (0.009)
6.137 (0.008)
6.317 (0.055)
Housetleecticc
1,639,424
6
-1.737 (0.006)
-1.743 (0.005)
-1.711 (0.010)
-1.774 (0.004)
Avg. Ranks
3.125 (0.125)
2.475
2.850 (0.150) 1.550 (0.172)
# Inducing points
1024
1024 1024
(K=50)/ / 1024
15
Table 2: Avg. test neg. log-likelihood values for the UCI classification datasets. The numbers in parentheses are standard errors. Best mean values are highlighted in bold face.
N
d
VSGP
SOLVE
SWSGP
IDSGP
MagieGamma
15,216
10
0.308 (0.004) 0.314 (0.005) 0.371 (0.005)
0.311 (0.002)
DefaultOrCredit
24,000
30
0.000 (0.000) 0.000 (0.000)
0.000 (0.000)
0.000 (0.000)
NOMAO
27,572
174
0.113 (0.004)
0.103 (0.004)
0.134 (0.004)
0.121 (0.004)
Barketrg
36,169
51
0,206 (0.001) 0.199 (0.001) 0,304 (0.021)
0,209 (0.002)
Miniboone
104,051
50
0.151 (0.001) 0.142 (0.001)
0.180 (0.007)
0.153 (0.001)
Skin
196,046
3
0.005 (0.000) 0.005 (0.000)
0.006 (0.001)
0.003 (0.000)
Crop
260,667
174
0.003 (0.000)
0.003 (0.000)
0.002 (0.000)
0.003 (0.000)
HTSensor
743,193 11
0.003 (0.001) 0.001 (0.000)
0.030 (0.009)
0.005 (0.001)
Avg. Ranks
2.425 (0.143) 1.775 (0.158)
3.175 (0.182) 2.625 (0.155)
# Inducing points
1024 1024/ 1024
(K=50) 1024
3
In these experiments we also measure the average training time per epoch, for each method. The results corresponding the UCI regression datasets are displayed in Table The results for the UCI classification datasets are found in Appendix D.4 They look very similar to ones reported here. We observe that the fastest method in terms of training time the proposed approach. Namely, IDSGP. Nevertheless, the speed-up obtained is impaired by the overhead of having to compute the output of the DNN and update its parameters. IDSGP also results in fastest prediction times than VSGP, SOLVE or SWSGP. See Appendix D.3 and D.4 for further results showing average prediction times.
Table 3: Average training time per epoch across the 5 splits for the UCI regression datasets. The numbers in parentheses standard errors. Best mean values highlighted.
Kin40%
Protein
KegaDirected
KEGGU
3dRoed
Song
Buzz
HouscElectric
VSGP
(0.58) (1.12 932.7 (2.56) 1128.1 (3.78) 7880.9 (66 79) (42.84) 9901.0 (146.07) 32784.2 (190.18)
SOLVE
1739.3 (0.45)
2015.9 (0.66) 2357.3 (1.70)
2909.1 (1.19) 19567.1 (10.34) 23196.6 (98.35) 25769.5 (20112) 92214.9 (452.18)
SWSGP
875.7 (0.68)
1023.5 (0.35) 1220.6 (2.896(1.89) 1458.0 (5.57) 10203.4 (12.03) 12241.7 (62.01) 13371.5 (12.34) 46163.3 (427.23)
IDSGP
190.3 (0.75)
371.5 (1.25) 533.0 (1.73)
693.7 (5.77) 4070.1 (201.09) 4296,5 (25.03) 3640.4 (33.36) 16352.2 (90.15)
8

---

Under review as a conference paper at ICLR 2022
5.3 LARGE SCALE DATASETS
A last set of experiments considers two very large datasets. The first dataset is the Airlines Delay binary classification dataset, described in Hemmández-Lobato & Hernández-Lobate (016)) with N 2, 127, 068 data instances and d = 8 attributes. The second dataset is the Yellow taxi dataset, as described in Salimbeni & Deisenroth (2017) with N 1 billion data-points and d attributes. In each dataset we use a test set of I0, 000 instances chosen at random. The number of inducing points set to be equal to . = 50 in IDSGP. the other methods, we use the same number of inducing points as in the previous section. The mini- size is set to 100. Training is also performed on the same GPU as in the previous section. The ADAM learning rate is set to 0.001.
The average negative test log-like nigllikeldoood of each method is displayed in Figure 4 for each dataset. We report performance terms of the training time, in 1og scale. The results corresponding to the RMSE very similar the ones displayed here. They can be found in Appendix [D5] We observe that the proposed method IDSGP performs best on each dataset. In particular, it obtains a better performance in smaller computational time. We believe this is consequence of using smaller number of inducing points, and also because of the extra flexibility that the DNN provides for specifying the locations of the indueing points.
Yellow Taxi
Airplane Delays
5W557
2W557
V507
V507
100
Traning Sin i log28 scale
Traning line logg0
Figure 4: Negative log- likelihoooo on the test set for each method as a function of the training time in seconds, in logio scale, for the Yellow taxi and the Airline delays datasets. Best seen in color.
6 CONCLUSIONS
Gaussian processes (GPs) are flexible models for regression and classification. However, they have a cost of ((N*) per iteration with N the number of training points. Sparse approximations based on MM < N inducing points reduce such cost to ((M*). A problem, however, is that in some situations a large number of inducing points have to be used in practice, since they determine the flexibility of the resulting approximation. There is hence need reduce their training cost.
We have proposed here input dependent sparse GP (IDSGP), a method that can improve the training time and the flexibility of sparse GP approximations. IDSGP uses deep neural network (DNN) to output specific inducing points for each point at which the predictive distribution of the GP needs to be computed. The DNN also outputs the parameters of the corresponding variational approximation on the inducing values associated to the inducing points. IDSGP can be obtained under formulation that considers an implicit distribution for the input instance the DNN. Importantly, such a formulation is shown to keep intact the GP prior on the latent function values associated to the training points.
The extra flexibility provided by the DNN allows to significantly reduce the number M of induc- ing points used in IDSGP. Such a model provides similar or better results than other sparse GP approximations from the literature at a smaller training cost. IDSGP has been evaluated on several regression and binary classification problems from the UCI repository. The results obtained show that it improves the quality of the predictive distribution and reduces the training cost of sparse GP approximations. Better results are most of the times obtained in regression problems. In classification problems, however, the performances obtained are similar to those of the ofate-oftthe-ee the-art, although the training and prediction times always shorter. The scalability of IDSGP is also illustrated on massive datasets for regression and binary classification of up to 1 billion points. There, IDSGP also obtains better results than alternative sparse GP approximations a smaller training cost.

---

Under review as a conference paper at ICLR 2022
ETHICS STATEMENT
The authors acknowledge to have read and commit to adhering to the ICLR Code of Ethics We do see any direct potential negative societal impact, because this paper focuses on the development of new methodology. We believe these would be indirect through the particular application in which the proposed method is used. As one of the main advantages of GPs that they provide uncertainty estimates associated with the predictions made, we think the potential harm of these models in society could arise applications when this uncertainty estimation is critical. For example, an Al system in which the decisions made can have an influence on people's s life, such as autonomous vehicles or automatic medical diagnosis tools.
REPRODUCIBILITY STATEMENT
Most of the details needed to reproduce this paper's results are described in Section [ In Appendix B we give further details about the neural network architecture and initialization. We provide the code used to run the experiments in the supplementary material, with implementation of our method and the rest of methods of the experimental comparison. Regarding the pre- pre-proesssigg of the data, you can find all the information needed to reproduce this step in Appendix A
REFERENCES
M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghemawat, I. Goodfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefowicz, L. Kaiser, M. Kudlur, J. Levenberg, D. Mané, R. Monga, S. Moore, D. Murray, C. Olah, M. Schuster, Shlens. B. Steiner, Sutskever, K. Talwar, P. Tucker, V. Vanhoucke, V. Vasudevan, F. Viégas, O. Vinyals, P. Warden, M. Wattenberg, M. Wicke, Y. Yu, and X. Zheng. TensorFlow: Large- scale machine learning on heterogeneous systems, 2015. Software available from tensorflow.org.
M. Bauer, M. van der Wilk, and C. E. Rasmussen. Understanding probabilistic sparse Gaussian process approximations. In Advances Neural Information Processing Systems 29, pp. 1533-1541. 2016.
C.-A. Cheng and B. Boots. Variational inference for Gaussian process models with linear complexity. In Advances in Neural Information Processing Systems, pp. 5184-5194, 2017.
L. Csató and M. Opper. Sparse on-line Gaussian processes. Neural Computation, 14:641-668, 2002.
D. Dua and C. Graff. UCI machine learning repository, 2017. URL http://archive.ics
uei edu/ml
T. Evans and P. Nair. Scalable Gaussian processes with grid-sttrctured eigenfunctions (GP-GRIEF).
In International Conference on Machine Learning, pp. 1417-1426. 2018.
Y. Gal. Uncertainty in deep learning. PhD thesis, University of Cambridge, 2016.
Gardner, G. Pleiss, R. Wu, Weinberger, and A. G. Wilson. Product kernel interpolation for scalable Gaussian processes. In International Conference on Artificial Intelligence and Statistics, pp. 1407-1416, 2018
Xavier Glorot and Yoshua Bengio. Understanding the difficulty of training deep feedforward neural networks. In Proceedings of the thirteenth international conference on artificial intelligence and statistics, pp. 249-256. JMLR Workshop and Conference Proceedings. 2010.
R. B. Gramacy and D. W. Apley. Local Gaussian process approximation for large computer experi- ments. Journal of Computational and Graphical Statistics, 24:561-578, 2015.
R. Henao and O. Winther. Predictive active set selection methods for Gaussian processes. Neurocom- puting, 80:10-18, 2012.
10

---

Under review as a conference paper at ICLR 2022
J. Hensman, N. Fusi, and N. D. Lawrence. Gaussian processes for big data. In Uncertainty in Artificial Intelligence, pp. 282-290, 2013.
J. Hensman, A. Matthews, and Z. Ghahramani. Scalable variational Gaussian process classification. In International Conference on Artificial Intelligence and Statistics, pp. 351-360, 2015a.
J. Hensman, A. G. Matthews, M. Filippone, and 7. Ghahramani. MCMC for variationally sparse Gaussian processes. In Advances in Newral Information Processing Systems 28, pp. 1648-1656. 2015b.
D. Hernnndez-Lobato and J. M. Hemánnezz-Lobat.. Scalable Gaussian process classification via expectation propagation. In Artificial Intelligence and Statistics, pp. 168-176, 2016.
D. Kingma and J. Ba. ADAM: a method for stochastic optimization. In Insernational Conference on Learning Representations, pp. 1-15, 2015.
D. Kingma and M. Welling. Auto-encoding variational Bayes. In International Conference on Learning Representations, 2014.
N. Lawrence, M. Seeger, and R. Herbrich. Fast sparse Gaussian process methods: The informative vector machine. Neural Information Processing Systems, pp. 609-616, 616, 2003.
A. Naish Naish-Guzman and S. Holden. The generalized FITC approximation. Advances in neural information processing systems, 20: 20:1057-1064, 2007.
F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Pretten- hofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825 28300, 2011.
Quifionero-Camtne Candela and E. Rasmussen. A unifying view sparse approximate Gaussian process regression. Journal of Machine Learning Research, 6:1939 119999-1999, 2005.
C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning). The MIT Press, 2006.
H. Salimbeni and M. Deisenroth. Doubly stochastic variational inference for deep Gaussian processes. In Advances Neural Information Processing Systems, pp. 4588-4599, 2017.
H. Salimbeni, C.-.. Cheng, B. Boots, and M. Deisenroth. Orthogonally decoupled variational Gaussian processes. In Neural Information Processing Systems, pp. 8725-8734, 2018a
H. Salimbeni, S. Eleftheriadis, and J. Hensman. Natural gradients in practice: Non-conjugate variational inference Gaussian process models. In International Conference on Artificial Intelligence and Statistics, pp. 689 0099-697420 2018b.
M. W. Seeger, C. K. I. Williams, and N. D. Lawrence. Fast forward selection to speed up sparse Gaussian process regression. Proceedings of the Ninth International Workshop on Artificial Intelligence and Statistics, pp. 261-261, 2003.
J. Shi, M. Titsias, and A. Mnih. Sparse orthogonal variational inference for Gaussian processes. In International Conference on Artificial Intelligence and Statistics, pp. 1932-1942, 2020.
Jiaxin Shi, Mohammad Emtiyaz Khan, and Jun Zhu. Scalable training of inference networks for gaussian-process models. In International Conference on Machine Learning, pp. 5758-5768. PMLR, 2019.
R. Shu, H. H Bui, Zhao, M. Kochenderfer, and Ermon. Amortized inference regularization. In Advances in Neural Information Processing Systems, pp. 4393 4402, 2018.
E. Snelson and Z. Ghahramani. Sparse Gaussian processes using pseudo-inputs. In Advances in Neural Information Processing Systems, pp. 1257-1264, 2006.
Z. Snelson, E.d Ghahramani. Local and global sparse Gaussian process approximations. In International Conference Artificial Intelligence and Statisties, pp. 524-531, 2007.
11

---

Under review as a conference paper at ICLR 2022
Shengyang Sun, Guodong Zhang, Jiaxin Shi, and Roger Grosse. Functional variational bayesian neural networks. In International Conference on Learning Representations, 2019.
M.K. Titsias. Variational learning of inducing variables in sparse Gaussian processes. In International Conference on Artificial Intelligence and Statistics, pp. 567-574, 2009.
G.-L. Tran, D. Milios, P. Michiardi, and M. Filippone. Sparse within sparse Gaussian processes using neighbor information, 2020. 201 1.. 1.05041 arXiv stat.ML.
V. Tresp. A Bayesian committee machine. Neural Computation, 12:2719-2741, 2000.
C. Villacampa-Calye B. Zaldívar, E. C. Garrido-Merchnn, and D. Hernández-Lobato Multi-class gaussian process classification with noisy inputs. Journal of Machine Learning Research, 22:1-52, 2021.
A. G. Wilson and H. Nickisch. Kernel interpolation for scalable structured Gaussian processes (kiss-gp). In International Conference on International Conference on Machine Learning, pp- 1775-1784, 2015.
A DATASETS PRE- PRE-PROCESSING
All the datasets are publicly available. The UCI repository datasets can be downloaded from the repository (Dua & Graff (017). Yellow taxi dataset was preprocessed following (Salimbeni & Deisenroth (0077) and downloaded from htt:ps: wtwp:/wwww.........wooo/o
gov wwt.llggoosssitiiggeeeeese /tlcc///ec//ee/e /wwotttb/t/t//ooo//o///..oo..ttt
Eic-trip record- data page, where have used data records from year 2015. Similarly, the Airlines Delay dataset was preprocessed
& Hernández-Lobato/ (2016)
and was downloaded from htttps:/ccommuityy.....m wtt.y...tttttt........ dataexpo/ dat aetpo keeping only the records from January 2008 to April 2008. All datasets have been standardized using scikit-learn's 's built-int StandardScaler class (Pedregosa et al.) 2011), which removes the mean and scales to unit variance.
nt www.:oinnntcsssssseeooeeese on/
B NEURAL NETWORK ARCHITECTURE
About the choice of architecture for the DNN we have tried to keep it small in order to take more advantage of the computational gain of the amortized scheme. In particular, we used 2 hidden-layer with 50 hidden units network for the toy problems Section [5.]] network for the UCI datasets in Section 552 and a 2 layer with 25 hidden units for the large scale datasets in Section 5.3 We used Relu activation functions. Keeping the network small reduces the number of parameters to optimize making the optimization process easier. In all problems using fully-connected layers with batch normalization and no skip layers. Regarding the initialization of the weights, all were initialized using the Glorot initialization (Glorot & Bengio 2010). In experiments we did not exhaustively explore the DNN architecture. This choice architecture and initialization was based some preliminary tests done before running the experiments. This does not mean that this is the best possible configuration. We did not optimize the architecture of the neural network. In practical applications, we suggest a configuration that performs well. The main suggestion, however, is to keep the network small as the input dependence will make the model expressive enough to still get very good results.
layer with 50 hidden units
some preliminary tests in order choose
C CHOOSING THE NUMBER OF INDUCING POINTS
We have observed that our proposed method IDSGP performs well in general with fairly small number of inducing points, much smaller than the number of inducing points used in the other methods. Namely, SOLVE, VSGP and SWSGP. This is probably related to the extra flexibility of having input-dependent inducing points in IDSGP. In very large datasets we recommend using around M = 50 inducing points. In medium-size regression datasets M 15 inducing points seem enough. In medium-size binary classification datasets, however, a smaller number of inducing points is enough M = 3. We believe the reason is that binary classification problems require less complicated latent functions. We did not optimize the number of inducing points. In practical applications, we suggest to run some preliminary tests in order to choose number of inducing points that performs well.
12

---

Under review as a conference paper at ICLR 2022
C.I KL-DIVERGENCE MINIMIZATION
In this section we show that maximizing (x, u) and p(x, u particular,
effectively minimizes the KI-divergence between
/ qx,, log
ply, X, f. u)
KL(g(X, f,
dxdfdu + const.
9(x. f, u)
/ qxx,, f. u) log
p(y f p) (
dxdfdu const.
p( (ulfll)gee (ux))(
/ ((x,, f,) log
p(y p(y(f) dxdffuu + const.
q( x)
/ q(,, u) log dog
/ quu,, x) log p(u x
dxdu + const.
g(u x)
Egg Bly + ( ([(x)])]((((
+ const.,
where we have used that the posterior is equal to the joint p(x, f. y) divided by normalization constant, i.e., the marginal likelihood. Moreover, L is simply the lower bound defined in (4). There- fore, maximizing C effectively leads to the minimization of the KL-divergence between g f. u) and p(8, f. uly).
D EXTRA EXPERIMENTAL RESULTS
In this section, we include some extra results that do not fit in the main manuscript. Namely, the RMSE in the test set results and prediction times for the UCI regression datasets, and the accuracy in the test set, training and prediction times for the UCI classification datasets. In both cases, the setup the same as described in Section 5 and the results are similar that the ones obtained in terms of the negative test log likelihood and training times in that section. Finally, we include similar plots to those in Section 5.3 but in terms of the test RMSE for the Yellow Taxi dataset and in terms of the test classification error for the Airline Delays dataset.
D.1 TOY REGRESSION DATASETS
Our method looks more and more similar to the full GP as number of inducing points increases. However, with small number of inducing points, it gives similar results to those of the full GP and similar to the results obtained when more inducing points are considered, which does not happen in the other methods. This is probably due to the extra flexibility of the neural network. The figures below (Figures to 8) show the results of each method on the toy regression problem as we increase the number of inducing points M. For M 128 IDSGP gives almost the same results VSGP.
13

---

Under review as a conference paper at ICLR 2022
(number
Indacing points
(number Inducing paints point
Figure 5: Toy regression example by varying number of inducing points M2 (2. 2, 16, 32, 64, 128) with location of initial and final inducing points for an arbitrary selected point r from training sets. The mean and standard deviation of full GP prediction are shown with dashed blue and brown lines, respectively. The blue lines and the dashed red lines the mean and standard deviation of IDSGP.
14

---

Under review as a conference paper at ICLR 2022
SOLVE Orducing points inducing
SOLVE (indacigg inducing
SOLVE Orducing points
Inducing
SOLVE Onducing Inductig
Ronder
SOLVE Oedacing poinne U: 12K Inducing polntcc $28)
Figure 6: Toy regression example by varying number of inducing points M,, Me (2,4,8, 16, 32, 64, 128) with location of initial and final inducing points. The mean and standard deviation of full GP prediction are shown with dashed blue and brown lines, respectively. The blue lines and the dashed red lines are the mean and standard deviation of SOLVE.
15

---

Under review as a conference paper at ICLR 2022
VSGP Ostreber induiing points 33)
VSGP (sumber Intlcing polnta 04)
VSGP (number induding potnts 126)
Figure 7: Toy regression example by varying number of inducing points M (2,4,8, 16, 32, 64, 128) with location of initial and final inducing points. The mean and standard deviation of full GP prediction are shown with dashed blue and brown lines, respectively. The blue lines and the dashed red lines are the mean and standard deviation of VSGP.
16

---

Under review as a conference paper at ICLR 2022
painte complis
pointe
SWSGP (number desest intuding points 32, total number induding points 128)
swiee (sumber of closest indoring points 65. futal numsher of Induaing points $38)
Figure 8: Toy regression example by varying number of the neighbor inducing points Mc {2, 4, 8, 16, 32, 64, 128) and total number of inducing points M = 128, with location of initial and final inducing points. The mean and standard deviation of full GP prediction are shown with dashed blue and brown lines, respectively. The blue lines and the dashed red lines are the mean and standard deviation of SWSGP.
D.2 EXTRA RESULTS FOR THE TOY REGRESSION EXPERIMENT
Here we run the toy regression experiment of Section 5.1 using the closed-form solution approach Titsias (2009) for finding g. More precisely, this method is exactly the same as the SVGP method we compare results with, but where the approximate distribution is not optimized at all. The reason for this is that it is possible to find a closed-form solution for However and importantly, the resulting method does not allow for mini- batch training. Since SVGP* does not allow for stochastic optimization, the batch size is set equal to the number of training points (N 200). Figure 9 shows the fit obtained for an increasing number of inducing points M. The results are very similar to the ones of SVGP in Figure] 7
17

---

Under review as a conference paper at ICLR 2022
VSGD Insmber of Indudig paints
VSGP Cumbbr of Induding poiece
Intucing palnts
Indacing paints
Weatee tertiis
VSOP Ctumber helucing polnts
32
VSOP Tnumber Indacing points 64)
Figure 9: Toy regression example by varying number of inducing points M {2, 16, 32, 64, 128) with location of initial and final inducing points. The mean and standard deviation of full GP prediction are shown with dashed blue and brown lines, respectively. The blue lines and the dashed red lines are the mean and standard deviation of VSGP*
D.3 UCI REGRESSION DATASETS
Table 4: Test Root Mean Squared Error (RMSE) values for the UCI regression datasets. The numbers in parentheses standard emror.. Best mean values highlighted.
VSGP
SOLVE
SWSGP
IDSGP
Kin40k 32,000
0.198 (0,002)
0.157 (0.001)
0.215 (0.002)
0,050 (0.002)
Protein 36,584
4.161 (0.011)
4.062 (0.011)
4.133 (0.008)
3.756 (0.019)
KeggDirected
42,730
0.032 (0.001)
0.079 (0.032)
0.024 (0.000)
0.022 (0.001)
KEGGU
51,686
26
0.024 (0.000)
0.020 (0.000)
0.022 (0.000)
0.014 (0.000)
3dRoad
347,899
9.641 (0.063)
10.020 (0.095)
11.726 (0.327)
7.250 (0.069)
Song
412,276
90
8,966 (0.022)
8.925 (0.020)
9.013 (0.029)
9.068 (0.011)
Buzz
466,600
77
175.076 (15.021)
173.352 (14.957)
160.744 (13.467)
166.78 166.784 (18.040)
HouseElectric
1,639,424
0.035 (0.000)
0.034 (0.000)
0.036 (0.001)
0.032 (0.000)
Avg. Ranks
3.075 (0.126)
2,400 (0.138)
3.025 (0.170)
1,500 (0.151)
# Inducing points
1024
1024 1024
(K=50) 1024
15
18

---

Under review as a conference paper at ICLR 2022
Table 5: Average prediction time per epoch across the 5 splits for the UCI regression datasets. The numbers in parentheses are standard errors. Best mean values are highlighted.
Kin40k Protein KeggDirected KEGGU 3dRoad
Song
Buzz
HouseElectric
VSGP 0.9(0.00) 1.1(0.0)
1.4(0.01)
1.7(0.01) 11.6(0.07) 14.5(0.08) 18.0(0.08)
48.3(0.12)
SOLVE 2.4(0.00) 2.8(0.0)
3.2(0.00)
4.0(0.00) 27.0(0.04) 31.8(0.19) 37.0(0.03)
127.7(0.43)
SWSGP 1.3(0.00) 1.5(0.0)
1.8(0.01)
2.1(0.01) 14.8(0.03) 17.6(0.02) 21.3(0.10)
66.2(0.88)
IDSGP 0.3(0.00) 0.5(0.0)
0.7(0.00)
1.0(0.02) 5.7(0.26) 5.6(0.05) 8.1(0.08) 22.1(0.14)
D.4 UCI CLASSIFICATION DATASETS
Table 6: Test Accuracy values for the UCI classification datasets. The numbers in parentheses are
standard errors. Best mean values are highlighted.
N
d
VSGP
SOLVE
SWSGP
IDSGP
MagicGamma
15,216
10
0.876 (0.001) 0.877 (0.002)
0.867 (0.002)
0.877 (0.002)
DefaultOrCredit
24,000
30
1.000 (0.000) 1.000 (0.000)
1.000 (0.000)
1.000 (0.000)
NOMAO
27,572
174
0.956 (0.002) 0.960 (0.001)
0.961 (0.001) 0.955 (0.001)
BankMarn
36,169
51
0.906 (0.001) 0.907 (0.001)
0.897 (0.001)
0.905 (0.001)
Miniboone
104,051
50
0.941 (0.001) 0.945 (0.001)
0.938 (0.000)
0.937 (0.001)
Skin
196,046
3
0.999 (0.000) 0.999 (0.000)
0.999 (0.000) 0.999 (0.000)
Crop
260,667
174
0.999 (0.000) 0.999 (0.000)
0.999 (0.000) 0.999 (0.000)
HTSensor
743,193
11
0.999 (0.000) 1.000 (0.000)
0.989 (0.003)
0.999 (0.000)
Avg. Ranks
2.475 (0.127) 1.975 (0.152)
2.900 (0.187) 2.650 (0.168)
# Inducing points
1024
1024 1024 (K=50)/1024
Table 7: Average training time per epoch across the 5 splits for the UCI classification datasets. The
numbers in parentheses are standard errors. Best mean values are highlighted.
Magic
DefaultOrCredi
NOMAO
BankMarket Minihoone
Skin
Crop
HTSene
VSGP 3105(459)
4759(516)
4445(549)
6231(862)
18447(1279)
49962(9292) 37835(7063)
(1516)(17511)
SOLVE 5154(1061)
7554(1039)
6718(1088)
8949(1635) 37022(7902) 64606(13314 88819(18864)
168709(21194) 194)
SWSGP 1547(145)
2354(182)
2728(188)
3682(351)
10040347)
20283(2796)
21770(3038)
67687(5880)
IDSGP 1143(100)
1293/90)
2026(94)
2987(354)
7654(134)
15700(1988) 21378(2561)
53895(5652))
Table 8: Average prediction time per epoch across the 5 splits for the UCI classification datasets. The numbers parentheses are standard errors. Best mean values are highlighted.
MagicGamma DefauliOrrreiir NOMAO BankMarketing Miniboone
Skin
Crop
HTScer
VSGP
3.6(0.56)
4.1(0.59)
4.4(0.74)
9.5(4.39)
(7.9(1.84)
19.7(11.15) 58.7(12.22)
(41.66)))67)
SOLVE
4.0(0.85)
5.4(0.73)
3.4(0.76)
4.9(1.2))
49.2(17.83) 48.8(16.73) 86.8336.26 26) 93.5(1573))
SWSGP
3.0(0.43)
3.910.38)
4.3(0.51)
5.4(0.72)
16.4(1.8)
36.008.00) 33.9(6.25)
88.4(9.16)
IDSGP
2.5(0.24)
2,5(0.21)
3,510.39)
4.8(0,54)
13.9(0.75) 26.1(4.92) 37.5(4.6)
83.4(8.23)
D.5 LARGE SCALE DATASETS
Yallow
Airplane Delays
1050
550
Tear
30
6.20
6.30
30
Traking time 1030 scale
Trxking fine nlgg
Figure 10: (left) Test RMSE for each method as a function of the training time in seconds, in log scale, for the Yellow taxi dataset. (right) Prediction error on the test set for each method as function of the training time in seconds, in log scale, for the Airlines Delays dataset. Best seen in color
19

---

Under review as a conference paper at ICLR 2022
D.6 NEURAL NETWORK TRAINED VIA MAXIMUM LIKELIHOOD
In this subsection we show extra experiment results on the UCI datasets when using a neural network trained via maximum likelihood. The architecture of the neural network is the same as the one of the network used the proposed method IDSGP. Training is done using ADAM. The learning rate used 0.001. The mini-batch size the same for the GP-based methods. regression, we use the neural network to predict the mean and variance of the Gaussian predictive distribution. In classification, we use sigmoid activation function. The average test negative log-likelihood obtained in each problem is shown in Table and Table The results are high-lighted in bold-face when the NN performs worse than any other GP based method. The tables show that, in the case of regression problems, most of the times the neural network performs worse than the GP based methods. In the case of classification problems, the performance of the neural network is worse in those problems in which the error is higher according to Table 6 By contrast, in those problems in which the accuracy is almost equal to 100%, there are no differences or performs slightly better.
Table 9: Avg. neg. test log-likelihood values for the UCI regression datasets for the neural network. The numbers in parentheses are standard errors.
Kin40k Protein KeggDirected KEGGU 3dRoad
Song
Buzz
HouseFlectrie
0.099(0.03) 2.794(0.02) -2,407(0.14) -5.124(0.18) (.661(0.02) 3.363(0.01) 27.840(33.44) -2.1 -2.TT0(0.03)
Table 10: Avg. neg. test log-likelihood values the UCI classification datasets for the neural network.
The numbers in parentheses are standard errors.
MagicGamma CefaultOrCredit NOMAO BankMarketing Miniboone Skin
Crop
HTssnsor
0.315(0.02) 0.000(0.00) 0.119(0.01) 0.232(0.00) 0.152(0.00) 0.002(0.00) 0.002(0.00) 0,000(0.00)
20


Reviews:
[{"review_id": "zFovgyk19Hp", "reviewer": "Reviewer_rpPD", "summary": "", "questions": "", "limitations": "", "rating": 3, "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": "The main contribution of the paper is using a deep neural network (DNN) to find inducing points' locations and variational distribution parameters. In general form, it can reduce the computational cost of the original methods (SWSGP and SVGP) by training the variational parameters as the output layer of a DNN. But such a reduction in computational cost happens when it uses a small number of inducing points for each data point and small network architecture. However, there are concerns about the paper:\n\n1- Novelty and Contribution: The paper has been inspired by the SWSGP method Tran et al. (2020). SWSGP samples a mini-batch from the input data set $\\mathcal{D}$ and for each data point in the mini-batch, it finds H nearest inducing points. Instead of finding the nearest inducing points, IDSGP uses a deep Neural network and amortized variational inference Kingma and Welling (2014) and Shu et al. (2018) (available method) to find a small number of inducing points for each data point in mini-batch. The other parts of the algorithm in IDSGP are the same as the proposed algorithm in Tran et al. (2020). Indeed, using DNN to estimate variational parameters is not a novel approach. For instance, this network architecture is used in the encoder part of the conventional variational auto-encoder (VAE) (and also variational inference with the normalizing flow) where the parameters of the variational distribution (and flow parameters in normalizing flow) are estimated via DNN.  \n    \n 2-  Limitations of the Method: The limitations of the method are not properly discussed. The objective function has severely pathological properties which have not been described in the paper. For instance, consider the $n=M=1$ case, where n is the batch size, and a stationary kernel with unit variance. The optimized bound is exactly equal to a sum of GPs with one point each i.e. $\\sum_i N(y_i, | 0, 1 + \\sigma^2)$, and this does not depend on the prior at all. In fact,  as all priors are equivalent to the white noise prior. The SWSGP method does indeed have a similar problem, considering the $n=M=H=1$, but in the paper, it is explicitly described how the effective prior is modified. \n    \n3- Fall back to the Prior:  In regions of the space without training data the extrapolation is determined by the DNN and there is no guarantee the predictions fall back to the prior. This is actually apparent in Figure 2 where VSGP is guaranteed to give the prior when there is no data or inducing points, whereas IDSGP visibly extrapolates in a way that doesn't seem sensible.\n    \n4- Performance: It is not clear why the IDSGP method can provide better prediction quality than SWSGP or SVGP while it only uses a small number of inducing points. Why is the DNN is better than the nearest point sampling in SWSGP, given that they maximize the same lower bound? SWSGP uses some nearest inducing points for each data point while IDSGP uses a small set of inducing points that have been estimated for each data point. Why does the second one have better performance? Since H-inducing points in SWSGP are selected based on input points, SWSGP is also an input-dependent SGP method. There is no reason to accept the inducing point set (in IDSGP) can lead to better predictive posterior and its improvement in prediction quality is ambiguous. For instance, in Appendix D.1, when the number of nearest inducing points (H) increases, SWSGP presents better results and its uncertainty is close to the standard GP. \n   \n 5- Computational cost: The paper claims that since it uses a fewer number of inducing points, computational cost can be improved. However, it is not always true and the paper has not presented a related analysis about that. For instance, let consider the toy example in Section 5.1 and Figure 2. In this small data set, SWSGP is a little faster than IDSGP based on the computational cost discussions in both papers and experiment settings (i.e. the number of the layers and hidden units in DNN, number of inducing points in IDSGP, and the value of H and M in SWSGP). Besides, for a non-diagonal variational covariance matrix S, the output layer of DNN has $2M + \\frac{M^2+M}{2}$ variables (M for inducing points, M for mean and $\\frac{M^2+M}{2}$ for Cholesky matrix). When it needs more inducing points for each input point, the computational cost of DNN increases because of the dimension of the output layer and also network parameters (number of layers and hidden units).      \n\n", "summary_of_the_paper": "Paper considers the sparse approximation Gaussian processes (GPs) which tries to scale GPs for large data sets. The authors in this paper propose a novel method (IDSGP) to improve the computational cost of sparse GPs. With the proposed method, the number of the inducing points is reduced drastically and only a small number of inducing points are estimated for each input point. The inducing points locations, as well as the parameters of the variational posterior distribution, are the outputs of a Deep neural network (DNN). \n\nFirst, they sample a mini-batch of the training data. For each data point in the mini-batch, their solution produces a small set of inducing points (and associated variational approximation) that are specific to each input data point.\n\nIt is an extension for Sparse within a Sparse GP (SWSGP) method Tran et al. (2020). For M inducing points, the SWSGP operates on a subset of H inducing points for each input point, with $H \\ll M$, while maintaining a sparse approximation with M inducing variables. Here, instead of H nearest inducing points, the model estimates M inducing points and variational distribution parameters in a DNN. The claim is the proposed IDSGP method can provide competitive predictions while it improves the time complexity of the state-of-the-art sparse GPs baselines.", "summary_of_the_review": "In general, the paper is well-written. However, the proposed model has not been explained appropriately in Section 3. \n\nSection 3.2 is the main contribution of the paper which proposes the computational graph that estimates inducing points and variational parameters. Here it is necessary that the authors address the major issues of their method: why maximizing the lower bond (Equation (4)) in the proposed DNN leads to better results than the available baselines (e.g. SWSGP and SVGP which maximize the same objective function)?  Why does a small number of inducing points provide a better posterior than other baselines which find more inducing points (for instance SWSGP provides more inducing points but uses only H nearest inducing points for each input point)? Indeed, it is not clear how the batch permutation-invariance is achieved.  Something like a sum over the batch dimension is clearly not a good idea because the encoder is sensitive to the ordering of the batch. However, the paper does not mention how this issue is considered.\n\nBesides, here are some omissions of related work that are not much discussed but could be mentioned, e.g. \nhttps://arxiv.org/pdf/1303.0383.pdf, \nhttps://mlg.eng.cam.ac.uk/zoubin/papers/aistats07localGP.pdf, \nhttps://www.dbs.ifi.lmu.de/~tresp/papers/bcm6.pdf,\nhttp://proceedings.mlr.press/v84/salimbeni18a.html.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "d_F9sbd9wgc", "reviewer": "Reviewer_spvN", "summary": "", "questions": "", "limitations": "", "rating": 3, "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": "1) Strengths\n - a) The manuscript is rather well written and simple to understand.\n - b) The building blocks are established and well-understood.\n - c) The evaluation is done on relevant public datasets.\n - d) Code is available.\n\n2)Concerns\n - a) It remains open which properties of the GP framework used in the motivation for the work remain intact after all.\n - b) The empirical evaluation is lacking at least three simple baselines.\n     i) A neural network predicting the class probability rather than the variational parameters with sigmoid activation in the last layer. Along the same line, a regression network predicting mean and standard deviation (or log standard deviation). Both networks should be trained with maximum likelihood.\n     ii) A sparse GP using FITC approximation.\n     iii) A neural network predicting the inducing points rather than the variational parameters which is then used in combination with FITC or similar.\n\n3) Remarks\n - a) The Abstract states that \"for some learning tasks a large number of inducing points may be required\" and that the proposed method is going to address this point. However, the conclusion on the \"drastically reduced\" number of inducing points is a little unfair as one could argue that all the complexity is shifted into the network for inducing point prediction. One should rather worry about number of parameters instead.\n - b) I have the impression that a direct neural network prediction of the class probability or the mean/std (possible mean/log(std)) would do a similar job as IDSGP and that the training would be much faster.", "summary_of_the_paper": "The manuscript proposes a variant of a sparse variational approximation in order to render Gaussian process inference scalable. The central idea of the work is to use variationa inference and have a neural network predict the parameters of the variational distribution directly from the input data rather than optimising them using ELBO. Experiments on several substantial regression and classification datasets conclude the paper. The available code was not checked for this review.", "summary_of_the_review": "The proposed algorithm is certainly an option that has not been explored before but the practical value is at least questionnable due to the missing baseline experiments.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "ZjtCQTF--b4", "reviewer": "Reviewer_QKXW", "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": "The paper has made a trade-off between training speed (pros) and throwing away the readily available gradient from the ELBO (cons).\n\nStrengths:\n* Optimizing the ELBO for Sparse GP using deep neural network is interesting for improving the running time (computational cost)\n\nWeaknesses:\n* Training a deep neural network may require a large amount of data as opposed to the standard sparse GP.\n\n* It could be a drawback from the design in high-dimensional input setting. This can be inefficient to learn a mapping from one input (in high-dimensional) to multiple output (again in high-dimensional) of the inducing points. Learning such mapping may require a significantly large amount of training data.\n\n* It is unclear about how to make use of the ELBO gradient in learning the neural network. On the one hand, existing variational inference techniques for sparse GP make use of such gradient when taking partial derivative with respect to the variational parameters. On the other hand, with deep neural network, one can utilize the automatic differentiation to perform back-propagation to learn the network (as a black-box). Note that optimizing the ELBO as a black-box (with deep neural network) using automatic differentiation may be inefficient if we have the gradient signal readily available from the ELBO.\n\nNovelty:\n* The paper follows Tran et al 2020 to derive the lower bound on the log-marginal likelihood.\n\n* Note that existing works has considered using deep neural network to optimize the ELBO, (but not in the setting of sparse GP), such as [1]\n\nWriting and presentation:\n* The paper is well written and nicely presented.\n* The related work is great.\n\nSuggestions:\n* One of the key advantages of GP against other supervised learning approaches is the uncertainty quantification. Therefore, additionally to presenting the classification and regression performance on UCI datasets, it is recommended to study/analyse/compare the uncertainty estimation. Without such uncertainty estimation, people can simply perform standard deep learning method to learn these classification and regression.\n\n\nMinor points:\n* in the last paragraph of the introduction, the paper claims to have different sets of inducing points associated to each input location. However, given the network, it seems that each input will only produce a single set of inducing point, not a set of including points?\n\n* what is the dimension of Cholesky factor of the covariance matrix S=LL^T and the number of inducing points per data point Z(x_i)\n\n* In Fig 1, there is a typo at the first column, last row which should be h_P^(1), instead of h_P^(L)\n\n\n[1] Blundell, C., Cornebise, J., Kavukcuoglu, K., & Wierstra, D. (2015, June). Weight uncertainty in neural network. In International Conference on Machine Learning (pp. 1613-1622). PMLR.\n\n", "summary_of_the_paper": "This paper proposes to combine sparse GPs with a neural network architecture to compute the\ninducing points locations associated to each input point. In particular, the paper employs a neural network to carry out amortized VI to compute the parameters of the approximate variational distribution\nq modeling the posterior distribution associated to the values of the outputted inducing points. The inducing points are given by a mapping from the inputs provided by a neural network. The paper demonstrates improvement on training and prediction time while the proposed method is able to perform similar or better than the standard sparse GP approaches.\n", "summary_of_the_review": "The paper proposes to improve the training time and flexibility of sparse GP approximation using deep neural network. While it makes sense to optimize the ELBO using deep neural network, the proposed framework may have certain limitation, such as learning the mapping from one (high-dim) input to multiple (high-dim) output can be inefficient, especially when we have limited data. Another suggestion is to demonstrate the uncertainty estimation of the proposed approach.\n\nThe paper is well written and presented. The related work section is well covered.\n\nThe paper is currently in the border-line.\n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "IOic1amo1bF", "reviewer": "Reviewer_m21i", "summary": "", "questions": "", "limitations": "", "rating": 5, "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": "\n## Significant Comments and Questions\n*  Suppose I am interested in making a prediction of the joint distribution at a pair of inputs $x,x'$. How would I do this with you method? If I understand correctly, I would feed both $x$ and $x'$ into a neural net, which would output two distinct sets of inducing inputs, as well as distributions over them. Can you clarify how these can be combined to approximate, for example, the posterior covariance between $x$ and $x'$?\n* Relatedly, it seems to me as though the posterior predictive distribution $f$ changes between train and test time. In particular, the posterior distribution depends on $p(\\tilde{x})$ which differs at train and test time. Do you foresee this leading to any issues?\n* Matthews et al 2015 established that maximizing the lower bound in Titsias, 2009 is equivalent to minimizing a KL-divergence between stochastic processes. This is essential to justifying maximization of the lower bound with respect to $Z$ as a form of variational inference. Given the augmented model, it is no longer clear to me that a similar argument applies. I am curious to hear if the authors have any thoughts as to whether it does or does not. If not, does the method risk over-fitting by optimizing $Z$? Relatedly, can the quality of inference (as measured by KL-divergence to the posterior) be shown to be monotonically increasing in $Z$?\n* The authors should consider adding the closed form solution for approximate Gaussian process regression in experiments (Titsias, 2009). While more computationally expensive per iteration that the method of Hensman et al 2013, using the analytic form of the optimal posterior often improves inference quality signficicantly, and convergence of the bound proposed in Hensman et al 2013 can be slow.\n\n## Minor Comments and Questions\n* The authors repeatedly refer to the \"training cost\" when they mean the \"cost per iteration of training\". The training cost should incorporate how long it takes for the training of the method to converge.\n* At several points you state that at least several hundred inducing points are needed for complicated problems. It seems as though several thousands of inducing points might be a more appropriate statement. Generally using several hundred inducing points would not represent a major computational obstacle. The lower bounds on the number of inducing points needed  for the KL-divergence between approximate posterior and prior to not be large in Burt et al, 2020 might also be of interest.\n* On page two, $k(x,x') = \\mathbb{E}[f(x)f(x')]$ is only correct if the mean is $0$. While you state that you generally assume this, it should be more clear this is needed in this definition.\n* The last paragraph in section 2.1 is copied almost word for word from an earlier paragraph.\n* In experiments, when results are compared, are model parameters selected via ELBO maximization (approximate maximum marginal likelihood)?\n\n## References\n``On Sparse variational methods and the Kullback-Leibler divergence between stochastic processes'' Matthews, Hensman, Turner, Ghahramani. 2015.\n\n``Convergence of Sparse Variational Inference in Gaussian Processes\nRegression'', Burt, van der Wilk, Rasmussen.", "summary_of_the_paper": "The paper proposes a variant of the method in Tran, 2021 for approximate Gaussian process regression. The major changes with respect to the classic Titsias, 2009 paper on variational Gaussian process regression are:\n* The inducing inputs are data-dependent, in the sense that $Z=Z(x_i)$. \n* The inducing inputs, posterior mean, and posterior covariance are learned via amortized variational inference through learning a neural network.\n\nThe authors show some promising empirical evidence on commonly used regression and classification benchmarks.\n", "summary_of_the_review": "## Summary of Review\nThe method is interesting and the empirical results seem promising. I do have some concerns that the posterior distribution defined over $f$ is not consistent. Perhaps this is not essential if only marginal predictions are needed but it does seem to be a significant drawback in cases when covariances or samples from the approximate posterior are needed.  This also makes me slightly concerned that the objective function is not quite as well founded as that of Titsias, 2009 (see the discussion in Matthews et al 2015). I think some of these critiques apply equally to Tran et al, and I am curious to hear the authors thoughts on these limitations. \n\nThe writing is generally reasonably clear, though aspects could be improved. While I have rated the paper as marginally below the acceptance threshold, I think it is quite close to the threshold.\n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "FP74Abuh3kX", "reviewer": "Reviewer_oyvg", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I liked the paper a lot once I understood what was going on.  But even with some familiarity with variational GPs, it took me a while to process this.  Most of my comments are individually minor, but they share a common theme that they are not easy to process without either reading out of order or knowing something about variational GPs already.  I think the idea here is quite interesting, and hope thinking about these points will make the paper more accessible!\n\nDetailed comments:\n\n- Top of p2: \"This reduces the training cost to O(M^3)\" -- perhaps say \"O(M^3) per iteration\"?\n\n- The kernel is introduced as k(x,x') = E[f(x) f(x')], but perhaps it would be better to write E[(f(x)-m(x)) (f(x')-m(x'))] since the mean function is only \"often set to zero\" at this point?\n\n- On page 2, we have \"a popular covariance function k(.,.) is the squared exponential\", but a Matern 3/2 kernel is used later in the experiments.  Perhaps it would be better to replace \"squared exponential\" with \"Matern family\" or the like?\n\n- \"In Titsias (2009), they do optimize q(u) in closed-form\" -> \"In Titsias (2009), they optimize q(u) in closed form\"\n\n- \"The ELBO can be expressed as a sum...\" -- there is no prior explicit reference to the evidence lower bound, though of course the bound itself has already been introduced.  It's probably worth a word or two here.\n\n- \"To achieve this, we consider a meta-point xtilde\" -- I still don't completely understand the role of the meta-points (or what is meta about them)!  Some intuition would be useful.  Also, it is unfortunate that the meta-points are quickly given the same notation as the data points.  It seems like they end up coinciding in the algorithm as implemented, but are conceptually distinct from the data points?\n\n- In Algorithm 1, I think standard fonts (vs the math spacing) would be better for ELBO, KL_div, and Log_lk.\n\n- In related work: \"it is introduced a new interpretation\" -> \"the authors introduce a new interpretation\"; \"that allow to consider\" -> \"that allow one to consider\"; \"it is also described a mechanism\" -> \"a mechanism is also described\".\n\n- As a general comment: the approximation for the pointwise posterior variance is clear, and this is the most important part of the uncertainty in many cases.  But it's unclear to me how one would compute the approximate posterior covariance between two points, which is also sometimes useful.", "summary_of_the_paper": "The paper describes a method that combines variational Gaussian processes with a neural network that produces inducing points and associated variational parameters depending on the evaluation location.  The approach reduces training time compared to other sparse GP approaches, and has excellent performance in regression experiments (and good performance in classification experiments).  The method also does a good job at respecting the (exact) GP prior structure on the latent functions.", "summary_of_the_review": "Though some points of the presentation could be improved, I think this is an interesting idea illustrated by good experience.  I would like to see it appear in the conference.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}]

Author rebuttal:

Meta-review: Three reviewers recommend Reject. Two reviewers recommend Accept although do not champion the paper. I believe the paper develops an interesting idea to better estimate the location of the inducing inputs in sparse GP models. However, I still think the paper would benefit from another careful revision and therefore I do not recommend Acceptance at this stage. I agree with reviewers that 1) currently the method is unable to estimate the covariance between two data points. This is important in applications of GPs for uncertainty quantification such as Bayesian optimisation. For example, including a BayesOpt example would clearly strengthen the paper. 2) empirical evaluation lacks simple baselines, e.g. Titsias (2009). The authors claim that Titsias (2009) does not scale and that's why they don't care for it. Even if this is true, including an example that helps to better compare against this method at a different scale might strengthen the model proposed here.
Decision: Reject