New Bounds for Sparse Variational Gaussian Processes
Authors: Michalis Titsias
Venue: ICML 2025 spotlightposter (2025)
Keywords: Sparse variational Gaussian process
new collapsed bound
Primary area: probabilistic_methods->gaussian_processes
TL;DR: It presents new collapsed and uncollapsed  bounds for sparse variational Gaussian processes using inducing points.

Abstract: Sparse variational Gaussian processes (GPs) construct tractable posterior approximations to GP  models. At the core of these methods is the assumption that the true posterior distribution over training function values ${\bf f}$ and inducing variables ${\bf u}$ is approximated by a variational distribution that incorporates the conditional GP prior $p({\bf f} | {\bf u})$ in its factorization. While  this assumption is considered as fundamental, we show that for model training we can relax it through the use of a more general variational distribution $q({\bf f} | {\bf u} )$ that depends on $N$ extra parameters, where  $N$ is the number of training examples. In GP regression, we can analytically optimize the evidence lower bound over the extra parameters and express a tractable collapsed bound that is tighter than the previous bound. The new bound is also amenable to stochastic optimization and its implementation requires minor modifications to existing sparse GP code. Further, we also describe extensions to non-Gaussian likelihoods. On several datasets we demonstrate that our method can reduce  bias when learning the hyperparameters and can lead to better predictive performance.

