Tighter sparse variational Gaussian processes
Authors: Thang D Bui
Matthew Ashman
Richard E. Turner
Venue: Accepted by TMLR ()

Abstract: Sparse variational Gaussian process (GP) approximations based on inducing points have become the de facto standard for scaling GPs to large datasets, owing to their theoretical elegance, computational efficiency, and ease of implementation. This paper introduces a provably tighter variational approximation by relaxing the standard assumption that the conditional approximate posterior given the inducing points must match that in the prior. The key innovation is to modify the conditional posterior to have smaller variances than that of the prior at the training points. We derive the collapsed bound for the regression case, describe how to use the proposed approximation in large data settings, and discuss its application to handle orthogonally structured inducing points and GP latent variable models. Extensive experiments on regression benchmarks, classification, and latent variable models demonstrate that the proposed approximation consistently matches or outperforms standard sparse variational GPs while maintaining the same computational cost.

Published in Transactions on Machine Learning Research (05/2025)
Tighter sparse variational Gaussian processes
Thang D. Bui
thang. bui eduu
School of Computing Australian National University
Matthew Ashman
m.ca3990com.coo uk
Deportment Engineering
University of Cambridge
Richard E. Turner
ret26@com.cc uk
Department of Engineering University of Cambridge
Reviewed OpenReview: https: //ppre/veeeerveeeeuum net/ forum? id= L33DSu3320
Abstract
Sparse variational Gaussian process (GP) approximations based on inducing points have become the de facto standard for scaling GPs to large datasets, owing to their theoretical elegance, computational efficiency, and ease of implementation. This paper introduces provably tighter variational approximation by relaxing the standard assumption that the conditional approximate posterior given the inducing points must match that in the prior. The key innovation is to modify the conditional posterior to have smaller variances than that of the prior at the training points. We derive the collapsed bound for the regression case, describe how use the proposed approximation in large data settings, and discuss its application handle orthogonally structured inducing points and GP latent variable models. Extensive experiments on regression benchmarks, classification, and latent variable models demonstrate that the proposed approximation consistently matches outperforms standard sparse variational GPs while maintaining the same computational cost.
1 Introduction
provide a powerful framework for modelling prob- ability distributions over functions, offering principled uncertainty quantification and ease of use. Their flexibility in encoding domain knowledge such as smoothness, peridocity, domain-specific structure has led to widespread adoption across scientific and engineering applications. Exact inference in GP models poses significant computational challenges, requiring ((N3) time and O(N2) space complexity for N obser- vations. A suite of approximations have been developed to address these limitations. Most notably, sparse variational Gaussian processes (SVGP: Titsias 2009 Hensman et al. 2015 Matthews et al. 2016) address the poor computational complexity through the use of an approximate posterior distribution parameterised by a small set of inducing points.
Gaussian processes (GPs) (Rasmussen & Williams 2006
The standard SVGP framework employs a structured variational approximation that factorises the posterior distribution over the unknown function into two components: q(f) p(fæu u)g(u). Here, p(fzulu) rep- resents the GP prior distribution conditioned on the function values at inducing locations z, u= fz). The second term, q(u), is modelled as multivariate Gaussian distribution. Improved variational approximations have been developed such as SOLVE-GP (hi et al. 2020 g(u).
which use more sophisticated distributions for
This paper introduces a novel approach to improving SVGP approximations by modifying the conditional GP prior distribution at observed inputs, rather than focusing solely on the inducing point distribution. For

---

Published in Transactions on Machine Learning Research (05/2025)
Gaussian likelihoods, our approach yields new and improved collapsed lower bound the log marginal likelihood that involves no additional variational paramelers. Furthermore, we show how the uncollapsed form of our bound facilitates the use of stochastic mini-batch optimisation and extends naturally to non- Gaussian likelihoods through a single additional variational paramcter. We demonstrate the versatility of our method by integrating it with SOLVE-GP and extending to sparse variational approximations in the GP latent variable model (GPIVM; Lawrence 2005; Damianou et al. 2016). Our results demonstrate that by targeting our improved lower bound, our approach consistently improves the predictive performance and log marginal likelihood estimates across range of regression, classification, and latent variable modelling tasks. An implementation is made available in this repository https://g' tiigthu...... com conm/thannbuu/tightte sparse_gP We note that the concurrent work of Titsias (2025 presents an identical approximation but offers a different exposition to this paper.
2 Background
This section provides a concise introduction to pseudo-point based sparse variational Gaussian processes (SVGP: Titssas 2009 Hensman al. 2015 Matthews et al. 2016 Bui et al., 2017). Consider GP regression
with Gaussian observation noise:
= ((()))=(((f;..k)) x NNy; f(x),
(1)
(2)
where E NxD and y RN are the training inputs and corresponding noisy outputs, f denotes the unknown function mapping from input to output, ky is the covariance function governed by hyperparameters %, and o² is the observation noise. These hyperparameters, denoted collectively as 8, can be found by maximising the log marginal likelihood:
N 2
((0) = log(2T) yKK(Kr
(Kir + +0 y -ー log K + o'Il.
(3)
where Kr is the covariance between training function values f fxx) This objective takes O(N3) compute and is thus computationally prohibitive for large N. To sidestep this, use an approximate posterior judiciously parameterised by small set of pseudo-points inducing points as follows:
q(S) = p(. sillffuvyii (tllff:n)(((ll))) (())()),
(4)
where f() RM and RMxD are the inducing outputs and imputs, respectively, and M < N. The conditional ((fu)) in the approximate posterior is chosen to match that in the prior, leading to the following variational objective,
Folg(u), 0) (oog p( p())p(lls.m) {(f)
log log f: pappppppppll plf ulppffttgal
9S
v(5)
Σ
(5)
(.)())))
showed that when the likelihood Gaussian, an analytic optimal form for g(u) can be found, ((u) x ((u)N(p (y;; Kr, Kidu, I²I,, and that a collapsed bound is also analytically available,
Titsias
2009)
(1) N
1 log Qu + o²Il 202 trace
log(2x) 2 (Qr + -'I)-"y
(6)
2
where Qi Kn K Kut and DØ Kit Qr. Crucially, the bound above be computed in O(NM²). The non-collapsed bound in eq.
(5) is amenable to non-Gaussian likelihoods and data mini-batch sottings
(see e.g., Hensman et al. 2015), further reducing the training computational complexity to O(BM² + M
where B is the mini- batch size. Due to this small complexity and the ease implementation, the above variational approach has arguably become the go-to sparse approximation in the GP literature. In this work, we will revisit its core assumption of matching prior and posterior conditionals and show that relaxing this assumption results in tighter and more performant approximation.
2

---

Published in Transactions on Machine Learning Research (05/2025)
3 A tighter variational approximation
The variational approximation in eq. 4 is chosen such that the conditional of (u) identically matches the prior conditional p(flu). Instoad, we propose using the following variational posterior,
a()) - p(fff ()
(7)
where qffu)) = (f:: Mn > 0. Note that the mean of the prior conditional u(f(u - N(f: K kudu; D#)) is retained in ((flu). The resulting variational bound is,
D
1/2 MD T2), M is a diagonal matrix, M = diag((m1, m.,..., mxl) and
( log loggllf+rulf u (f
Fe(a(u), 0,
plf+sall u]g (
9(s)
-Kll((u))lbuu] ((u))LLa(f(le (p( (f)))
=I
a (en)) (())q(f(z,) (u) log log
(8)
Due to the structure of the variational distribution, the middle term can be simplified to,
1 [ ((u))LL(((())))l(f(f = trace(M) + log M
1
N 1 Σ 1 log(mm) - mnl
2
Collapsed bound and optimal M In the regression case, similar to the Titsias' bound, we can obtain the optimal form for q(u) x Pyu) Kr.k u, I²I), leading to the following collapsed bound,
N
1 Σ
m, 1 log(mm)) + mn
Fa M)
(Qn + + o (2T( 2 (Qu
log Qu + 0²1
2
2
0²
Setting the partial derivatives of [0 M) wrt mn to O, we arrive at mn = dntoz and the following bound,
o²
Fi(0) 2
N (2g(2 2 y (Qu + orT)-'y log Qu 0
Σ
log
1+
dn 02
where dn is the n-th element in the diagonal of Dr, d, = kjasn - kuf.
Comparison to Titsias' bound When M is the identity matrix, that is mn Vn, the approximation eq. (7) becomes the Titsias' variational approximation in eq. 4 and the bound in Fs(0) becomes the Titsias' bound (0) in eq. 66. We note that that (a)0) is tighter than (0) due the inequality log(1+-2 < for all x > 1. Our solution improves upon the solution to the Titsias' bound by allowing the marginals of the conditional approximate posterior, ((I(z,))))) to have smaller variances than that of the conditional prior, since the optimal m, - < 1. Intuitively, this reduces the strength of the coupling between q(u) and q()), enabling v()) to more freely model the data whilst allowing q(f) to be close to the prior elsewhere.
It is also worth noting that the middle term of our bound is always non-positive. One might think that adding this term to the bound would give a poorer approximation, yet, the improvement in the expected log-likelihood (due to the smaller predictive variances at the training points see predictions below) can yield larger improvement to counteract.
The new bound can also handle data mini-batching, yielding unbiased est istttooor of the uncollapsed bound in eq. as follows,
Stochastic mini-batch settings
(ig(u),0,,)) - [LL +
Σ [1 + log(ms) -m.
2B
b
N B
Σ I ((w)g(((==)(1) log pog fxx)))
(10)
B ((s(ss)
3

---

Published in Transactions on Machine Learning Research (05/2025)
Con-Gaussian likelihoods and mn parameterisation One parameterise m,'s to satisfy their positive constraint and optimise them directly at the cost of having N extra parameters. However, the optimal form for m. in the Gaussian likelihood setting suggests more efficient parameterisation my = 3/(d, + B) with 3 > shared across all data points. We will use the latter parameterisation for all of our experiments.
Predictions
The predictive mean and variance of the predictive distribution at a test input 2, are
m, k.u Mu:
(11)
U, k.. kun k.uKKlSyKKit
(kir Q-rk) - Qr-),
(12)
where Va = D (1 (1- M)D 1/2 Note that (i) we can compute the predictive mean at the same cost as previous sparse approximations, and (ii) the prodictive variance a training point can be approximated by More generally, the variance al new input that is not training inducing input is expensive due to the presence of Dæ in the last term. One path to address this could be to approximate D# by its diagonal matrix or to use only subset of training points for this computation. However, we find that simply ignoring the last term at test time does not impact the predictive performance while substantially reducing the prediction cost (see section 6.2
U, =,,d, + /, K
Connections to existing bounds collapsed bound:
We can use the log-sum inequality to bound the last term of our
En (1 [x)
Σ log 1 1 02)
dn
trace K Qui)
S N log
N
1
(13)
No²
Thus looser collapsed bound can be obtained:
N 2 (2g( 2 (Qr -"I) 'y 2
1 log QQ 0'Il
N
traco(K Q+1)
Fo(
log 1+
No²
This bound was derived by Artemev et al. 2021 based on bounds of the quadratic and log-determinant terms in the exact log marginal likelihood. This is also tighter than the Titsias' bound, that is Fa Z (0) 2 (1)...
One also view the proposed variational approximation
instance of the sparse orthogonal approach
of Shi et al. (2020) in which there are two sets of inducing points and with f, m, := 0 and S. D MD However, this view does not suggest new insights or potential improvements. We will next discuss how to use the proposed variational approximation to improve the sparse orthogonal approach and in the latent variable settings.
4 Application to sparse orthogonal variational GPs
The sparse orthogonal approach (SOLVEGP) of Shi et al. (2020) can be viewed as a structured approximation with two sets of pseudo-points u and v,
g(S) = p(S=s u, f)u(flu v u), quu, v) (N(;; ma, (w; K K add u mp, Sw)
.V ( 0
mu KvnK
S Kv
S, Kvn Kvn K kmnk)
Su
ma
E 2, For non-negative numbers 01 0g. .g.....a, and by bn. En a 1o ae > (Ea, a) log
with equality itf
2 b,
an/b constant.

---

Published in Transactions on Machine Learning Research (05/2025)
where Sm) and Sm) are the mean and covariance variational parameters. This approximation brings computational benefits over naively using single set of pseudo-points with cardinality M = M + Mn while matching the latter's performance. We will show that the same trick used for sparse variational GPs GPs-relaxing the conditional matching assumption q(f(u,v) = p(flu, w)ー can improve SOLVEGP. Ky [ut, o 1/2 K - Kr uv. M is diagonal matrix, M diag((m) "2,....,,N). and mn > 0. The resulting variational bound is
particular, similar to SVGP, we will use {f u, v)
MD'
Ir' T/2) where D
log p u, 0 o)p(f u, v)p(u, v)p(uls,
Fe(g(u,,),0,M) M)
plfasmaoo f. u, e) () u, w v)
9(3)
-KLLg(())lp(())) [LL
1 + 2 [1 + (og(mm)) mml
Σ
n
vl, |u, v) log p(yn)s((N))) + +Σ (.v.((en))
(14)
where q(v) = NN;;; my So), p(v) = Nw; O, Cvv), and Cvv = Kvv Kvu Kuv. Note that the predictive distribution a training point can be approximated efficiently, ((f(x,)) ~N(f(xn) m,y Un) with
m& = kk.kkkkkuu + Cf. C
Un = mm - CJ. Cwy Cwf.)) + k. Kul K.ukku + + S.
where Cab = kub - kan Kuukuu Similar to SVGP, the predictive variance at new test point is expensive due to the dependence on all training points. However, similar to the tighter approximation in section we found that simply ignoring this difficult term works well practice.
5 Application to Bayesian GP latent variable models
Consider a GP latent variable model (GPLVM; Lawrence 2005 with Gaussian observation noise:
p(x) = N(x;0,11
p(fh) = GP(f:0,k.
flylf,z) = NWy: o(x)
Both the posterior p(fly) and marginal likelihood p(y) are intractable. Instead, we introduce an approximate posterior of the following form:
v(fxx) = q(v))(() (f ()(()) xf Nif: Kuu D 2 M()) /3)
where M(x) diag([m; (I1), (2(2)))(.)..)... my (xx))) I) and mz
Note that q(flu, x) depends on æ through Kin; D# and M(x), and that when M is the identity matrix, that is mn = 1, we obtain the variational approximation of Damianou et al. 2016) We can bound the log marginal likelihood
F(g(f.z),0) = (L(g(()))p(a))) - KLLg(w)llp(u)) + 2 Σ (1 log(m(xx)) - m man))))
Σ/ ((z.)
( (f 2 u) Dog((u (((ln((((w))) (xx)).
We can obtain the collapsed bound by noting that the optimal form for q(u) is given by
qu)) p(u) exp (log ((y K K- u, o o²T g(e))
5

---

Published in Transactions on Machine Learning Research (05/2025)
SGPR T-SGPR
SVGP T-SVGP
90-1
1.0
0.5
negatve lower bound
0000
7x
10a
-05
-10
-15
1.00
SGPR of- of10.126
T-SGPR: of /3000330 =0.11.15
-2.5
20
40
600
Iteration [1bbggs]
hteratio [adam]
Figure Left and middle: Optimisation traces for SGPR, T-SGPR, SVGP and T-SVGP on the Snelson dataset with inducing points. Right: Predictive means and uncertainties. The stronger shade is for noiseless predictions.
Note also that
J. u ((u)g(fl() logg(yl))) = J, q(v) log Ny:: Kruk o²1) - Σ
My )
20
Together with Jensen's inequality, we arrive at the collapsed bound
F(a(x)) = -L(( - Σ
ma (a(xn)d 1 log mn (xn) + mu((n)
20
glan)
+ log ( Clo u. D)sco ((ec))
Setting derivatives w.rtt. m, (x) to gives
o2 dn + 02
(ma (mn(ln))((a))
g(xn)
which is satisfied by m, (x,) T,toz mz(#n) not need to (approximately) integrate out In to find mn.
a²
0²
The former is easier to implement
dutoz/
(((n)
6 Experimental results
We validate the utility of the proposed variational posterior in suite of experimental settings. We switch the variational objective with the proposed approximation in each setting, keep all other configurations unchanged, and measure the two's predictive performance. Implementations based on GPytorch and GPflow are released in this repository https: /tttps://githb. www. htmm/hhannnu//iggheee.. http://ttttterrprrreeeeee
6.1 Toy 1-D regression
To build intuition about the proposed method's behaviour, we first evaluate it on 1-D regression problem used by [Snelson & Ghahramani (2005). We compare (i) Titsias" collapsed bound in eq. 5 [SGPR] with the proposed collapsed bound in eg. (T-SGPR]. and (ii) Titsias' uncollapsed bound [SVGP] with the
proposed uncollapsed bound in eq. T-SVGP). Figure | illustrates the optimisation trajectories of these methods and the final fits for both SGPR and T-SGPR using five inducing points. The final values for both

---

Published in Transactions on Machine Learning Research (05/2025)
leggoiracted N-48827, Dw20
seggundirected
N--0000,
N-45730,
D-27
0.766
-1.94
1.00
Test Log-ileeliheed
1.00
-4.00
-0.97
100
-0.00
0.00
-5.00
-0.00
-0.05
-2.00
-2.00
1.00
0.000
5ang
houseelectric
N::34874 D=3
N=-15335,,
N-583250.
N=-0492800 D=
0.650
1.22
0.435
Log-ilkelhhod
-1.90
1.4
-1.32
1.38
-0.750
0.670
1.38
-1.99
-0.03
134
130
souero
Tovue
750tUroP
sowrc
Nowsrr
whoro
M-512
M-1024
Tighter
Tighter
Figure 2: Test ogg-ikelihhood for various sparse approximations eight regression datasets and various numbers of pseudo-points. For SOLVEGP and T-SOIVEGP, M is evenly split for u and v. Higher is better. Best viewed in colour.
uncollapsed and collapsed versions of the proposed bound appear lighter than that of the Titsias' bound in practice. The learned hyperparameters reveal that T-SGPR prefers smaller observation noise (0.115) and larger kernel variance (0.107) compared to that of SGPR (0.126 and 0.087, respectively).
6.2 Efficient predictive variances
A key practical consideration is the computational cost predictive variance in eq. The exact compu- tation requires Di which scales poorly with the training set size. We evaluate simplified variant that omits the term that involves Dr, the last term in eq. (12), and compare it to the exact variance calculation across three small benchmark datasets: wine, solar, and pumadyn32nm. Table 1 presents detailed comparison between the exact and approximate versions. We can see a pattern across all datasets: the simplified variant consistently matches the full model's performance while offering substantial computational savings. For this reason, we will be using the simplified version for all remaining experiments. The improvements in predictive performance sections 6.3 and 6.4 are therefore solely due to better estimation of the hyperparameters and different q(u).
Dataset
N/D
eq. 12
RMSE
Log- lokeliiedd
Time (s)
w. last term
0.47 0.01
-0.66 0.01
0.15 0.00
wine
1599/11
last term
0.47 0.01
-0.66 0.01
0.03 = 0.00
last term
0.93 0.07
-1.57 0.20
0.07 + 0.00
solar
1066/10 wo. last term
0.93 = 0.07
-1.56 + 0.20
0.03 + 0.00
8192/32 w. last term
1.00 + 0.01
-1.42 1. 0.01
21.12 1 0.06
pumadyu32mm
wo. last. term
1.00 0.01
-1.42 0.01
0.00
Table 1: RMSE, lieelioo,,ood,, and run time for two variants of predictive variance computation.

---

Published in Transactions on Machine Learning Research (05/2025)
-0.2
-0.4
SVGP,
SOLVEGP, M-512+512
-0.6
T-VGPP, M-S12
TSOLVEGP, M-5124512
0.2
0.0
Test IMSE
0.15
50
100
150
200
250
300
Epoch
Figure 3: Variational bound and test performance for various approximations trained the kin40k dataset.
Best viewed in colour.
6.3 Large-scale regression benchmarks
We next compare lour methods, SVGP, T-SVGP, SOLVEGP, and the SOLVEGP variant in eq. 14 [T- SOLVEGP), across three inducing-point configurations (M = 256, 512, 1024), on eight medium to large regression datasots. The datasots range from 40K to 2M data points with varying input dimensionalities (Yang al. 2015). We use the Matern-3/2 kernel and repeat each experiment 10 times, each employing random train/test split. The comparison results shown in figs. 2 and g We note that (i) both T-SVGP and T-SOLVEGP consistently match or slightly outperform (on 8/8 datasets), significantly outperform (on 3/8 datasets) their base counterparts, (ii) the performance improvement is also consistent across varions inducing-poinl configurations, (iii) the improvements (on 3/8 datasets) are also consistent across training runs and iterations, as shown in fig. S for the kin40k dataset, and (iv) SVGP and T-SVGP (and similarly SOLVEGP and T-SOLVEGP) have almost identical run time so the improvements here not come at any cost. We also quantitatively compare the hyperparameters provided by these approximations, including the kernel variance, kernel lengthscale, and observation noise in figs. to 11 respectively. Similar the 1D example above, observe that the proposed approximations prefer smaller observation noises and larger kernel variances.
6.4 MNIST classification
To evaluate the performance of the proposed approximation on non-Gaussian likelihoods, we run an ex- periment the MNIST digit classification task with 256, 512, 1024, and 2048 inducing points, using the SVGP, T-SVGP, SOLVEGP, and T-SOLVEGP variational objectives. F shows that both the proposed approximations achieve substantial performance gains in all metrics compared to their base versions.
6.5 Comparison to exact GP regression
We further evaluate the proposed approximation by comparing it exact GP regression small dataset, following the sct-up in Titssas 2009). We the Boston housing datasct vary the number of inducing

---

Published in Transactions on Machine Learning Research (05/2025)
N=254
M-512
M-1024
S0-244
Tighder better
-...
Tigfter NOTE
-4.270
1.70
135
6.300
-242
5,055
6.30
1.73
-0.235
Powor
Pawse
Figure 4: Log marginal likelihood approximations and test performance the MNIST way classification task. Best viewed in colour.
points, and the collapsed bounds in eq. 6 [SGPR] and eq. (9) [T-SGPR]. The hyperparameters fixed to that of the exact GP regression model so that the only variational parameter is the inducing inputs. We use four metrics for this comparison: the variational bound, test error, test log-likelihood, and the KL divergence between the predictive distributions given by exact GP regression and that provided by SGPR or T-SGPR. The average results across five repeats are included in fig. 5 confirming that the proposed tighter approximation consistently gives predictive distributions closer the exact GPR predictions than that of SGPR. The test metrics that only use the marginal predictive distributions for T-SGPR are also better, consistent with the trend that, for fixed hyperparameters, tighter variational approximation is better.
(exact pred. || approx pred.)
-350
-375
15-
Train ELBO
-400
4425
10
450
Exact GPR
5
445
SGPR
T-SGPR
-500
-0.45
0.340
Test Log- likelihood
0.335
0.50
Test RMSE
0.330
0.325
-0.55
0.320
-0660
0.315
10a
10a
10a
10a
M
M
Figure 5: Comparing the Titsias' collapsed bound and the proposed collapsed bound to exact GP regression, with the model hyperparameters fixed to those of exact GP. Best viewed in colour.
In addition, we evaluate the impact of the hyperparameters on the bound and predictive performance when the variational approximation fixed. Specifically, we optimise the uncollapsed SVGP bound on the Boston

---

Published in Transactions on Machine Learning Research (05/2025)
housing dataset until convergence then subsequently switch to the tighter bound while keeping q(u) and the inducing imputs fixed. The results included in fig. G demonstrate that the hyperparameters given by the improved bound yield marginally better predictive performance.
-360
SVGP
0.35
-370
T-SVGP
Train ELBO
-380
Test RMSE
0.34
-390-
0.33
-400-
0.32
4410
10a
102
10a
102
M
M
Figure 6: Comparing the SVGP bound and the proposed bound while keeping the variational distribution fixed. Best viewed in colour.
7 Summary
We build upon the standard sparse variational Gaussian process (SVGP) approximate posterior distribution through simple modification to the conditional GP prior distribution at observed inputs. Using proposed posterior approxination, bounds to the log marginal likelihood, and an uncollapsed form which facilitates its application with non- Gaussian likelihoods and is compatible with stochastic mini-batch optimisation. Furthermore, we show how our approach can be used to improve non-standard SVGP posterior approximations, such as SOLVE-GP (Shi et al. 2020).
derive a collapsed bound which improves upon existing SVGP lower
Our empirical results demonstrate consistent improvements in both predictive performance and log marginal likelihood estimates across diverse applications, including regression, classification, and latent variables mod- The proposed posterior approximations can be easily applied to other settings such as deep
elling tasks.
GPs and convolutional GPs (Van der Wilk et al. 2017 Blomqvist et al. 2020 Sun et al. 2021 Bui et al. 2016 Salimbeni & Deisenroth 2017).
Acknowledgments
We thank Martin Jankowiak for pointing out parameterisation issue in earlier version of this paper. We also thank the TMLR editor and reviewers for the constructive foedback.
References
Artem Artemev, David R Burt, and Mark der Wilk. Tighter bounds the log marginal likelihood of Gaussian process regression using conjugate gradients. In International Conference Machine Learning, pp. 362 372, 2021.
Christopher M Bishop and Gwilym D James. Analysis of multiphase flows using dual-energy gamma densito- metry and neural networks. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Delectors and Associated Equipment, 327(2-3):580 593, 1993.
Kenneth Blomqvist, Samuel Kaski, and Markus Heinonen. Deep convolutional Gaussian processes. European Conference Machine Learning and Knowledge Discovery in Databases, pp. 582-597, 2020.
10

---

Published in Transactions on Machine Learning Research (05/2025)
Thang Bui, Daniel Hernandez-Lobate Jose Hernandez-Lobato, Yingzhen Li, and Richard Turner. Deep Gaussian processes for regression using approximate expectation propagation. In International Conference on Machine Learning, pp. 1472 1481, 2016.
Thang D. Bui, Josiah Yan, and Richard E. Turner. A unifying framework for Gaussian process pseudo-point approximations using power expectation propagation. Journal of Machine Learning Research, 18(104); 1-72, 2017.
Andreas C. Damianou, Michalis K. Titsias, and Neil D. Lawrence. Variational inference for latent variables and uncertain inputs in Gaussian processes. Journal of Machine Learning Research, 17(42):1-62 62, 2016.
James Hensman, Alexander Matthews, and Zoubin Ghahramani. Scalable variational Gaussian process classification. In International Conference Artificial Intelligence and Statistics, pp. 351 360, 2015.
Vidhi Lalchand, Aditya Ravuri, and Neil D. Lawrence. Generalised GPIVM with stochastic variational inference. In International Conference on Artificial Intelligence and Statistics, pp. 7841 7864, 2022.
Neil D. Lawrence. Probabilistic non- con-lincar principal component analysis with Gaussian process latent variable models. Journal of Machine Leurning Research, 6:1783-188 1816, 2005.
Alexander de G Matthews, James Hensman, Richard Turnrr, and Zoubin Ghahramani. On sparse varia- tional mothods and the Kullbac--eebler divergence between stochastic processes. In International Con- ference on Artificial Intelligence and Statistics, pp. 231 239, 2016.
Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian Processes for Machine Learning. The
MIT Press, 2006.
Hugh Salimbeni and Mare Deisenroth. Doubly stochastic variational inference for deep Gaussian processes. Advances in Newral Information. Processing Systems, 30, 2017.
Jiaxin Shi, Michalis Titsias, and Andriy Mnih. Sparse orthogonal variational inference for Gaussian processes. International Conference on Artificial Intelligence and Statistics, pp. 1932-1942, 2020.
Edward Snelson and Zoubin Ghahramani. Sparse Gaussian processes using pseudo-inputs. Advances in Newral Information Processing Systems, 18, 2005.
Shengyang Sun, Jiaxin Shi, Andrew Gordon Gordon Wilson, and Roger B Grosse. Scalable variational Gaussian processes via harmonic kernel decomposition. In International Conference on Machine Leorning, pp. 9955-9965, 2021.
Michalis Titsias. Variational learning of inducing variables in sparse Gaussian processes. In International Conference on Artificial Intelligence and Statistics, pp. 567-574, 2009.
Michalis K. Titsias. New bounds for sparse variational Gaussian processes, 2025.
Mark Van der Wilk, Carl Edward Rasmussen, and James Hensman. Convolutional Gaussian processes. Advances in Neural Information Processing Systems, 30, 2017.
Zichao Yang, Andrew Wilson, Alex Smola, and Le Song. A la Carte Learning Fast Kernels. In International Conference on Artificial Intelligence and Stalislics, pp. 1098 1106, 2015.
A An even tighter but expensive approximation
We consider a more general form for the conditional covariance of qf(l)) as follows:
a(f) p(f4s ( ())
<ffu) Krnkku C),
11

---

Published in Transactions on Machine Learning Research (05/2025)
Again, also obtain the optimal form for q(u) Km o leading to the following collapsed bound
1
1
Fo(() - log(2n (T) 2 (Qr (Qu o² - y log Qu o'1 -
- traace((o 2 D
log C DD
2
We can derive the optimal C, C-1 =D DFF to 2I and the bound becomes:
NV Fs(0) log(2
1
1
y (Q + 1
2 log QQ + + log + o D
2
log(2
1 2y (Qtr + 0 1 log QQ + o
Σ
log(1 + 2 (Dit))
2
2
2
2
where (X) is the n-th eigenvalue of X. The bound above is expensive the original log marginal
likelihood.
B Exploring alternative parameterisations for the conditional posterior
Instead of the general as above or the form considered in the main text = 2 MD we consider two other parameterisations that might allow efficient collapse/uun-colaasse un-loollppeed bounds and predictions. We first rewrite the uncollapsed bound and the predictive mean and variance here for clarity,
F Funcollapsed KKLq - ((u))L]g((fuulll(ffu)) + EL
m. kmnk Kin
v. = k,, k.n K Kulkno kuo + kon Ken SiSkKKKukkn Küdkn. - (k.t Q.F) C) - Qr.),
We first consider C D While this allows efficient exact predictive marginal distributions training points, the middle term in the bound is costly compute due to the presence of Dæ:
J (f
N traceeD " D log M +
2
Another special case of the parameterisation presented in the main text is C mD iee. single m is shared across all training points. This conveniently leads to tractable exact predictive variances at training points, Un = mdn + f k KKKKS Kutkuf... The middle term in the bound can be simplified to,
J. (W)K((f((l)llll(TIl))
N + log(m) - -m].
-
In the regression case, this leads to the optimal m 0² 02/(N~1 En dn + 02) and the following collapsod bound,
N F0) log(2r)) / y (Qi 1211--12 log Qu 0²1
N log 1+ En dn No²
(15)
2
This is identical to the bound derived by Artemev et al. (2021). As shown in the main text, this bound is looser than the collapsed bound in eq. (9 due to the Jensen's inequality log(1 + +En Cn/N) 2 N-1 Σn log(1 x.)).
C Additional experimental results
C.1 Large-scale regression benchmarks
In addition to the test log likelihood results presented in the main text, we also compare the approximations using the test root mean squared error (RMSE) in fig. 7 and the variational lower bound in fig. The
12

---

Published in Transactions on Machine Learning Research (05/2025)
proposed approximations (T-SVGP and T-SOLVEGP) tend to give smaller RMSEs and better log marginal likelihood approximations, indicating lighter bounds leading to better predictions.
We also quantitatively compare the lyperparameters provided by these approximations, including the kernel variance, kernel lengthscale, and observation noise in figs. 9 to respectively. Overall, the proposed approximations prefer smaller observation noises and larger kernel variances.
kineDk
protein
leggdinected N=44827,
kaggundrected N=63608
0,030
0125
0.21
6.22
3.7
0.002
0.223
3,1
0.20
0.32
2.62
0.000
0321
0.23
Idroad
song
buzz N=583250, D
housedlectric
434874 D=3
N=515345, D=90
N-2049260, D=
0001
0.54
0.000
0.52
0.79
0.200
0.092
0,50
0.797
0.25
0.259
0.533
0.99
0.25
0.053
0791
0355
0.057
saueco
SoUSSOP
Monerta
900
SOURO
Tapro
HOWE
M-256
M-512
M-1024
Tighter better
Tighter
Figure 7: Test root mean squared errors (RMSE) for various sparse approximations eight regression datasets and various numbers of ascudo-points For SOLVEGP and T-SOLVEGP, M is evenly split for u and v. Lower is better. Best, viewed in colour.
C.2 My parameterisations for non- Gaussan likelihoods
We next assess the difference between two parameterisations of may m, B/(n, + B) with 3 shared across all data points, and {m with each m, specific to training point,
simple binary classification datasct. The results are included in fig. 12 We observe that (i) the former parameterisation is tighter than vanilla SVGP and the latter parameterisation is even tighter, and (ii) the difference between the predictive probabilities given by these parameterisations is small. For small and medium datasets, it is thus potentially beneficial to have N additional parameters. For large datasets, the former parameterisation is preferred to keep the number of parameters manageable.
C.3 GPLVM on the oil flow dataset
Finally, we demonstrate the proposed method's applicability to latent variable models through experiments with Bayesian GPIVM on the oil flow dataset. The multi-phase oil flow dataset consists of 1000, 12- dimensional data points belonging to three classes which correspond to the different phases of oil flow in a pipeline (Bishop & James 1993). Figure 13 compares the standard variational BGPLVM (Damianou et al. 2016 Lalchand et al. 2022) [VBBGLLMM against the proposed approximation in section TV-BGPLVM averaged across five repeats. The optimisation trajectories show that TV-BGPLVM achieves a lower final
negative ELBO (roughly -5.5 versus -522), indicating a more accurate posterior approximation.
13

---

Published in Transactions on Machine Learning Research (05/2025)
kr-kk
protesh
keggefrected
legoundirectdd N-63600.
N=45730,
N-46827,,
0.66
-1.0
133
-0.20
0.33
-1.0
-0.25
0.00
3droad N-434974,D=3
5010
Buzz
houseelactric
D-50
N-583250,
-0.72
-0.74
-1211
132
-0.16
-1222
66.00
-0.78
-6.00
-0.20
-2255
-0.25
laurcr
Pae
soutoo
Tipoo
M-25 56
M-512
M-1024
Tighter is botter
Tighter is worse
Figure 8: Log marginal likelihood approximations (ELBO) for various sparse approximations on eight re- gression datasets and various numbers of pseudo-points. For SOLVEGP and T-SOIVEGP, M is evenly split for and Higher is better. Best viewed in colour.
protein N-45730,
kegadirected
logoundirected N-63608,
N-40000,
X-49827,
0.40
0.52
Kemel vanaacee
6.55
0.34
0.3
0.2
322
s0ng
D-30
N-583250,
N-20R9200,
0.055
1.4
0.085
6.1
0.2
Kanar
Fonannc
Doonse
M52
M-512
M-1024
Tighter gives higher values
Tighter gives lower values
Figure 9: Kernel variances for various sparse approximations on eight regression datasets and various numbers of pscudo-points. For SOLVEGP and T-SOLVEGP, M is evenly split for and v. Best viewed in colour.
14

---

Published in Transactions on Machine Learning Research (05/2025)
kinA0k N=40000,
prosein
leggdirected
vaguundrected N-63008, D=27
N=43827,
04151
0.032
0.08
3/5
0,0150
0.353
02152
0.0000
1A3
0.151
0.008
0.0000
9A1
0.0022
2.2
0.0100
3droad X=434874, D-3
DIP? N-583250, D-
hoursddecttrcc N-204920,, D-11
N-515345, D-90
0.52
0.042
0.64
0.77
0845
0.060
0.611
0.009
0.003
Sauaco
veaupor
Filese
saure
Eaqueror
Tensuroo
M=1024
Tighter gives higher
Figure 10: Observation noise variances for various sparse approximations on eight regression datasets and various numbers of pseudo-points. For SOLVEGP and T-SOLVEGP, M is evenly split for and Best viewed in colour.
kin40k
protein
leggdirected 458277,
veggundrected N-63608.
7.8
Kernel lengthscale
21-
7.6
3.5
3.4
2.2
3.3
2.0
1.2
3droad No-430874, Dm3
DUz2
housoslectric N=2045280, D=11
N=515340, De90
R=583500, Dr
6.20
16.5
28.2
15.5
18.3
24.5
0.20
578
174
Soutte
Tdvor
M=1024
Tighter gives higher
Tighter
Figure 11: Kernel lengthscales for various sparse approximations eight regression datasets and various numbers of pseudo-points. For SOLVEGP and T-SOLVEGP, M is evenly split for and Best viewed in colour.
15

---

Published in Transactions on Machine Learning Research (05/2025)
T-SMGP ELBO -15,38
Defference
ELB0: 35,43
BLBO: -3533
T-VVP T-SMPP with
15
14
1.3
23
8
6
65-
5.
63
-16
-1.
-1.0
-10
SOTA End0 E.RT C.0M 120L
6.3 0.2 04
3.6 3.0
1.0
Probublity
Figure 12: comparison mn parameterisations on a two-dimensional binary classification task. First: predictions when using the standard SVGP approximation. Second: predictions when using the proposed tighter SVGP bound with m, 3/dn, + B) with 3 shared across all data points. Third: predictions when using the proposed tighter bounds with N additional parameters {m,) J Fourth: the difference between the predictive probabilities given by the two mh parameterisations.
60
V-BGPLVM TV-BGPLVM
-35
50
-4.0
40
Negative ELBO
-45
30
5500
20
5.5
10
5000
7500
10000
12500
15000 17500 20000
0
2500
5000
7500
12500 15000 17500 20000 Iterations
10000
Optimisation traces for variational Bayesian GPIVM the oil flow dataset. Best viewed in colour.
Figure 13:
16


Reviews:
[{"review_id": "99gCEc5qMy", "reviewer": "Reviewer_c3F7", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The main strength of the paper is that it attempts to improve the approximation of sparse variational Gaussian processes (SVGP), which is a key method in scaling up GPs. The extra GPLVM and SOLVE-GP applications solidify the technique's applicability. Furthermore, the experiments are convincing in benchmark datasets. The idea is well explained, and the paper is well written.\n\nHowever, the motivation for the introduced conditional posterior $q(\\mathbf{f} \\mid \\mathbf{u})$  could be clarified. The paper would benefit from a more explicit exploration of the interplay between hyperparameters and the proposed approximation, perhaps through controlled experiments isolating the effects. Moreover, explicitly demonstrating the constant-parameter approximation's validity in the non-Gaussian likelihood setting and discussing recent related works would strengthen the presentation.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The paper revisits the popular variational approximation for sparse variational Gaussian processes (GPs). It introduces an adjusted variational approximation that replaces the conditional posterior (traditionally set equal to the prior) with a more flexible variational distribution. The new approximation leads to a tighter variational bound. Furthermore, the approach naturally integrates into stochastic mini-batch optimization, enabling scalability to large datasets. Experiments validate improvements across multiple tasks, including regression, classification, and latent-variable modelling. The paper also demonstrates how the proposed method could benefit existing algorithms such as SOLVE-GP and GPLVM, highlighting its broad applicability.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "## Requested Critical changes\n\n- The motivation for $q(\\mathbf{f} \\mid \\mathbf{u})$ is abrupt. The start of section three goes straight into defining the new approximation. It would help the reader understand what goes wrong with the normal approximation. Some simple toy examples contrasting the predictive variances here would again help similarly as in [1]. I suppose this is the purpose of Figure 1, but it introduces a full optimization of hyperparameters and variational approximation, making it hard to tell the exact contribution. You could reuse this example but zoom in on the predictive variances, hold hyperparameters, and inducing points constant.\n\n- Related to this, it's unclear how the hyperparameters interplay with the new approximation. The paper states 'solely due to better estimation of the hyperparameters and a different $q (\\mathbf{u})$, but which one can you quantify this statement? What is the interplay with hyperparameters? \n\n- Non-Gaussian likelihood approximation. The definition of $\\beta$ needs some more explaining and testing. You are assuming that we can use the exact Gaussian approximation in the non-Gaussian case to save on blowing up the parameter count. I am okay with that, but it would be great to see a test. We know from [2] that the likelihood approximation now looks similar to a heteroskedastic regression, so a constant approximation would potentially adjust the term accordingly. However, the more expressive form could improve the non-Gaussian approximation. It's not clear, I don't think.\n\n- I found the writing of the Predictions section confusing. Why can we ignore the last term? I think it comes back to the interplay of variational parameters and hyperparameters.\n\n- A recent paper [3], appears to have come out concurrently. It would be good to mention and contrast the differences between the two. \n\n- The $\\mathbf{c}$ that appears before equation (9) is not defined. I understand it's a constant, but it would be good to say that.\n'We not that $F_4(\\theta)$' the F is not in the right latex format with eqn (9).\n\n\n## General (non-critical) improvements\n\n- Could you have $q(\\mathbf{f} \\mid \\mathbf{u})$ be in a different color to contrast equation (7) with equation (4) to help the reader to distinguish the change quickly?\n\n- The $\\mathbf{c}$ that appears before equation (9) is not defined. I understand it's a constant, but it would be good to say that.\nIt would be interesting to plot the various terms in the new collapsed bound during optimization.\n\n- The GPLVM on the oil flow example seems redundant as it only plots the ELBO. Can you not visualize the model's result in some way to see a change?\n\n\n[1] Bauer et al, Understanding Probabilistic Sparse Gaussian Process Approximations (2016)\n\n[2] Adam et al, Dual parameterization of sparse variational Gaussian processes (2021).\n\n[3] Titsias, New Bounds for Sparse Variational Gaussian Processes  (2025).", "broader_impact_concerns": "None", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "qJ0PWf1CAv", "reviewer": "Reviewer_Zy61", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n- The proposals are well motivated and easy to follow. The authors also compare their approximation to two other bounds in the literature.\n- A small ablation study was added to justify the cheaper predictive variances, thus allowing experimental validation of claims that cannot be theoretically validated.\n\nWeaknesses:\n\n- The section on experiments for GPLVM feels incomplete. Only the oil flow dataset is used, and the evaluation metric is limited to the ELBO value. This contrasts with other experiments that use predictive performance evaluation. Compare this with the exposition in Damianou et al. (2016), which uses nearest neighbor prediction as a proxy for class separation and additionally tests GPLVMs in various tasks such as image reconstruction.\n- The authors could expand more on their proposal, which is currently split between the main text and appendices A and B. In particular, a discussion on $\\boldsymbol{M}$ controlling the singular values of the conditional kernel, perhaps including how it might affect the spectral decomposition of the kernel, would be valuable.\n- None of the metrics compare the quality of the approximation to the true distribution. Since the authors tested their method on small-scale datasets in section 6.2, it should be possible to replicate Titsias's (2009) analysis and compute the moment-matching KL between the approximate predictive distribution and the full GP predictive distribution, demonstrating that the tighter ELBO corresponds to a better approximation.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "Starting from the variational framework for sparse Gaussian processes, the authors propose a different parametrization of the variational distribution $q(\\\\boldsymbol{f}|\\\\boldsymbol{u})$ that leads to a tighter lower bound on the model evidence. They replace the standard variational distribution $q(\\\\boldsymbol{f}|\\\\boldsymbol{u}) = p(\\\\boldsymbol{f}|\\\\boldsymbol{u}) = \\\\mathcal{N}(\\\\boldsymbol{K}\\_{fu}\\\\boldsymbol{K}\\_u^{-1}\\\\boldsymbol{u}, \\\\boldsymbol{K}\\_f-\\\\boldsymbol{Q}\\_f)$, where $\\\\boldsymbol{Q}\\_f = \\\\boldsymbol{K}\\_{fu}\\\\boldsymbol{K}\\_u^{-1}\\\\boldsymbol{K}\\_{fu}^\\\\intercal$, with a more general one $q(\\\\boldsymbol{f}|\\\\boldsymbol{u}) = \\\\mathcal{N}(\\\\boldsymbol{K}\\_{fu}\\\\boldsymbol{K}\\_u^{-1}\\\\boldsymbol{u}, \\\\boldsymbol{C})$ with different parametrizations of $\\\\boldsymbol{C}$. The authors highlight an approximation where $\\\\boldsymbol{C}$ is has the same singular value decomposition to $\\\\boldsymbol{K}\\_f-\\\\boldsymbol{Q}\\_f$ except that the singular values multiplied by positive free parameters $\\\\{m\\_i\\\\}\\_{i=1}^N\\\\subset\\\\mathbb{R}^N$.\n\nHowever, the approximate posterior of the highlighted approximation still requires inverting matrices of rank $N$, which has been a roadblock for previous attempts at these reparametrizations. The authors propose effectively using the original prediction distribution as if $q(\\\\boldsymbol{f}|\\\\boldsymbol{u})=p(\\\\boldsymbol{f}|\\\\boldsymbol{u})$ while using the hyperparameter/variational parameter estimations from their method. Experiments on UCI datasets demonstrate this approach doesn't negatively impact predictive performance.\n\nThe authors extend their work beyond Gaussian process regression/classification to the Gaussian process latent variable model (GPLVM), a building block of the more expressive deep Gaussian process model, and to SOLVE-GP, another reparametrization of variational sparse GP.\n\nThese proposals are validated by proving that the ELBO of this new parametrization is tighter than standard approaches and by comparing predictive metrics and ELBO across methods. In most cases, the new method achieves either better or comparable ELBO/predictive metrics.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "- [Strengthen] It would be nice to have appendices A and B incorporated into section 3.\n- [Critical] I believe section 6.5 is currently weak and should either be improved or moved to the appendix. For improvement: the text doesn't specify how many repetitions were performed for the oil experiment — if only one was conducted, the experiment should be repeated multiple times as in section 6.3; As mentioned earlier, additional metrics should be used, following Damianou et al. (2016), since a lower ELBO alone doesn't clearly communicate the downstream benefit of this methodology; Including another dataset or task for GPLVMs would also better illustrate the advantages of this approximation.", "broader_impact_concerns": "I believe there are no obvious broader impact concerns with this paper.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "MFqFMX3r8N", "reviewer": "Reviewer_uPL1", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n- The paper is well-written and easy-to-follow.\n- The proposed tighter bound is a novel contribution, and provides a new, orthogonal perspective to improve upon existing variational GP methods. \n- The experiment validation is comprehensive and convincing, which shows consistent improvement by leveraging a tighter variational bound. \n\nWeakness:\n-  While the paper focuses on the predictive performance, there is a lack of discussion on the implication of such tighter bounds in other aspects, e.g. uncertainty quantification. For example, Bauer et al 2016 shows that various sparse GP methods can overestimate/underestimate the noise parameter. It remains unclear whether the introduced term $m_n = \\frac{\\sigma^2}{d_n+\\sigma^2}$, which appears to function as a shrinking operator, offers any interpolation effects in uncertainty quantification between the SVGP and some limiting class of GP approximation models. Expanding on this aspect could offer valuable insights for using GP in downstream applications beyond predictive performance. \n\nReference:\n- Bauer et al Understanding Probabilistic Sparse Gaussian Process Approximations. Neurips 2016.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper presents a sparse GP learning method based on the inducing point framework and formulates a tighter variational bound. This tighter bound is achieved by allowing the variational posterior at the training inputs q(f|u) to have weaker dependency on the inducing points q(u). Similar to standard inducing point methods, the paper a collapsed bound and a non-collapsed bound, where the former applies to Gaussian likelihoods and the latter is amenable to non-Gaussian likelihoods and mini-batch optimization. The paper also discusses the application of the proposed method to handle orthogonally structured inducing points (Shi et al) and GP latent variable models. \n\nExperiments on regression, classification, and latent variable models show that the proposed method were able to achieve a better training ELBO compared to standard inducing point methods, and match or outperform their test-time predictive performance.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "A discussion over the weakness point raised above would be appreciated.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper significantly improves one of the most popular scalable GP methods. Its experiments are thorough and broad, providing substantial evidence for this method. These comprehensive results should lead to the adoption of this method, provided that the authors follow through on pull requests to popular GP libraries.

The biggest concern with this paper was organization and presentation, as the original exposition of the method was not easy to follow. The authors worked with the reviewers to improve the flow and clarity of the paper, though they have not yet combined an appendix despite promising a reviewer that they would do so.
Decision: Accept (poster)