New Bounds for Sparse Variational Gaussian Processes
Michalis K. Titsias'
Abstract
time cost and O(N2) storage where N is the number of training examples. This has originated a large body of research on scalable or sparse GP methods expanded in sev- eral decades; see e.g.. Chapter 8 in Rasmussen & Williams (2006) for an early review and Heaton et al. (2018); Liu et al. (2020); Leibfried et al. (2022) for recent treatments. An important class of methods bases an approximation on a small set of M < N inducing points (Csato & Opper, 2002; Lawrence et al., 2002; Seeger et al., 2003; Snelson & Ghahramani, 2006;
Sparse variational Gaussian processes (GPs) con- struct tractable posterior approximations to GP models. At the core of these methods is the as- sumption that the true posterior distribution over training function values f and inducing variables u is approximated by a variational distribution that incorporates the conditional GP prior pf(lu) its factorization. While this assumption is considered as fundamental, we show that for model training relax it through the use of more general variational distribution (f) that depends on N extra parameters, where N is the number of train- ing examples. In GP regression, we can analyti- cally optimize the evidence lower bound over the extra parameters and express a tractable collapsed bound that is tighter than the previous bound. The bound is also amenable stochastic opti- mization and its implementation requires minor modifications to existing sparse GP code. Further, we also describe extensions to non-Gaussian like- lihoods. On several datasets we demonstrate that our method can reduce bias when learning the hyperparameters and can lead to better predictive performance.
SuliooneroCCandel & Rasmussen, 2005; Banerjee et al., 2008; Finley et al., 2009; Titsias, 2009; Hensman et al., 2013; Bui al., 2017: Burt et al., 2020) that reduce the time complexity to O2) and the storage to O(N
Among inducing point methods, the sparse variational Gaus- sian process (SVGP), introduced for standard regression (Titsias, 2009), applies variational inference to obtain a pos- terior approximation and selects hyperparameters and induc- ing points by maximizing an evidence lower bound. Unlike the prior approximation framework Rasmussen, 2005), SVGP leaves the GP prior unchanged and instead it reduces the cost to ((N by imposing a special structure on the variational distribution. This frame- work has been extended stochastic gradient optimization (Hensman et al., 2013) and non-Gaussian likelihoods (Chai, 2012; Hensman et al., 2015; Lloyd al., 2015; Dezfouli & Bonilla, 2015; Sheth et al., 2015). Also, it has been ex- plained as KL minimization between stochastic processes (de G. Matthews et al., 2016).
Caliionero-Candele &
1. Introduction
Gaussian processes (GPs) are nonparametric models for learning functions using Bayesian learning. Thanks to their flexibility and ability quantify uncertainty, GPs have found many applications in machine learning (Rasmussen & Williams, 2006), spatial modeling (Cressie, 1993), computer experiments (O' Oogaan 1978; Gramacy, 2020), Bayesian optimization (Jones et al., 1998; Garnett, 2023), robotics and control (Deisenroth & Rasmussen, 2011), 1), unsupervised learning (Lawrence, 2005) and others.
An important aspect of the SVGP method is that it uses a special form for the variational distribution. It approximates the exact posterior distribution p(f, uly) over the training function values and the inducing variables (see Section 2 for precise definitions) by variational distribution of the form qf, u) = ((f),)g(u), where q(u) is some optimiz- able distribution over the inducing variables, while (fflu) is the conditional GP prior. This special form of the varia- tional approximation seems be fundamental, and it has been applied also to more complex GP models, such as those with multiple outputs (Álvarez et al., 2010; Nguyen & Bonilla, 2014; Yousef et al., 2019), uncertain inputs (Titsias & Lawrence, 2010; Damianou et al., 2016) and multiple layers (Damianou & Lawrence, 2013; Salimbeni & Deisen- roth, 2017). However, an open question regarding the SVGP
Despite the numerous applications, GPs suffer from (
'Google DeepMind. Correspondence to: Michalis K. Titsias < www..e/ooogggo.......
Proceedings of the 42nd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

New Bounds for Sparse Variational Gaussian Processes
framework is whether this particular form of variational dis- tribution is really necessary obtain scalable computations. The answer we give in this paper is that "it is not", since at least for training GP model it can be relaxed.
To learn the hyperparameters (0, we can maximize the log marginal likelihood which is analytically available,
logp(y) (an00--pssyppii p(ylDp(()d (yyo, Kr + oo)) (2)
To this end, we derive new variational bounds for train- ing sparse GP regression models by replacing p(f(u) in the variational distribution with more general conditional distribution (f(u). This ((('u) depends on N additional pa- rameters (on top of the parameters of p(f(u)), i.e., as many as the training examples, and it is constructed to enable bet- ter covariance approximation of the underlying true factor pff|u, y). We show how to analytically optimize over the N parameters and obtain a better posterior approximation together with tighter collapsed evidence lower bound. The new bound is also amenable to stochastic gradient optimiza- tion, and its simple form suggests that it can be implemented with minor modifications to existing sparse GP code. We also describe extensions of the method to non-Gausse Gan-Gaassinn like- lihoods. Furthermore, we point out the concurrent work of Bui et al. (2025) who derived similar sparse GP approxima- tions and variational training objectives by using the same form for the (((lu) distribution.
After training we can perform predictions at test inputs X. by first computing the posterior over the corresponding test function values f.:
p(L - f
p(t. (ff()))
(3)
(f., (r. 'f Kr.s. Kr,, -K...)
and then writing the predictive density as p(y,ly)
ʃ [y,,,,, o² (p((([y)dd) which is the same as the above Gaussian but with I added to the covariance.
While the log marginal likelihood and predictive density have closed-form expressions, they require the inversion of Kr - I which costs O(N3) and is prohibitive for large datasets. Next we review methods using inducing points and particularly the variational approach (TTitsiss, 2009) that our method in Section 3 improves upon.
The remainder of the paper is as follows. Section 2 pro- vides an overview of GPs and the variational approach to sparse GPs using inducing points. Section 3 derives the new evidence lower bounds for training. Section 4 discusses connections with previous works. Section 5 presents experi- ments using several datasets showing that the new bounds can reduce underfitting bias and can lead to better predictive performance. Section 6 concludes with discussion and suggestions for future work.
2.1. Sparse Variational Gaussian Process (SVGP)
The idea of inducing points is base GP approximation on smaller set of M << N function values; see e.g.. Csato & Opper (2002); Seeger et al. (2003); Snelson & Ghahramani (2006); Sariionero-Candell & Rasmussen (2005). Snelson & Ghahramani (2006) introduced pseudo inputs by instanti- ating extra GP function values u = [f(zz)) ] evaluated at locations Z = {zm J that can be optimized freely with gradient- based methods. However, the GP prior modifi- cation procedure (Quininero-Candela & Rasmussen, 2005; Snelson & Ghahramani, 2006) does not result in rigorous approximation to the GP model. An alternative variational inference method (Titsias, 2009), next referred to as SVGP!, does not modify the GP prior but instead it augments the model with extra function values u:
2. Background
GP is a distribution over functions specified by a mean function m(x) and a covariance or kernel function k( x'), where the kernel function is parametrized by 0. By assuming that m(x)= 0 we denote GP draw as
f(x)
p f. (() p(f]u)= NV Kr u, Kir Km Kur)
augmented joint (4)
For a finite set of inputs X { the distribution over the function values f = [ (stored as N 1 vector with fn := f(x,)) is the multivariate Gaussian p(f) = (f]o, Kn) where the N x N covariance matrix K has entries [K]jy = (z))
cond. GP
(5)
p(u)= -u K.n)
inducing GP prior (6)
where Kmn is the M x M covariance matrix on the induc- ing inputs Z, Kr, is the N x M cross covariance between points in X and Z, while K,t K SVGP approximates the exact posterior pff. uly) by a variational distribution q(f, u) through the minimization of [L[g(f, (f]] y.y)].
We consider standard GP regression where we are given set of training inputs X and corresponding noisy outputs = 1 J where yn E R. Conditionally on the latent values f. these outputs follow a factorized Gaussian likelihood, p(ylt)- = IIA (y, If, o²) = Nyfff, o²1). The joint distribution over outputs y and latent values f is
Another common name for this method is Variational Free Energy (VFE); see Bui et al. (2017); Bauer al. (2016); Liu al. (2020).
(fy)f)p(f= = o Kr).
(1)
2

---

New Bounds for Sparse Variational Gaussian Processes
A critical assumption is the following choice for y:
where p(f, = J p(t, (fffff) comes from the GP consistency. For completeness, in Appendix A we include further details about SVGP such as a derivation of the col- lapsed bound and the Gaussian form of the optimal (u).
(f,, u) - (())
where (u) is an optimizable M-dimensional variational distribution, while p(flu) is the same conditional GP prior from Equation (5) that appears in the joint in (4). The KL minimization is expressed as the maximization of an evi- dence lower bound (ELBO) on the log marginal likelihood,
We conclude this review of SVGP for regression with a couple of remarks that will be useful next.
Remark 2.1. The approximation becomes exact when K Qr and the collapsed bound matches the log marginal like- lihood in (2). However, to obtain good approximations we may need sufficiently large number of inducing points (Burt et al., 2020). Otherwise the bound will cause underfitting. For instance, as studied by Bauer al. (2016) and Titsias (2009) the SVGP bound tends to overestimate the noise variance 0².
/ (( log p(y p(yI))(F dfdu p(f
logp(y) 2
f a(u) log exp(ʃ p(f(u) log (ddr du.
(())
If we optimize over g(u) and obtain the optimal choice (u) x exp {J p(fu) log df) p(u), then we can sub- stitute this g" (u) in the last line above and express the called collapsed bound, having the general form
Remark 2.2. SVGP approximates [.f y) in the exact posterior in (10) by the conditional GP (ff in the variational posterior in (7), while g(u) is treated opti- mally by KL minimization. If p(f(u,) = p(f]u) then [L ()ll uly)] = 0 and the approximation becomes exact, meaning (ff.ly) = p(f.. [y) for any f. = fx.)
feep exp
{ p(f(u) log logr(sif)ar) p(u)du
(F
which for the standard GP regression model takes the form
3. Proposed Method: Tighter Bounds
F = log.Nyyl0 Qr + 1
20 tr (Kn- Qi).
(8)
Remark 2.1 suggests that it would be useful to tighten the collapsed bound in order to reduce underfitting bias and match better exact GP training. Remark 2.2 suggests that one way to tighten the bound is to replace p(f[u), in the variational approximation in (7), with another distribution that can better approximate p(f Next we develop a method that does this while keeping the cost unchanged.
DTC log Lik
trace term
where Oir = K Kur is the M-rank Nystróm matrix. The first term in the bound is the deterministic training condi- tional (DTC) log likelihood (Seeger et al., 2003; Quiñonero- Candela & Rasmussen, 2005) while the second a regular- ization term which, since tr(Kn Qir) 2 0, promotes Qr to stay close to Kr. The inducing points Z can be learned as variational parameters by maximizing the bound jointly with the hyperparameters (0, o²), which requires M2) operations per optimization step. Hensman et al. (2013) fur- ther reduced the operations to O 3 per optimization step by applying stochastic minibatch training for maximizing the uncollapsed version of the bound; see Section 3.2.
Let write the exact form of p(f(u, y). By noting that this quantity is the exact posterior over f in a GP regression model with joint ply [f) u) we conclude that this is
1 57 1
p(f y) = =N flm(y, u), (Kn
where m(y, u) E[fuu] Kiy to -y- R[ with E[fu] = Km, and Kit = Kr Qu. Note that under this notation, p(fup)=(()))))) = (f Km). We will construct a new (f(u)) that keeps the same mean E[[[u] as p(flu) but it replaces Kr with closer approximation to the covariance (K + I)) of p(flu. y). We first write this matrix as
To obtain the form of the GP posterior over any test function values f, we can first write the exact form
p(f.ly)= = -f - p(t./ff.) u)p(f, uly)dduu
(9)
where p(f, u) is the conditional GP of f, given training function values f and inducing values u while p(f,ulyy is the posterior over (f, u) written also as
1 o2
(K₁ + 0² 1 = K (I
K
K
(12)
p(t, uly) = p(flu, y)l(()y).
(10)
Then we approximate the inverse (I + Kn) by di- agonal matrix V = diag("1,..., UN) of N variational pa- rameters > 0. In other words, in the initial q(t. u) =
The SVGP method approximates p(f, uly) by q(f, u) and therefore by plugging in this q into (9) we obtain
()) we will replace uf by
g(C.ly)= =
p(E. [f. (((((())(())((e)))) / (u)dd,
(f(u) = Kr, u, (Kr Qx) V(Kn +m) (13)
(11)
3

---

New Bounds for Sparse Variational Gaussian Processes
The ELBO now is written as
(Kr - Qrr) in Equation (12). Also note that in the ordering of positive definite matrices it holds V S I, from which it follows that q(flu) has smaller covariance than p(f u) and more accurately approximates the covariance of p( u, y). The latter as implied by Equation (12), has also smaller covariance than p((u).
ʃ (flu))( log p(yl)p((u))p()))
dfdu
((
GeloX18g e log ply p(u)
[Llg(f du
(u)
and the challenge is to see whether (L(g(f()) [ and Byg(r) p(y(f) are computable in O(N time. We have the following results (proofs are in Appendix B).
3.1. Predictions
To perform predictions we will be using the same pre- dictive posterior from Equation (I1), i.., (f,, = J pf. (u))(u)dd where the optimal q"(u) (see Ap- pendix A) is exactly the same as in the standard SVGP method. The alternative expression (and strietly speaking more appropriate since our variational approximation is (())) is given by
Lemma 3.1. (l_gg(f))]ll((iiv)
= _ - log 0.- -1).
Lemma 3.2. Let us denote the diagonal elements of Kr Qi as has qu for ,,...,NN. N. Then
E Bygm
N vz (kis qa). (14)
= =log Km. kum o2)
2
Thigh cost (f. (f.ly)= / f (17)
By combining the two lemmas the full bound is written as
But this is expensive since has cost O(N3). The reason is that f p(f, [f, (ff does not simplify anymore since (fflu) not the conditional GP, which means that p(f. u) and (ffuu) are not consistent under the GP prior. Never- theless, q(f,, ly) and Dhigh (f, y) have exactly the same mean, since ((flu) and p(flu) have the same mean, but the tractable will give higher variances than qhigh_oost.
f (((u) log
Ny K u, o pp) du
q(u)
N
kis - 4si
1).
Ci
1
log 0s -1 -
(15)
σ²
1=1
Proposition 3.3. Maximizing the bound in (15) with respect to q(u) and each v; results in the optimal settings q" (u) xc
3.2. Stochastic Minibatch Training
N KyyKKa u, p'u) and vi = (1+ ku=gu
By
The initial SVGP method (Titsias, 2009) does the training in a batch mode where all data are used in each optimiza- tion step. Stochastic optimization using minibatches was proposed by Hensman et al. (2013). Here, we apply our new approximation to this stochastic method.
substituting these values back to (15) we obtain
1 N
hist qis
Fnew log o
Eioo log ( [1+
(16)
We start from Equation (15), and substitute only the optimal values for each v without using the optimal setting for g(u). This results in the uncollapsed bound
The first term is the DTC log likelihood as in the original bound in (8), but the regularization term makes the bound tighter, i.e., log p(y) 2 Fnew 2 F, due to the inequality log(a + 1) S a. Also since log(a + 1) < a for all > O, if Kr * Qir (so there is at least one hüs qa > 0), then F > F. This means that Fnew is strictly better than F unless both bounds match exactly the log marginal likelihood.
Betwa (log,
NV MoggNii (uilk (k Kmi oz)
1 log [
kis
9sa
}
(L(g())(lP(u)),
(18)
2
o²
where ku is the 1 x M vector of all kernel values between the training input x, and the inducing inputs Z, while the expectation under quu) in the first line is analytic; see Hens- man et al. (2013). The above bound strictly better than the previous uncollapsed bound in Hensman et al. (2013), since 2 (R& ga) - - log (1 kurg" Based on the above we can apply stochastic gradient methods to op- timize qu) and the hyperparameters by subsampling data minibatches to deal with the sum over the N training points,
Clearly, Frew has M2) cost and its implementation requires a minor modification to the initial bound. The opti- mal q² (u) is the same as in the initial SVGP method, while an interpretation of the optimal ey values is the following. Remark 3.4. The diagonal matrix V (with the optimal of values in its diagonal) is the inverse obtained after zeroing out the off dii-oiagonal elements of I + 1 (Krr - Qx), i.c.,
y diag(I + 1 (Kn Qm))- which approximates (I +
4

---

New Bounds for Sparse Variational Gaussian Processes
i.e., at each iteration we use the stochastic ELBO:
where we used the fact that (L [(f[u)] is obtained from Lemma 3.1. The above bound is not computationally efficient since the marginal g(fi) = ((()))g(u))) has O(N9) cost. This is because the marginalization q(s u) = ((flu)df= cannot be trivially expressed, due to the complex structure of the covariance (Kr QQr) V(Kr Qr) à in q((')), To overcome this, we will use a simplified version of qfflu), in which we choose spherical V vl with v > 0. Then, things become tractable.
N
MoogN(yiilk,na u, o²]]
B
iEB
+ kis - qii
=mm) }} (L((()))))))) (19)
log
x²
where B denotes minibatch.
The most common parametrization of q(u) is q(u) =
N ()) where the mean vector m and covariance matrix are variational parameters. Another popular parametrization, for instance used as the default in GPflow (de G. Matthews et al., 2017), is the whitened parametriza- tion that we consider in experiments. For any choice of q(u), the new bound is always tighter than its correspond- ing previous uncollapsed bound and requires minor mod- ifications to existing implementations, i.e., to replace the
Proposition 3.5. Let q( = (fu u, ((Kr Q())
for > 0. Then (22) is computed in O time
N Σ Egcssllogp(ullA) =
N (v 2
-loo v 1)
- -KL[(())lp(u)),
(23)
previous term 20 (his qa) with 1 log 1 + ku1
2
where the marginal is q(fi) = ( k ky m, ukk
qx) + ky.k un "SK. kuf.).
3.3. Mon-Gaussian Likelihoods
The parameter v multiplies the term kn - qu inside the variance of g(fx), and it also appears in the regularization N (v - log u - 1). If = 1 the bound in (23) re- duces to (20), while by optimizing over v it can become a tighter bound. The optimization of v is done jointly with the remaining parameters m, S, Z, 0 using gradient- based meth- ods. Stochastic gradients can also be used by subsampling minibatches to deal with the sum Σ E (ogral)) and reduce the complexity to O(M)).
Consider a factorized likelihood p(ylr) H1 fA)) where p(v.] ff) is non- Gaussian, e.g., Bernoulli for binary outputs or Poisson for counts. In this non-conjugate set- ting the sparse variational GP approximation imposes the same form for the variational distribution, i.e., q(f, u) = (()) where p(flu) is the conditional GP prior. As shown in several works (Chai, 2012; Hensman et al., 2015; Lloyd et al., 2015; Dezfouli & Bonilla, 2015; Sheth al., 2015), this leads to the bound
term - 2
Li=1
The above framework can be extended to non-conjugate models having multiple functions, such as multi-class GP classification, by introducing separate v parameter per GP function. In our experiments, we consider only single- function non-conjugate GP models and we leave the experi- mentation with more complex models for future work.
N
Σ B p(u:.ff) - KL(g(u)llx()),
(20)
where g(f)) f p( is the marginal over J : f(x:) with respect the approximate posterior q(f, u). Given that q(u) is Gaussian with mean and co- variance S. ((f)) can be computed fast in O 2M2) time (after precomputing the Cholesky factorization of K as follows
4. Related Work
Several recent works on sparse GPs focus on constructing efficient inducing points, such works that place inducing points grid (Wilson & Nickisch, 2015; Evans & Nair, 2018; Gardner et al., 2018), construct inter-domain Fourier features (Lázaro- Gredilla & Figuciras- Vidal, 2009; Hens- man et al., 2018), 8), provide Bayesian treatments to inducing inputs (Rossi et al., 2021) or use nearest neighbor sparsity structures (Tran et al., 2021: Wu al., 2022). There exist also algorithms that allow to increase the number of induc- ing points using the decoupled method (Cheng & Boots, 2017; Havasi et al., 2018) 8) and the related orthogonally de- coupled approaches (Salimbeni et al., 2018; Shi et al., 2020; Sun et al., 2021; Tiao et al., 2023). Our contribution is orthogonal to these previous methods since we relax the conditional GP prior assumption in the posterior variational
((f.) [N ( m. kis q: + kj.uu K...uK"SKK.. kufe).
(21) For the discussion next it is useful to observe that the effi- ciency when computing ((ft) comes from p(flu) being conditional GP prior, so expressing p(fa u) is trivial.
Suppose now that we wish to impose the more struc- tured variational approximation q(f, u) = ((((() where q(u) = N(ul,, j, S) and (f(u)) is given by Equation (13). The bound can be written
1 N
E E
Σ (v: log Us- 1)
2
KL ((()))))()))
(22)
5

---

New Bounds for Sparse Variational Gaussian Processes
approximation. This means that our method could be used to improve previous variational sparse GP approaches, as the ones mentioned above as well earlier schemes that select inducing points from the training inputs (Cao et al., 2013; Chai, 2012; Schreiter et al., 2016).
5. Experiments
5.1. Illustration in 1-D Regression
In the first regression experiment we consider the 1-D Snel- son dataset (Snelson & Ghahramani, 2006). We took a subset of 40 examples of this dataset and we fitted the exact GP with the squared exponential kemel k(x,z')
Zhu et al. (2023) proposed inducing points GP approxima- tions that change the conditional GP p(flu) in the variational approximation to a modified conditional GP that uses differ- ent kernel hyperparameters in its mean vector. Note that our method differs since our (fflu) directly tries to construct better approximation to the exact posterior p(f]u, y), using the extra V variational parameters, without changing the kernel hyperparameters: see Section 3. More importantly, our method has ((N cost, while the ELBO in Zhu et al. (2023) (see Section 3.1 and Appendix A.1 in their paper) has cubic cost ((N3) since it depends on the inverse of K Qrr (denoted as K,,, in their paper).
(opp- We also fitted sparse variational GPs of exp(- 2 with either the standard collapsed bound (Titsias, 2009) from Equation (8) (SGPR) or the new collapsed bound from Equation (16) (SGPR-new). Both sparse GP methods use seven inducing points initialized at the same values as shown in Figure 1. All methods are initialized to the same hyper- parameter values; see Appendix D.
Figure 1 shows the results. Note that both SGPR and SGPR- new find similar inducing point locations. But SGPR-new, as a tighter bound (see panel (e)), is able to reduce some bias when estimating the hyperparameters since finds a noise variance o² closer to the one by exact GP (see panel (f)) This results in better predictions that match better the exact GP, as shown by the comparative visualization in panel (d). From panel (d), observe that both the mean and variances of SGPR-new are closer to the exact GP than SGPR.
Artemev et al. (2021) derived an upper bound on the log determinant log Kr + o²I] in the exact GP log marginal likelihood and obtained the following tighter upper bound to the initial trace regularization term 2rz tr (Kr- Qm):
N
log (1 1 2
tr( (k Qn)
(24)
No²
5.2. Medium Size Regression Datasets
To further investigate the findings from the previous section, we consider three medium size real-world UCI regression datasets (Pol, Bike, and Elevators) with roughly 10k train- ing data points each, and for which we can still run the exact GP. We choose the ARD squared exponential kernel kt x²) di exp( - Σ (w,-=()²² We run all three previous methods (Exact GP, SGPR, SGPR-new) five times with different random train-test splits; see Appendix D for experimental details. We also include in the comparison fourth method (discussed in Related Work) which is the Artemev et al. (2021)'s bound (SGPR-artemev) that does training using the collapsed bound from Equation (36) in Appendix B.4. All sparse GP methods use M = 1024 or M 2048 inducing points initialized by k-means. Figure (in the first two lines) shows the objective function and the noise variance 6² across 10k optimization steps using Adam with base learning rate 0.01 and for M = 1024. For SGPR-new, the third line in Figure 2 shows histograms of the estimated final values of the optimal e; variational parameters. Figure 4 in Appendix D. shows the correspond- ing plots for M = 2048. We observe that for Pol and Bike, SGPR-new matches closer the exact GP training than SGPR and SGPR-artemev. Specifically, SGPR-new gives higher ELBO and estimates the noise variance with reduced under- fitting bias. For the Elevators dataset, M = 1024 inducing points were enough for sparse GP methods to closely match exact GP training. This happens because in this case Qir
Our bound is tighter since from Jensen's inequality holds
log 1 (r( S E log 1 1 hus
a²
Further, the above regularization term can be interpreted as a restricted special case of our method, obtained through a (ffu) from Equation (13) where the diagonal matrix V constrained to be spherical V = vI; see Appendix B.4. Finally note, that unlike (24) (where the sum is inside the logarithm) our bound allows to apply stochastic optimiza- tion as described in Section 3.2.
Finally, Bui et al. (2017) used power expectation propaga- tion that minimizes a-divergence and derived an approxi- mation to the log marginal likelihood that interpolates be- tween the FITC (a 1) log marginal likelihood (Snelson & Ghahramani, 2006; Muinonero-Candell & Rasmussen, 2005) and the standard collapsed variational bound in (8) (a - 0). This approximation uses the regularization term
N 1 a
1 a log
Kus - qts
(25)
20
g²
This is different from ours since there is no value of a such that the two regularization terms will become equal. For example, note that for a - O. Equation (25) reduces to tr tr (Kn- - Qn) as discussed in Bui al. (2017).
6

---

New Bounds for Sparse Variational Gaussian Processes
++ Ho4
(a)
(b)
(c)
-20
0.30
ELBO or log- -l-llellhhodd
-25
Noise Vartance
0.25
0.20
-30
-35
Exact
0.15
SGPR
0.10
-40
SGPR-new
200
400
600
800
1000
0
200
400
600
800
1000
Iteration
Iteration
(c)
(1)
Figure 1. First row shows posterior predictions (means with 2-standard deviations) after fitting the exact GP (a), and the sparse GPs with either the standard collapsed SGPR bound (b) the proposed SGPR-new collapsed bound (c). In panels (b).(c) the seven inducing points are initialized to the same random locations (shown on top with crosses) while the optimized values are shown at the bottom. Panel (d) superimposes all predictions in order provide marginal likelihood for the exact GP) values across optimization steps while (f) shows the corresponding values for the noise variance ?
comparative visualization. Finally, panel (e) shows the ELBO exact log
5.3. Large Scale Regression Datasets
Table 1. Average test log- tiglliiklihoods for the medium size regres- sion datasets. The numbers in parentheses are standard errors.
We consider 8 regression datasets, with training data sizes ranging from tens of thousands to millions. We implemented the stochastic optimization versions of the two scalable sparse GP methods: (i) the one that trains using the previous uncollapsed bound from Hensman et (2013) (SVGP) and (ii) our new bound from Equation (18) (SVGP-new), We denote these stochastic optimization versions by SVGP to distinguish them from the corresponding SGPR methods that use the more expensive collapsed bounds. We run the SVGP methods with M = 1024 and 2048 inducing points, Matern3/2 kernel with common lengthscale, minibatch size 1024, Adam with base learning rate 0.01 and 100 epochs. These experimental settings match the ones in Wang et al. (2019) and Shi al. (2020) further described in Ap- pendix D.2. Table 2 reports the test log likelihood scores for all datasets. In the comparison we also included two strong baselines from Table 2 in Shi et al. (2020), i.e., SOLVE-GP and ODVGP (Salimbeni et al., 2018).
Pol
Bike
Elevators
Exact GP
1.089(0.011)
(.100().022)
0.386(0.001)
M 1024
SGPR
0.821 (0.008)
2.176(0.020)
0.387(0.001)
SGPR-artemev
0.859(0.007)
2.199(0.024)
-0.387(0.001)
nGPR-new
0.920(0.006)
2.326(0.026)
-0.387(0.001)
= 2048
SGPR
0.958(0.008)
(.337().030)
00.87(00.00)
SGPR-artemev
0.976(0.008)
2.356(0.029)
0..87((0.001)
SGPR-new
0.998(0.008)
(0.11(
-0.887(0.001)
accurately approximates Krr, i.c., the elements ky qu get close zero. For this latter dataset, observe that since the ki 4s values are close to zero the corresponding y; values are concentrated around one as shown by the corresponding (right-most) histogram in Figure 2.
From the predictive log likelihood scores in Table 2 and also the corresponding Root Mean Squared Error (RMSE) scores reported in Table 4 in Appendix D.2, that training with the new SVGP-new variational bound provides a clear improvement compared to training with the previous SVGP bound. Note that this improvement requires no change in the computational cost, and fact there
conclude
Table reports test likelloodd predictions which show that SGwR-new outperforms SGPR and aGPR-artemev.

---

New Bounds for Sparse Variational Gaussian Processes
1e4
1e4
- -0.4 led
Elevators
ELBO log oog-likellood
0.5
2
- 0.6
0.0
Exact GP
SGPR
0
- 0088
-0.5
SGPR-new
-0.9
SGPR-artemev
-1.0
-1.0
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.8
1.0
0.0
0.2
0.6
0.8
1.0
Iteration
Iteration
led
Iteration
1.01e-1 1e-1
1.0
1.0
1e-2
3.0 1e-1
Noise variance
0.8
2.5
0.6
0.6
2.0
0.4
0.4
1.5
0.2
0.2
0.0
1.0
0.0
0.0
0.2
0.4
0.6
0.8
1.0
0.2
0.4
0.8
1.0
0.0
0.2
0.4
0.6
1.0
Iteration
le4
Iteration
1e4
Iteration
1e4
1e2
1e3
1e4
1.5
1.0
0.8
0.6
0.4
0.5
0.2
8.0
0.4
Vi 0.8
0.6
1.0
0.0.0
0.2
0.4
0.6
1.0
0.0.0
0.2
0.4
0.6
0.8
1.0
V
Figure 2. The three plots in each column correspond to the same dataset: first row shows the ELBO (or log-likelihood) for all four methods (Exact GP, SGPR, SGPR-new and SGPR-artemev) with number of iterations, and the plot in the second row shows the corresponding values for o², SGPR methods use M 1024 inducing points initialized by k-means. For these two first lines plot the mean and standard error after repeating the experiment five times with different train-test dataset splits; see Appendix for further experimental details. For one of the runs of SGPR-new, the third line shows histograms for the estimated values of the variational parameters U; = (1 ki
only minor modification needed to be done in existing SVGP implementation in order to run SVGP-new.
are generated using Poisson intensities given by A(x) 3.5 + 3 sin(z).. We train the GP model with the SVGP bound and the proposed SVGP-new bound using 6 inducing points initialized to the same values for both methods; see Appendix D.3. Figure 3(left) shows the observed counts together with the predictions obtained by SVGP, SVGP-new and non-sparse variational GP (Full GP). From this figure and from the ELBO values, we observe that SVGP-new remains closer Full GP.
5.4. Poisson Regression
We consider non-Gaussian likelihood example where the output data counts modeled by Poisson likelihood p(y(f) = i= w.T follow a GP prior. For such case the new variational approx- imation includes single additional variational parameter denoted by which is optimized together with the remain- ing parameters; Section 3.3. We will compare training with the new ELBO from Equation (23) (we denote this method by SVGP-new) with the standard ELBBO that is obtained by restricting v = (SVGP).
eli -e' where the log intensities values
Secondly, we consider a real dataset (NYBikes) about bicy- cles crossings going over bridges in New York City². This dataset a daily record the number of bicycles crossing into or out of Manhattan via one the East River bridges over a period months. The data contains 210 points and
Firstly, we consider an artificial example of 50 observations with 1-D inputs placed in the grid [-10, 10] where counts
"This datasct is freely available from https: www. www.:kagglc..o com wawpaasttt/nooyooo-ccee nyc-east- -aw.aaa----ee---ooooee--eeeeeeee
8

---

New Bounds for Sparse Variational Gaussian Processes
Table Test lng-likelihoods for the large scale regression datasets with standard errors in parentheses. Best mean values highlighted.
Kin40k
Protein
KegeDirected
KEGGU
3dRoa
Song
Buzz
HouscElectric
N
25,600
29,267
31,248
40,708
278,319
329,820
373,280
1.311,539
8
20
27
90
77
9
From Shi et (2020)
ODVGP
1024 + 1024 0.137(0.003)
0.956(00005)
0.119(0.0067) 0.105(0.033)
-0.664(0.003)
-1.193(0.001)
-0.078(0.001)
1.317(0.002)
1021 + 8096 0.144(0.002)
-0.946(0.005)
-0.136(0.063) 0.109(0.033)
-0.657(0.003)
-1.193(0.001)
-0.079(0.001)
1.319(0.004)
SOLVE-GP
1021 + 1024 0.187(0.002)
-0.943(0.00S)
0.973(0.003)
0.680(0.003)
-0.659(0.002)
-1.192(0.001)
-0.071(0.001)
1.333(0.003)
SVGP
1024
0.108(0.002)
-0.969(0.006)
1.042(0.009) 0.699(0.005)
-0.704(0.003)
-1.192(0.001)
-0.069(0.002)
1.383(0.002)
2048
0.237(0.002)
- 0.961(000006)
1. 1.050(0000)
0.703(0.005)
-0 630(0.003)
-1.190(0.001)
-0.063(0.00)
1.419(0.000 002)
SVGP-new
1024
0.152(0.003)
-0 (.965(0.006)
1.044(0.009)
0.699(0.005)
-0.701 (0.00))
-1.192(0.001)
-0.065(00.00)
1.387(0.003))
2048
0.286(0.002))
0.938(0.000 006)
1.051(0.009) 0. 0.703(0.005)
-0.65((0.004)
11.10000.001)
0.060(0.001)
(0.00(
1e4
-100
-10
120
-1.0
Count
10
130
1.5
144
-2.0
Full GP
SVGP SVGP-new
150
-2.5
-10
-5
0 Input
10
160
500
1000 Itcration
1500
2000
-3.0
2
3
Iteration
1e4
Figure 3. (left) shows the predictions (means with 2-standard deviations) over counts (black dots) in the artificial data example after fitting the Full GP, and the two SVGPs. This plot superimposes all predictions in order to provide comparative visualization. (middle) shows the ELBO across optimization steps for the artificial data example. (right) shows the ELBO for the NYBikes dataset and M 16.
we randomly choose 90% for training and 10% for test. We apply GP Poisson regression for the Brooklyn bridge counts where the input vector æ is taken to be two-dimensional con- sisted of maximum and minimum daily temperatures. We train the sparse GPs with either SVGP or SVGP-new and with M = 8, 16, 32 inducing points initialized by k-means. Since the dataset small we also run the non-sparse Full GP. The ELBO across iterations in Figure 3 (right) and the test log likelihood scores (Table 5 in Appendix D.3) indicate that SVGP-new provides better approximation than SVGP.
Acknowledgements
am grateful to the reviewers for their comments. Also, wish to thank Jiaxin Shi and Francisco Ruiz for their invaluable advice and generous help during this project.
Impact Statement
This paper presents work whose goal is to advance the field of Machine Learning. There are many potential socictal consequences of our work, none which we feel must be specifically highlighted here.
6. Conclusions
We have presented method that relaxes the conditional GP assumption in the approximate distribution in sparse varia- tional GPs. This leads to tighter collapsed and uncollapsed bounds, that maintain the computational cost with the previ- ous bounds and can reduce training underfitting. For future work an interesting topic is to apply our method to more complex GP models, such as those with multiple outputs, with uncertain inputs and deep GPs. For the Bayesian GP- LVM, where the collapsed closed form bound has strong similarities with the previous GP regression collapsed bound in Equation (8), deriving new collapsed bound is tractable as described in Appendix B.5. Finally, it might be useful to investigate whether theoretical convergence results on sparse GPs (Burt et al., 2020; Wild et al., 2023), can be improved given the new collapsed lower bound.
References
Artemev, A., Burt, R., and van der Wilk, M. Tighter bounds the log marginal likelihood of gaussian process regression using conjugate gradients. In Meila, M. and Zhang, T. (eds.), Proceedings of the 38th International Conference on Machine Learning (ICML), Proceedings Machine Learning Research, pp. 362-372. PMLR, Jul 2021.
Banerjee, S., Gelfand, A. E.,, Finley, A. O., and Sang, H. Gaussian predictive process models for large spatial data sets. Journal of the Royal Statistical Society Series B: Statistical Methodology, 70(4):825-888, 07 2008
Bauer, M., van der Wilk, M., and Rasmussen, C. E. Under-

---

New Bounds for Sparse Variational Gaussian Processes
standing probabilistic sparse gaussian process approxima- tions. In Lee, D., Sugiyama, M., Luxburg, U., Guyon, I,, and Garnett, R. (eds.), Advances in Neural Information Processing Systems, volume 29. Curran Associates, Inc., 2016.
de G. Matthews, A. G., van der Wilk, M., Nickson, T., Fujii, K., Boukouvalas, A., León-Villagrá, P., Ghahramani, Z.,, and Hensman, J. Gpflow: A gaussian process library using tensorflow. Journal of Machine Learning Research, 18(40):1-6, 2017.
Bui, T. D., Yan, J., and Turner, R. E. A unifying framework for Gaussian process pseudo-point approximations us- ing power expectation propagation. Journal of Machine Learning Research, 18(104):1-22 2017 2017.
Deisenroth, M. and Rasmussen, C. E. PILCO: model- based and data-efficient approach to policy search. In International Conference on Machine Learning, pp. 465- 472, 2011.
Bui, T. D., Ashman, M., and Turner, R. E. Tighter sparse variational gaussian processes. Transactions on Machine Learning Research, 2025.
Dezfouli, A. and Bonilla, E. V. Scalable inference for gaus- sian process models with black-box likelihoods. In Cortes, C., Lawrence, N., Lee, D., Sugiyama, M., and Garnett, R. (eds.), Advances Neural Information Processing Systems, volume 28. Curran Associates, Inc., 2015.
Burt, R., Rasmussen, C. E., and van der Wilk, M. Con- vergence of sparse variational inference in gaussian pro- cesses regression. Journal of Machine Learning Research, 21(131):1-63, 63,, 2020.
Evans, T. and Nair, P. Scalable Gaussian processes with grid- structured eigenfunctions (GP-GRIEF). In International Conference on Machine Learning, pp. 1416-1425, 2018.
Cao, Y. Brubaker, M. A., Fleet, D. J., and Hertzmann, A. Efficient optimization for sparse gaussian process regres- sion. In Burges, C., Bottou, L., Welling, M., Ghahra- mani, Z., and Weinberger, K. (eds.), Advances in Neural Information Processing Systems, volume 26. Curran As- sociates, Inc., 2013.
Finley, A. O., Sang, H., Banerjee, S., and Gelfand, A. E. Improving the performance of predictive process model- ing for large datasets. Computational Statisties and Data Analysis, (388)28873-888 2009.
Chai, K. M. A. Variational multinomial logit gaussian pro- cess. Journal of Machine Learning Research, 13(56): 1745-1808, 2012.
Gardner, J., Pleiss, G., Wu, R., Weinberger, K., and Wilson, A. Product kernel interpolation for scalable Gaussian processes, In Artificial Intelligence and Statistics, pp. 1407-1416, 2018.
Cheng, C.-A. and Boots, B. Variational inference for Gaus- sian process models with linear complexity. In Advances in Neural Information Processing Systems, pp. 5184- 5194, 2017.
Garnett, R. Bayesian Optimization. Cambridge University Press, 2023.
Gramacy, R. B. Surrogates: Gaussian Process Modeling, Design and Optimization for the Applied Sciences. Chap- man Hall/CRC, 2020.
Cressie, N. Statistics for spatial data. John Willey & Sons, New York, NY, 1993.
Csato, L. and Opper, M. Sparse online Gaussian processes. Neural Computation, 14:641-666 2002.
Havasi, M., Hernándzz-Lobtto, J. M., and Murillo-Fuentes, J. J. Deep Gaussian processes with decoupled inducing inputs. arXiv preprint arXiv: 02939, 2018
Damianou, A. and Lawrence, N. D. Deep Gaussian pro- cesses. In Carvalho, C. M. and Ravikumar, P. (eds.), Proceedings of the Sixteenth International Conference on Artificial Intelligence and Statistics, volume 31 of Pro- ceedings of Machine Learning Research, pp. 207-215, Scottsdale, Arizona, USA, 29 Apr-01 May 2013. PMLR.
Heaton, M., Datta, A., Finley, A., Furrer, R., Guinness, J., Guhaniyogi. R., Gerber, F., Gramacy, R., Hammerling, D., Katzfuss, M., Lindgren, F., Nychka, D., Sun, F., and Zammit-Mangion, A. A case study competition among methods for analyzing large spatial data. Journal of Agri- cultural, Biological and Environmental Statistics, 24, 12 2018.
Damianou, A. C., Titsias, M. K., and Lawrence, N. Vari- ational inference for latent variables and uncertain inputs in gaussian processes. Journal of Machine Learning Re- search, 17(42):1 12, 2016.
Hensman, J., Fusi, N., and Lawrence, N. D. Gaussian processes for big data. arXiv preprint arXiv:1309.68)) 2013.
de G. Matthews, A. G., Hensman, J., Turner, R. E.,, and Ghahramani, Z. On sparse variational methods and the kullback-leibler divergence between stochastic processes. In International Conference on Artificial Intelligence and Statistics, 2016.
Hensman, J., Matthews, A., and Ghahramani, Z. Scalable variational Gaussian process classification. In Artificial Intelligence and Statistics, pp. 351-360, 2015.
10

---

New Bounds for Sparse Variational Gaussian Processes
Hensman, J., Durrande, N., and Solin, A. Variational fourier features for gaussian processes, Journal of Machine Learning Research, 18(151):1-52, 2018.
Rossi, S., Heinonen, M., Bonilla, E., Shen, Z., and Filippone, M. Sparse gaussian processes revisited: Bayesian ap- proaches to inducing-variable approximations. In Baner- jee, A. and Fukumizu, K. (eds.), Proceedings of The 24th International Conference on Artificial Intelligence and Statistics, volume 130 Proceedings of Machine Learn- ing Reseurch, pp. 1837-1845. PMLR, 13-15 Apr 2021.
Jones, D. R.,. Schonlau, M., and Welch, J. Efficient global optimization of expensive black-box functions. Journal of Global Optimization, 13:455-492 1998.
Lawrence, N. Probabilistic non-linear principal component analysis with Gaussian process latent variable models. Journal of Machine Learning Research, 6(Nov):1783- 1816, 2005.
Salimbeni, H. and Deisenroth, M. Doubly stochastic vari- ational inference for deep Gaussian processes. In Ad- vances in Neural Information Processing Systems, pp. 4588-4599, 2017.
Lawrence, N. D., Seeger, M., and Herbrich, R. Fast sparse Gaussian process methods: the informative vector ma- chine. In Advances in Neural Information Processing Systems. MIT Press, 2002.
Salimbeni, H., Cheng, C.-A., Boots, B., and Deisenroth, M. Orthogonally decoupled variational Gaussian processes. In Advances in Neural Information Processing Systems, pp. 8711-8720, 2018.
Schreiter, J., Nguyen- Souyen-Tuong, D., and Toussaint, M. Efficient sparsification for gaussian process regression. Neurocom- puting, 192:29-37, 2016.
Gredill,,,,,,la, M. and Figuciras-Videe A. Inter-domain gaussian processes for sparse inference using inducing features. In Bengio, Y., Schuurmans, D., Lafferty, J., Williams, C., and Culotta, A. (eds.), Advances in Neural Information Processing Systems, volume 22, 2009.
Seeger, M. W., Williams, C. K. I., and Lawrence, N. D. Fast forward selection speed up sparse gaussian process regression. In Bishop, C. M. and Frey, B. J. (eds.), Pro- ceedings of the Ninth International Workshop on Artificial Intelligence and Statistics, volume R4 of Proceedings of Machine Learning Research, pp. 254-261. PMLR, 03-06 Jan 2003.
Leibfried, F., Dutordoir, V., John, S., and Durrande, N. A tutorial on sparse gaussian processes and variational inference. arXiv preprint arXiv:2012.13962 2022.
Liu, H., Ong, Y., Shen, X., and Cai, J. When gaussian process meets big data: A review of scalable gps. IEEE Transactions on Neural Networks and Learning Systems, PP:1-19, 01 2020. doi: 10.1109/TNNLS.2 MINNNLS.2019.299 2019.2957109
Sheth, R., Wang, Y., and Khardon, R. Sparse variational inference for generalized gp models. In Bach, and Blei, D. (eds.), Proceedings of the 32nd International Confer- ence on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pp. 1302-1311, Lille, France, 07-09 Jul 2015. PMLR.
Lloyd, C., Gunter, T., Osborne, M., and Roberts, S. Varia- tional inference for gaussian process modulated poisson processes. In Bach, and Blei, D. (eds.), Proceedings of the 32nd International Conference on Machine Learning, volume 37 of Proceedings of Machine Learning Research, pp. 18 1822, Lille, France, 07-09 Jul 2015. PMLR.
Shi, J., Titsias, M. K., and Mnih, A. Sparse orthogonal variational inference for gaussian processes. In Chiappa, S. and Calandra, R. (eds.), Proceedings of the Twenty Third International Conference Artificial Intelligence and Statistics, volume 108 of Proceedings of Machine Learning Research, pp. 1942. PMLR, 26-28 Aug 2020.
Nguyen, T. V. and Bonilla, E. V. Collaborative multi-output gaussian processes. In Proceedings of the Thirtieth Con- ference on Uncertainty in Artificial Intelligence, UAI'14, pp. 643-652, 2014.
Snelson, E. and Ghahramani, Z. Sparse Gaussian processes using cocudo-nippts... In Weiss, Y., Schölkopf, B., and Platt, J. C. (eds.), Advances in Neural Information Pro- cessing Systems, pp. 1257-1264. 2006.
O' Hagan, A. Curve fitting and optimal design for predic- tion. Journal of the Royal Statistical Society: Series B (Methodological), 40(1) 11:1114244 1978.
Quinonero-Candela J. and Rasmussen, C. E. A unifying view of sparse approximate Gaussian process regression. Journal of Machine Learning Research, 6:1939-1959, 2005.
Sun, S., Shi, J., Wilson, A. G. G., and Grosse, R. B. Scal- able variational gaussian processes via harmonic kernel decomposition. In Meila, M. and Zhang, T. (eds.), Pro- ceedings of the 38th International Conference on Ma- chine Learning, volume 139 of Proceedings of Machine Learning Research, pp. 9955-9965. PMLR, 18-24 24 Jul 2021.
Rasmussen, E. and Williams, C. K. Gaussian Processes for Machine Learning. MIT Press, 2006.

---

New Bounds for Sparse Variational Gaussian Processes
Tiao, L. C., Dutordoir, V., and Picheny, V. Spherical in- ducing features for aaateggonnllyyeecoupie Gaussian pro- cesses. In Krause, A., Brunskill, E., Cho, K., Engelhardt, B., Sabato, S., and Scarlett, J. (eds.), Proceedings of the 40th International Conference on Machine Learning, vol- ume 202 of Proceedings of Machine Learning Research, pp. 34143-34160, PMLR, 23-29 Jul 2023.
Zhu, X., Wu, K., Maus, N., Gardner, J., and Bindel, D. Varia- tional gaussian processes with decoupled conditionals. In Oh, A., Naumann, T., Globerson, A., Saenko, K., Hardt, M., and Levine, S. (eds.), Advances in Neural Informa- tion Processing Systems, volume 36, pp. 46191-46211 Curran Associates, Inc., 2023.
Álvarez, M., Luengo, D., Titsias, M. K., and Lawrence, N. D. Efficient multioutput gaussian processes through varia- tional inducing kernlls. In Teh, Y. W. and Titterington, M. (eds.), Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, vol- ume 9 of Proceedings of Machine Learning Research, pp. 25-32, Chia Laguna Resort, Sardinia, Italy, 13-15 May
Titsias, M. K. Variational learning of inducing variables in sparse Gaussian processes. In Artificial Intelligence and Statistics, pp. 567- 517-5777-755 2009.
Titsias, M. K. and Lawrence, N. D. Bayesian gaussian pro- cess latent variable model. In Teh, Y. W. and Titterington, M. (eds.), Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, vol- ume 9 Proceedings of Machine Learning Research, pp. 844-8 851, Chia Laguna Resort, Sardinia, Italy, 13-15 May 2010. PMLR.
2010. PMLR.
Tran, G.-L., Milios, D., Michiardi, P., and Filippone, M. Sparse within sparse gaussian processes using neighbor information. In Meila, M. and Zhang. I. (eds.), Proceed- ings of the 38th International Conference Learning, volume 139 of Proceedings of Machine Learn- ing Research, pp. 10369- 18-24 Jul 2021.
Machine
Wang, A., Pleiss, G., Gardner, R., Tyree, S., Wein- berger, K. Q., and Wilson, A. G. Exact Gaussian processes million data points. arXiv preprint 1rXiv:1903.08114 2019.
Wild, V., Kanagawa, M., and Sejdinovic, D. Connections and equivalences between the nyström method and sparse variational gaussian processes, 2023.
Wilson, A. and Nickisch, H. Kernel interpolation for seal- able structured Gaussian processes (KISS-GP). In In- ternational Conference on Machine Learning, pp. 1775- 1784, 2015.
Wu, L., Pleiss, G., and Cunningham, J. P. Variational nearest neighbor Gaussian process. In Chaudhuri, K., Jegelka, S., Song, I.,, Szepesvari, C., Niu, G., and Sabato, S. (eds.), Proceedings of the 39th International Conference Machine Learning, volume 162 of Proceedings of Machine Learning Research, pp. 14114-244100 PMLR, 17-23 Jul 2022.
Yousefi, F., Smith, M. T., and Álvarez, M. Multi-task learn ing for aggregated data using gaussian processes. In Wal- lach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc Bucc F,, Fox, E.,, and Garnett, R. (eds.), Advances in Neural Information Processing Systems, volume 32. Curran As- sociates, Inc., 2019.
12

---

New Bounds for Sparse Variational Gaussian Processes
A. Further details about the SVGP method
We give a brief overview of the derivation the standard collapsed bound in Equation (8). Some steps the derivation will also be instructive for proving the main results this paper in Appendix B.
Given the variational distribution g(f, p(flu) the lower bound is
/ ((flu))(()) log p(y())(Fa)p(u) dfdu
log py) 2
p(fHa)q(u)
/ p(t(u)q(u)lo) log
p(y(f)p(u)
dfdu
q(u)
/ a(u) { p(f (oggly/ff))ttteoo log p(u) qu))
du
f ( log exp(ʃ ((fu)) ply p) du.
(26)
(u))
The expectation f pf(ju) log p(y f)df can be computed as f [D((u)) log logp(y(f)df
p(f(u) logp(y)f)ef
f p(flu) {
1 - N log(2mo²)
20 yyy" 2yf ff df
N log(2no²)
20 tr [yy" 2y 2y u) + (Kra K (maaKy'
+ Kn
Qn]
= log [M(y Km Ku o 0D]]]
20 tr(K(r Qir)-
(27)
where we highlighted with blue term the third line contrast it with similar term when proving Lemma 3.2 in Appendix B.4. The ELBO Equation (26) is written as
log ply) 2 f q(() log Kr T))(( df 20
tr (Kr Qm)-
(28)
g(u)
By maximizing this bound wrt the distribution (u) obtain the optimal 9":
(u) N (y
()) N (y u, o I))pp [N Km u, o² o'))p(uddu
(29)
N(y (0, QH I2
= Nulo " (Km + o Kntknn) -aacKa))Kutt Kmn (Kuu + a 2 Kar "Kn),
(30)
where the expression the second line (obtained after some standard completion of square procedure) shows that g (u( can be computed in O(N2²) time. In fact, this optimal 4" (u) the same as the one obtained by the DTC (also known as projected process) approximation (Seeger et al., 2003; Quinonero-Candelaa Rasmussen, 2005). By substituting the expression in (29) into the bound in (28) obtain the well-known formula of the collapsed bound:
log p(y) 2 Ny Qu + o²1) 202 KH - Qr)-
(31)
Given the Gaussian form of q(u) = (ulu,, A) the posterior GP is given by q(f.) = J p({.. pfF,/u)g(u)dd:
(f.)) = (f, Kr. μ, Kr,t. Kr,n Kut. Kr.uK AK Kut.)
(32)
which further simplifies if substitute the optimal mean and covariance of (u):
y (L.)) -T, -' Kut 22 Ky.t. Kr.u Kut. + Kr.uA Kut.)
(33)
where A Knn + o - Knr
13

---

New Bounds for Sparse Variational Gaussian Processes
B. Further details about the proposed bounds
Here, we provide several proofs regarding the proposed bounds.
B.1. Proof of Lemma 3.1
(f(j)) and pflu) are Gaussian distributions having the same mean but different covariance matrices. Thus, the KI. divergence reduces to
Kt
Qn
N + tr( Qui) Qn) Qm) Qm) - 2m)
(L
log
2 [Kr
Qm V (Kr Qm)
log V N + +v())
(34)
2
where all the terms involving Kn Qi cancel out by using standard properties of the matrix determinant and trace. Now since
V is a diagonal matrix (with diagonal elements > 0) we conclude that (L)(f = Ei= (2, log 1).
B.2. Proof of Lemma 3.2
The derivation of f ((fju) log (f) is similar to the derivation in Equation (27) with small difference highlighted in blue:
/ p(((u) logp(y/Ndff
N
2 log(2To")
202 yyy -2y--2yyeee 2y(Kuum''u) u) + (Kuukm'u (aluk,u')(aaan (Kuu, uu? + (Krr Q 2 Qm)
log [W K Km u, 01)
20 tr( (V(V Qm)),
(35)
where we used that tr( (K Qm V(r 2m) = tr(v - Qn)). Now since V a diagonal matrix we have
tr(V(Kn Qm)) = Em Vi(ks (a) which completes the proof.
B.3. Proof of Proposition 3.3
The ELBO is written
(f = Satu exp(Bytu pog(yy)))ll (f log 2 (f) log ((f
qu))
and by using the results from the two lemmas this becomes
logp(y) 2 / q(u)) log ((y]Kan K o T)p( du
1 2
his
Ti σ²
1
log
quu)
Clearly, maximizing over (u) gives the same optimal distribution in Equation (29), and the first term in the bound is the DTC log likelihood. The second term that depends on the r.s is concave function over these parameters. Thus, by
differentiating and setting to zero we obtain the optimal values oi (1
ky gei
we plug these values back into the
bound we obtain the new tighter collapsed bound in Proposition 3.3.
B.4. Reinterpretation of Artemev et al. (2021)'s bound
We consider the following form of (f]ju):
(fluu) (f u(Kr Q(r))-
Then, (Lf]((u N (v log 1) and (, [lo p(y(f(
2
log (y) Kr Km u, T) 20 K Qm) and the bound written as
logp(y) > f ((()) log
N ( K o ()) du
log 1)}
tr((Kr Qn) + NNo
q(u)
14

---

New Bounds for Sparse Variational Gaussian Processes
By maximizing wrt e we obtain (1 lr(Ke Qm)) and by substituting this back into the bound we obtain Artemev No²
et al. (2021)'s tighter bound on the initial trace regularization term. Overall this collapsed bound has the form
N poy 2 log. Nyyl, On + o² log 1 tr tr(Krr Qi))
(36)
2
No
This collapsed bound is what the method SGPR-artemev is using in Section 5.2. Note that Artemev et al. (2021) propose also additional but more expensive bounds for the first DTC log likelihood term that require running conjugate gradients. We do not consider those in our comparisons (such bounds could be used all SGPR bounds since all share the same DTC log likelihood term) as they have higher cost.
B.5. Bound for the Bayesian GP-INM
Due to the strong similarity of the standard collapsed SVGP bound in Equation (8) with the collapsed bound in the Bayesian GP-LVM (Titsias & Lawrence, 2010), applying the new approximation to Bayesian GP-LVM seems to be simple and we discuss next.
Given observed data Y E BRN and latent variables X RNNQ we have the latent variable model
p(((X)((X)- = (MP()ax)) p(X).
D
where p() a Gaussian prior over the latent variables and p(ya) X) = (Vs.. +02). Given a Gaussian variational distribution q(X) over the latent variables, the initial form of the bound is
F =
9 (X) log((Y(X)) - (L[((X)|Ix(X))
Σ
((X) ((x))opp(xx
KL({(X)IIb(X) )
Σ Fa - (L(((X))xx(
(37)
d=1
where F2 - (()) P(Xxx/XXXXX The KL part will be tractable KL between two Gaussians, and thus the difficulty is to approximate Fg. Given that log p(ya has the same form with the log marginal likelihood in GP regression, we can lower bound it using inducing variables and exactly the same form of (fluu) as we did in the main paper. This gives
logp(y2(X)2 2 / (((u) log Ny al Km (X (X)K''u (')) du 2
Σ
N
Us (1
1
k( x.,) g(0;, xei)
log Us - 1
(u)
a²
where x, is the latent variable for the i-th data point and g( 2.) = k(x;) T k(x,) = tr(k (zx (<x;)) T Note that we write the cross kernel matrix as K₁, (X) to emphasize its dependence on the latent variables X, while K does not depend on X. Note also that we assume that each U; parameter does not depend on 2; and this is crucial to obtain closed form collapsed bound. Now follow the derivation in the initial Bayesian GP-LVM where we use the above bound to replace log (X) in f x(X) log p(ya| xlX)d )dX and do first the expectation over X, and then solve for the optimal (u)) This eliminates (u) and it gives the bound
log / e (eiis,)K (X) 2 Maxipp 2
N
(k(2,, x; (())) (q(2,, Bi)))(ma) )
log vi - }
0: 1
where we have used physics notation for expectation, i.e., (.). For v; = this is the previous collapsed bound used by Bayesian GP-LVM. By maximizing over each v, we obtain the new collapsed bound
Fa log J CCNisalRa (a) () p(u)du
2 log 1+ (k(x; a()) - (g(x;; x. (()))))))))
1
N
x²
which can substituted the overall Bayesian GP-LVM bound above. Again the implementation of the new bound requires minor modification to existing Bayesian GP-LVM code.
15

---

New Bounds for Sparse Variational Gaussian Processes
C. Learned hyperparameters in 1-D Snelson dataset
Table 3 provides the learned hyperaarameters for the 1-D Snelson dataset.
Table Hyperparameter values in I-D Snclson example.
a²
0j
(
Exact GP
0.0715
0.712
0.597
SVGP-new
0.087
0,485
0.615
SVGP
0.108
0.331
0.617
D. Further experimental details and results
For all regression experiments (apart from the toy Snelson 1-D dataset) we repeat the runs for five times using different random training and test splits. By following Wang et al. (2019) and Shi et al. (2020) we consider 80% / 20% training / test splits. A 20% subset of the training set is used for validation.
The training inputs and regression outputs are normalized to have zero mean. For the hyperparameters o², o². (or lz for ARD kernels) we use the softplus activation to parametrize the square roots of these parameters, i.e., to parametrize a, o l. For all experiments we use the initializations a 0.51, of 0.69, & 1.0. The inducing inputs Z are initialized by running maximum 30 iterations of k-means clustering with the centers initialized at random training data subset.
D.I. Medium size regression datasets
For all three datasets this section we can run Exact GP given the medium training size. For the Pol dataset the training size is N 9600 and imput dimensionality d 26. For Elevators is N 10623 and d 18. For the Bike dataset the initial train size (see e.g., Table 7 in Shi et al. (2020) is N 11122 (with d 17) but since Exact GP training gave out-of-memory етror when running in a V100 GPU, we had to slightly reduce the training size to N = 10600.
An mentioned in the main paper the standard squared exponential ARD kemne was used in all experiments in this section. For training, we perform 10000 optimization iterations using the Adam optimizer with base learning 0.01.
Figure shows the objective function values and noise variance parameter o² across iterations when the SGPR methods use M 2048 inducing points. Figure 2 in the main paper shows the result for M 1024.
D.2. Large scale regression datasets
The experimental settings are chosen match the ones from Wang et al. (2019) and Shi et al. (2020), where we used GPs with Matérn32 kernel (with common lengthscale). Following these settings, for all datasets train for 100 epochs using Adam with learning rate 0.01 and minibatch size 1024.
Table 4 reports RMSE scores, while test log likelihood scores are given in Table 2 of the main paper.
D.3. Poisson regression
Figure 5 shows the ELBOs across iterations for the NYBikes dataset for M = 8 and M = 32 inducing points, while the plot for M 16 is shown in the main paper. Table 5 presents test log- likeliood scores for the NYBikes dataset. Average values and standard errors are computed by repeating the experiment five times where at each repeat we randomly split the initial data into 90% for training and 10% for test.
Regarding the scalar value of for the toy Poisson regression the learned value was around v = 0.675, while or the real Poisson example in NYBikes, the value gets very small below 0.01.
16

---

New Bounds for Sparse Variational Gaussian Processes
1e3
Pol
1e4
Bike
-0.4 1e4
Elevators
ELBO or log pe-likllhhood
2-
0.5
1
Exact GP
0.7
SGPR
0-
SGPR-new
SGPR-artemev
-1
-0.9
-1.0
0.0
0.2
0.4
0.6
0.8
1.0 le4
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
Iteration
lteration
le4
Iteration
le4
1.0 1e-1
1.0
1e-2
3.0
1e-1
0.8
0.8
2.5
2.0
0.4
0.4
1.5
0.2
0.2
0.0
1.0
0.0
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
0.0
0.2
0.4
0.6
0.8
1.0
Iteration
Iteration
1e4
1e3
1e3
1e4
3
1.00
0.8
Count
0.75
2
0.6
0.50
0.25
0.2
0.00
0.2
0.4
Vi 0.8
0.6
1.0
8.0
0.2
0.4
0.6
0.8
1.0
0.8.0
0.2
0.4
0.6
0.8
1.0
W
Figure 4. The two plots in each column correspond to the same dataset: first row shows the ELBO (or log- sigllikelihood) for all four methods (Exact GP, SGPR, SGPR-new and SGPR-artemev) with the number of iterations and the plot in the second row shows the corresponding values for o² SGPR methods use M = 2048 inducing points initialized by k-means. For these two first lines we plot the mean and standard after repeating the experiment five times with different train-test datasct splits. For one of the of SGPR-new, the third line shows histograms for the estimated values of the variational parameters t,.
Table Test RMSE values large scale regression datasets with standard errors in parentheses. Best mean values are highlighted.
Kin40k
Protein
KeggDirected
KEGGU
3dRoad
Song
Buzz
HouseBleccric
25,600
29.267
31,248
40,708
278.319
329,820
373,280
1,311,539
d
9
20
27
90
77
9
From Shi et al. (2020)
ODVGP
1024
1024
0.183(0.001)
0.625(0.004)
0.176(0.012)
0.156(0.004)
0.467(0.001)
0.797(0.001)
0.263(0.001)
0.062(0.000)
1024 8096
0.1 0.180(0.001)
0.618(0.004)
0.157(0.009)
0.157(0.004)
0.462(0.002)
0.797(0.001)
0.263(0.001)
0.062(0.000)
SOLVE-GP
1024 + 1024
0.172(0..01)
0.618(0.004)
0.095(0.002)
0.123(0.001)
0.464(0.001)
0.796(0.001)
0.261(0.001)
0.061(0.000)
1024
(0.95(0.001)
0.086(0.001) 0.635(0.004)
0.122(0.001)
0.486(0.002) (0.00))
0.707(0.001)
(.261() 0.059(0.000)
2048
0.171(0.000)
0.086(0.00) 0.619(0.004)
001)
0.121(00001)
0.460 (0.00)
0.260(0.001) (.795(0.001)
0.057(0.000)
SVGP-new
1024
(1.02(
0.086(0.00) 0.601().004)
0010))
0.122(0.001)
(..04(0.001)
0..06)).001)
0.259().001)
0.058(0.000)
2048
(0.08))
(0.04)
0.086(0.001)
(0.20()
0.461(0.002)
(0.00)
0. 0.208)
0.057(0.000)
17

---

New Bounds for Sparse Variational Gaussian Processes
1e4
1e4
-0.5
-1.0
-1.0
ELBO
1.5
1.5
-2.0
Full GP
-2.0
Full GP
SVGP
SVGP
-2.5
-2.
SVGP-new
SVGP-new
-3.0
-3.0
1
2
3
4
1
2
3
4
Iteration
1e4
Iteration
1e4
Figure 5. The left panel shows the lower bounds across iterations when the sparse GP methods run with M = inducing points, while the right panel shows corresponding plot with 32 inducing points.
Table 5. Test log likelihoods on the NYBikes Poisson regression dataset with standard errors in parentheses. For the sparse methods we consider varying numbers inducing points, i.e. M 8, 16,32.
Full GP
-5.061(0.010)
SVGP
8
36.397(6.017)
SVGP
-16.557(43307)
SVGP
32
-8.556(0.728)
SVGP-new
8
-9.713(0.345)
SVGP-new
16
-9.301 (0.296)
SVGP-new
(..00().
18


Reviews:
[{"review_id": "svtKTacNHL", "reviewer": "Reviewer_kDa9", "summary": "The authors present an improvement on the standard SVGP approximation by departing from the standard conditional GP prior distribution. The approach introduces an additional $N$ variational parameters which modify the covariance matrix of the conditional distribution. This leads to an improvement on the resultant bound for the log marginal likelihood. The authors demonstrate that their method works in the stochastic minibatch optimisation setting, and can be extended effectively to non-Gaussian likelihoods by judicious constraints on the additional variational parameters. The theoretical findings are shown to hold in practice on a number of small to large scale regression experiments.", "questions": "", "limitations": "", "rating": 5, "confidence": null, "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": "Thank you for very accurately describing the contribution of the paper and for pointing to concurrent work. As also mentioned in the response to Reviewer Pyzb below, we plan to discuss the concurrent work in the next version of our paper.", "claims_and_evidence": "Yes---claims are theoretical and demonstrated through experiments.", "methods_and_evaluation_criteria": "Yes.", "theoretical_claims": "Yes.", "experimental_designs_or_analyses": "Yes.", "relation_to_broader_scientific_literature": "N/A.", "essential_references_not_discussed": "N/A.", "other_strengths_and_weaknesses": "## Strengths\n\n* The proposed approach is novel (together with https://arxiv.org/abs/2502.04750 which, coincidentally, was released at the same time), straightforward to both understand and implement, and effective in practice, demonstrating improved or on-par performance relative to baselines across almost all experiments. Collectively, the paper is very convincing.\n\n\n## Weaknesses\n\n* A key strength of the proposed method is in its generality---I suspect that it can be applied to improve a wealth of SVGP approximations such as SOLVE-GP, and different forms of approximations such those used for SVGP-LVMs and Deep-SVGPs. I believe that further extensions in the paper would improve it further, although the authors do touch upon this as future work.", "other_comments_or_suggestions": "N/A.", "questions_for_authors": "N/A.", "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": "", "extra_scores": {}}, {"review_id": "7lryb79jf5", "reviewer": "Reviewer_Rdie", "summary": "The paper introduces new evidence lower bounds (ELBOs) for sparse variational Gaussian processes (SVGP) by relaxing the traditional assumption that the variational distribution must factorize with the conditional GP prior p(f|u). Instead, the authors propose a more flexible variational distribution q(f|u), which allows for a tighter bound. Theoretical analysis shows that the new bound is provably tighter than previous SVGP bounds, and experiments on regression and non-Gaussian likelihood tasks demonstrate improved hyperparameter learning and predictive performance. However, the proposed methods between theoretical analysis and practical experiments are slightly different, which reduces the soundness of this paper. The practical method is computationally efficient, requires minimal code modifications, and is compatible with stochastic optimization.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "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": "Thank you for your comments. Below we provide some responses. \n \n> The key insight, replacing $p(f|u)$ with a diagonal-covariance $q(f|u)$, is sound. However, the gap between diagonal V and spherical V is not thoroughly assessed.\n\nPlease note that in the medium-size regression experiments reported in Table 1 and Figure 2 we do compare with the method that assumes spherical $V = v^* I$ which uses the optimal scalar value $v^* = \\left( 1 + \\frac{\\text{tr}(K_{f f} - Q_{f f})}{N \\sigma^2} \\right)^{-1}$. This is precisely, the method denoted in the experiments as \"SGPR-artemev\".  We will clarify this further in the paper. Note, that currently we briefly explain the spherical case $V = v I$  for GP regression and the connection with Artemev et al's bound in Related Work (Section 4) and also in Appendix B.4.  From Table 1 and Figure 2 we can observe that diagonal V does work better than the scalar $v$.  \n\n> The insight presented in this paper is promising, but the details may require further discussion. The key question is: to what extent does assuming a spherical V influence the results? The potential risks associated with this assumption are not explored, nor is there a theoretical or experimental analysis addressing its impact.\n\nWe agree with the reviewer that the fact that the non-Gaussian likelihood case requires spherical $V = v I$ is not ideal.  In seems that for such non-Gaussian likelihoods the only option that works is to use a spherical $V = v I$ and, as explained in Section 3.3, the diagonal $V$ (that works for GP regression) has cubic cost when trying to obtain each marginal $q(f_i)$. Notice also that if someone heuristically tries to use diagonal $V$ for the non-Gaussian likelihood case and approximate the marginal $q(f_i)$ by\n$$\nq(f_i)  = \\mathcal{N}(f_i | k_{f_i u} K_{u u}^{-1} m, v_i (k_{ii} - q_{ii}) + k_{f_i u} K_{u u}^{-1} S K_{u u}^{-1} k_{u f_i}),\n$$\n(which is not the correct marginal under $q(f|u) q(u)$ since the variance term  $v_i (k_{ii} - q_{ii})$  is wrong) then this creates an inconsistency in the variational distribution in the ELBO, since the $q(f_i)$ used to compute the expected log-likelihood term will be inconsistent with the $q(f|u) q(u)$ in the KL divergence term $KL[q(f|u) q(u) || p(f|u) p(u)]$), and the objective is not a rigorous ELBO anymore. \n\nWe can add further discussion to clarify that the spherical $V= v I$ is needed for non-Gaussian likelihoods in order to obtain a rigorous ELBO.", "claims_and_evidence": "The claims are supported by clear theoretical derivations (Lemmas 3.1–3.3, Proposition 3.5) and extensive experiments. The theoretical proofs are straightforward, and experiments on synthetic/real-world datasets (e.g., Snelson, UCI datasets, NYBikes) validate reduced bias in hyperparameters (e.g., noise variance) and better test log-likelihoods.", "methods_and_evaluation_criteria": "The methods are appropriate for scalable GP inference. The evaluation uses standard metrics (test log-likelihood, RMSE) and datasets (UCI, Kin40k), with comparisons to common baselines (SGPR, SVGP, SOLVE-GP). Experiments include multiple trials to report standard errors, ensuring statistical reliability.", "theoretical_claims": "The key insight, replacing p(f|u) with a diagonal-covariance q(f|u), is sound. However, the gap between diagonal V and spherical V is not thoroughly assessed.", "experimental_designs_or_analyses": "The comprehensive experiments cover regression (Gaussian/non-Gaussian) and varying dataset sizes.", "relation_to_broader_scientific_literature": "The work builds on SVGP, positioning the new bound as a tighter alternative. This applies to all of the sparse GP utilizing the inducing points.", "essential_references_not_discussed": "NA", "other_strengths_and_weaknesses": "NA", "other_comments_or_suggestions": "Typos: \"hyperpaparameters\" (Abstract), \"minibathes\" (Section 3.2).", "questions_for_authors": "The insight presented in this paper is promising, but the details may require further discussion. The key question is: to what extent does assuming a spherical V influence the results? The potential risks associated with this assumption are not explored, nor is there a theoretical or experimental analysis addressing its impact.", "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": "", "extra_scores": {}}, {"review_id": "ap1Gx67ttD", "reviewer": "Reviewer_Pyzb", "summary": "The paper revisits the widely-used variational approximation for sparse Gaussian processes (GPs). It proposes a refined variational formulation, introducing a more flexible conditional posterior distribution in place of the traditional assumption (where the conditional posterior matches the prior). This adjustment results in a tighter variational bound on the marginal likelihood, improving inference quality. Additionally, the method naturally supports stochastic mini-batch optimization, making it scalable to large datasets and practical for a broader range of applications.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "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": "We would like to thank the reviewer. We respond to the main comments below.\n\n> Yes, mainly. The proposed methods and evaluation criteria are appropriate for the problem at hand, and the experiments effectively demonstrate the benefits of the tighter variational bound. However, the datasets used in the evaluation seem somewhat simple, and the paper could have further strengthened its empirical validation by including a more complex benchmark, such as MNIST or a real-world dataset. This would provide additional evidence of the method’s scalability and effectiveness in practical applications.\n\nThank you for this comment. Here, we have focused mainly on the GP regression, large GP regression and we provide \nalso an experiment with Poisson regression which is a non-Gaussian likelihood example.   \n\n> The main weakness of the paper is that the results demonstrate only marginal improvements with the new bound across all evaluated tasks, without a clear example where the standard approach would fail without it. This makes it difficult to assess the practical necessity of the proposed refinement. Additionally, in the non-Gaussian likelihood setting, the paper opts for a simpler approximation rather than leveraging the more expressive variational form, seemingly due to computational constraints.\n\nPlease note that we do observe noticeable improvements in predictive performance in some experiments (Pol, Bike, Kin40k, Protein, Buzz). For some datasets, like Pol and Kin40k, the improvement is significant. We believe that a very practical feature of the new ELBO is that in GP regression it requires a minor modification to existing code. So a practitioner can easily train a GP with the new ELBO \nwhile keeping the computational cost the same as in the previous SVGP bound. \n \n> Tighter sparse variational Gaussian processes, Bui et al (2025), Under review TMLR.\n\nThank you for pointing to concurrent work. We will cite this work in the next version of our paper.  \n\n> Table 2 is confusing. The top methods are baselines but have nothing to do with your method I assume. The bottom methods both have ours but isn't your method only SVGP-new?\n\nWe agree with the reviewer.  We will modify the table to keep \"ours\" only for SVGP-new as suggested by the reviewer. The upper part of table gives two strong baselines from the literature (discussed in the first paragraph in the Related Work) that are based on different-type of extensions of the SVGP that allow to increase the number of inducing points. Indeed, they are unrelated to our method which replaces\n$p(f|u)$ by $q(f|u)$, but we have included them for comparison reasons. \n\n> When moving to the minibatch setting, the paper is essentially minibatching a KL term for the $q(f|u)$\n. How does this work exactly and how do you scale all the terms.\n\nStarting from the following expression for the ELBO (Equation 18 in the paper)\n$$ \n\\sum_{i=1}^N \\left(  E_{q(u)} [\\log \\mathcal{N}(y_i | k_{f_i u} K_{u u}^{-1} u, \\sigma^2 )]  - \\frac{1}{2}  \\log (1+\\frac{k_{ii} - q_{ii}}{\\sigma^2})  \\right)   - \\text{KL}[q(u) || p(u)],\n$$\nwe observe that each term in the sum (i.e., the full term inside the big brackets) depends only on a single data point $(x_i, y_i)$, which is what is needed for minibatch training. Based on this, we can obtain an unbiased ELBO (and gradient) using a minibatch of $b$ data as \n $$ \n\\frac{N}{b} \\sum_{i \\in minibatch} \\left(  E_{q(u)} [\\log \\mathcal{N}(y_i | k_{f_i u} K_{u u}^{-1} u, \\sigma^2 )]  - \\frac{1}{2}  \\log (1+\\frac{k_{ii} - q_{ii}}{\\sigma^2})  \\right)  - \\text{KL}[q(u) || p(u)],\n$$\nwhere note also that $E_{q(u)} [\\log \\mathcal{N}(y_i | k_{f_i u} K_{u u}^{-1} u, \\sigma^2 )]$ is analytic. We can include an appendix to describe the above details. \n\n> It’s unclear how the hyperparameters interact with the new approximation. Do they lead to better-calibrated predictive variances? Providing more details or empirical results on this aspect would help clarify the impact of the proposed method.\n\nAfter training with the new ELBO, we still use the previous standard SVGP predictive density as discussed in Section 3.1. \nSo the only difference is that the new ELBO can provide different hyperparameters (and inducing inputs $Z$).  Figure 1 and Table 1, 2 indicate that this can lead to better predictions in terms of test log-likelihoods. We can try to add an ablation to further study \nthe effect on the predictive variances.  \n\n| In the non-Gaussian case, the paper resorts to the standard SVGP prediction, as using the proposed approach directly would be computationally too expensive. This suggests that the main benefit of the method in this setting is primarily improved hyperparameter estimation rather than a fundamentally better approximation of the posterior. If so you should also mention some related work in this ara.\n\nJust to clarify here that for both Gaussian and non-Gaussian likelihoods we do predictions with previous SVGP predictive \nequations. We can clarify this further in the paper.", "claims_and_evidence": "Yes, the claims made in the submission are supported by clear and convincing evidence. The paper provides a well-structured theoretical foundation for the proposed variational approximation, including a derivation that justifies the tighter variational bound.", "methods_and_evaluation_criteria": "Yes, mainly. The proposed methods and evaluation criteria are appropriate for the problem at hand, and the experiments effectively demonstrate the benefits of the tighter variational bound. However, the datasets used in the evaluation seem somewhat simple, and the paper could have further strengthened its empirical validation by including a more complex benchmark, such as MNIST or a real-world dataset. This would provide additional evidence of the method’s scalability and effectiveness in practical applications.", "theoretical_claims": "I went over the proofs, and they appear to be correct. The derivations follow standard results from the literature, and the steps are well-structured and logically sound.", "experimental_designs_or_analyses": "NA", "relation_to_broader_scientific_literature": "Yes, the key contributions of the paper are well-situated within the broader scientific literature.", "essential_references_not_discussed": "There is another paper that was released at a similar time that essentially uses the same formulation and results in very similar findings. While I understand that this work can be considered concurrent, it would be beneficial for the authors to acknowledge and discuss this related research. Including a mention of this paper would provide readers with a clearer understanding of how the contributions of this work fit into the broader context and help differentiate any unique aspects of the approach.\n\nTighter sparse variational Gaussian processes, Bui et al (2025), Under review TMLR.", "other_strengths_and_weaknesses": "The main strength of the paper lies in its effort to improve the approximation of sparse variational Gaussian processes (SVGP), a crucial technique for scaling Gaussian processes. Additionally, the experimental results on benchmark datasets are compelling, providing strong empirical support for the proposed method. The paper is well-written, and the core idea is clearly explained.\n\nThe main weakness of the paper is that the results demonstrate only marginal improvements with the new bound across all evaluated tasks, without a clear example where the standard approach would fail without it. This makes it difficult to assess the practical necessity of the proposed refinement. Additionally, in the non-Gaussian likelihood setting, the paper opts for a simpler approximation rather than leveraging the more expressive variational form, seemingly due to computational constraints.", "other_comments_or_suggestions": "- Table 2 is confusing. The top methods are baselines but have nothing to do with your method I assume. The bottom methods both have ours but isn't your method only SVGP-new?", "questions_for_authors": "- When moving to the minibatch setting, the paper is essentially minibatching a KL term for the $q(f|u)$. How does this work exactly and how do you scale all the terms.\n- It’s unclear how the hyperparameters interact with the new approximation. Do they lead to better-calibrated predictive variances? Providing more details or empirical results on this aspect would help clarify the impact of the proposed method.\n- In the non-Gaussian case, the paper resorts to the standard SVGP prediction, as using the proposed approach directly would be computationally too expensive. This suggests that the main benefit of the method in this setting is primarily improved hyperparameter estimation rather than a fundamentally better approximation of the posterior. If so you should also mention some related work in this ara.", "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": "", "extra_scores": {}}, {"review_id": "kLE2i68u8p", "reviewer": "Reviewer_1vBe", "summary": "The authors introduce a tighter ELBO bound for inducing-point-based Gaussian process regression à la SVGP (Titsias 2009) and its mini-batchable extension (Hensman et al. 2013). The main idea is to use a more flexible ansatz for $q({\\bf f} | {\\bf u})$ than the conditional prior $p({\\bf f} | {\\bf u})$, in particular by introducing $N$ variational parameters into the covariance of $q({\\bf f} | {\\bf u})$, where $N$ is the total number of training data points. This leads to a small change in the form of the ELBO while retaining the same computational cost as SVGP and retaining mini-batchability in the case of the tightened Hensman-et-al-like objective (i.e. in the case where $\\bf u$ is not integrated out optimally). The authors provide a careful discussion of how the new bound differs from previous alternative bounds. The authors also consider how this construction can be applied to the setting with a non-gaussian likelihood. Experiments suggest that the new method can lead to higher ELBOs and better predictive distributions in practice.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "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": "Thank you for the insightful comments. \n \n> I would say that the evidence for 2 and 3 is convincing, while the evidence for 1 is somewhat weak. If this is to be a central claim of the submission, it should be supported with more empirical evidence, in particular in the simulated data setting where \"bias\" has a particularly straightforward meaning.\n\n> The authors note the tendency of SVGP to overestimate the observation noise and provide hints that the new objective lessens this bias. It would be great to see a more extensive study of this particular question, especially on simulated data, and reporting more detailed results for other kernel hyperparameters (e.g. the learned kernel scale).\n\nWe provide  here the learned hyperparameters for the 1-D Snelson dataset: \n\n            | Noise Variance  | Scale/amplitude  | Lengthscale|\n\nExact GP    | 0.0715    | 0.712    | 0.597 | \n\nSVGP-new | 0.087     | 0.485      | 0.615|\n\nSVGP         | 0.108       | 0.331 |  0.617 |\n   \nwhich show that SVGP-new has less bias than SVGP in this example. We plan to add further simulated regression examples \nin the appendix by following Reviewer's suggestion.  \n \n\n> I would love to get more details on the nature of the learned vs. Can the authors report a histogram and/or summary statistics \n\nYes, we can add histograms for the learned $v_i$s for some of the GP regression datasets (e.g. Snelson, Pol, Bike, Elevators) in the appendix. We give here some summary statistics (min, median and max values) of the learned $v_i$s  for the 1-D Snelson example (from  Figure 1): \n\n|min |  median | max |\n\n|0.172 |  0.952 |  0.9998 |\n\nThis indicates that there are few $v_i$s having quite small values but most of them are close to one in this example. We will update the Appendix to add histograms that visualize the final values of $v_i$s. \n\nRegarding the scalar values of $v$ for the non-Gaussian Poisson regression runs, we report here the\nlearned values.  For the toy Poisson regression the learned value was around $v=0.675$. Interestingly, for the real Poisson example in NYBikes, the value gets very small below $0.01$.  Perhaps this suggests that in this example the expected Poisson log-likelihood reconstruction term has a very strong influence in ELBO optimization,  so that the optimization prefers to set $v$ very small (which allows to reduce the posterior variance in $q(f_i)$ since the term $v( k_{ii} - q_{ii})$ becomes small). If the number of inducing points becomes sufficiently large and each $k_{ii} - q_{ii}$ becomes small, then of course the learned $v$ will become close to 1.     \n\n| Have you considered adding a new $N$-dimensional \"mean shift\" \n\nThis is a great suggestion. Yes, we tried to do this but we were unable to find to way to add a drift vector $\\Delta$ that gives a computational efficient ELBO of $O(N M^2)$ cost.  Here, is a short derivation of this. In the KL divergence $KL[q(f|u) || p(f|u)]$  and since the $q(f|u)$ and $p(f|u)$ do not have the same mean anymore the term $\\frac{1}{2} \\Delta^\\top (K_{f f} - Q_{f f })^{-1} \\Delta$ appears which has $O(N^3)$ cost. If we try to add some preconditioning to $\\Delta$, e.g., $(K_{ff } - Q_{ff})^{1/2} \\Delta$ then the $O(N^3)$ remains (although this time will appear in the expected log-likelihood term). We can add an appendix about these derivations, since they could be useful for future research.", "claims_and_evidence": "The main claims made by the authors are that their method can:\n\n1. reduce bias when learning the hyperparameters of the kernel\n2. can lead to better predictive performance\n3. result in tighter ELBOs\n\nI would say that the evidence for 2 and 3 is convincing, while the evidence for 1 is somewhat weak. If this is to be a central claim of the submission, it should be supported with more empirical evidence, in particular in the simulated data setting where \"bias\" has a particularly straightforward meaning.", "methods_and_evaluation_criteria": "While the empirical evaluation is not extraordinarily extensive, it is enough to convince me of the basic soundness and advantages of the proposed approach. Nevertheless, I would love to see the experiments and/or reported results extended in a few directions:\n- The authors note the tendency of SVGP to overestimate the observation noise and provide hints that the new objective lessens this bias. It would be great to see a more extensive study of this particular question, especially on simulated data, and reporting more detailed results for other kernel hyperparameters (e.g. the learned kernel scale).\n- I would love to get more details on the nature of the learned $v$s. Can the authors report a histogram and/or summary statistics for a few sample cases? Similarly in the non-gaussian case where $v$ is a scalar, what values of $v$ are you finding in practice? It's perhaps somewhat surprising to me that a single additional degree of freedom in the variational distribution can have the moderately large effect we can read off from Figure 3 (though in truth without predictive NLLs and the like it's a bit hard to measure the magnitude of the improvement/effect).", "theoretical_claims": "I have not checked the derivations in detail, but they are all intuitive/non-surprising and as such I have no reason to doubt their correctness.", "experimental_designs_or_analyses": "The setup of the experiments seems to follow best practices and otherwise looks sound. (Though personally I would never use a fixed Adam learning rate of $0.01$ without decaying it at least to some extent, as this is generically unlikely to result in particularly well-optimized objectives).", "relation_to_broader_scientific_literature": "The discussion of related work and how this new bound fits into the existing methodological landscape is extensive, very clear, and generally excellent.", "essential_references_not_discussed": "N/A", "other_strengths_and_weaknesses": "I would like to again stress the great clarity of the manuscript and thank the authors for taking the time to write a very clear and well-argued manuscript.", "other_comments_or_suggestions": "- typo in abstract: \"hyperpaparameters\"\n- typo in sec 3.2: whiten => whitened\n- typo: \"intiliazed\" in Fig 1", "questions_for_authors": "Have you considered adding a new $N$-dimensional \"mean shift\" variational parameter to your ansatz for $q({\\bf f } | {\\bf u})$? In particular this would mean using a mean of $\\bf K_{fu} K^{-1}_{uu} u + \\Delta$ for ${\\bf \\Delta} \\in \\mathbb{R}^N$. This would seem to be a natural companion to your adjustment of the covariance matrix. And the exact form for $p({\\bf f}|{\\bf u}, {\\bf y})$ in Sec. 3 suggests that there should be some slackness in the ELBO resulting from using the prior conditional mean without any modification. Granted this new parameter might conceivably be difficult to optimize (say because its entangled with the variational mean for ${\\bf u}$), or such a modification may have little effect on ELBO tightness or the learned hyperparameters or the learned predictive distribution, but from a methodological point of view it seems a very natural question to ask, so its omission is unfortunate. Why not satisfy the reader's curiosity?", "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": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: The sparse variational approach of Titsias (2009) and its siblings has enabled computationally tractable inference and training of Gaussian process models, allowing stochastic data mini-batching, non-Gaussian likelihoods, and efficient inference in hierarchical and latent variable settings. One key assumption in this approach is that the factor p(f|u) is retained in the variational approximation, q(f|u) = p(f|u). This submission proposes a new structured q(f|u) that results in tighter variational bounds whilst maintaining tractability. The proposal is an interesting contribution to the sparse GP literature, simple to implement and yields marginal but consistent practical gains. The reviewers were very positive and in agreement that this is a clear accept.
Decision: Accept (spotlight poster)