Sparse Gaussian Processes: Structured Approximations and Power-EP Revisited
Authors: Thang D Bui
Michalis Titsias
Venue: NeurIPS 2025 poster ()
Keywords: Gaussian processes
approximate inference
Primary area: probabilistic_methods

Abstract: Inducing-point-based sparse variational Gaussian processes have become the standard workhorse for scaling up GP models. Recent advances show that these methods can be improved by introducing a diagonal scaling matrix to the conditional posterior density given the inducing points. This paper first considers an extension that employs a block-diagonal structure for the scaling matrix, provably tightening the variational lower bound. We then revisit the unifying framework of sparse GPs based on Power Expectation Propagation (PEP) and show that it can leverage and benefit from the new structured approximate posteriors. Through extensive regression experiments, we show that the proposed block-diagonal approximation consistently performs similarly to or better than existing diagonal approximations while maintaining comparable computational costs. Furthermore, the new PEP framework with structured posteriors provides competitive performance across various power hyperparameter settings, offering practitioners flexible alternatives to standard variational approaches.

Sparse Gaussian Processes: Structured Approximations and Power-EP Revisited
Thang D. Bui
Michalis Titsias Google DeepMind mtitsias@google..oom com
Australian National University thang. bui@anu. edu.
Abstract
maducing-point--ee sparse variational Gaussian processes have become the stan- dard workhorse for scaling up GP models. Recent advances show that these meth- ods can be improved by introducing diagonal scaling matrix to the conditional posterior density given the inducing points. This paper first considers an extension that employs block-diagonal structure for the scaling matrix, provably tightening the variational lower bound. We then revisit the unifying framework of sparse GPs based on Power Expectation Propagation (PEP) and show that it can leverage and benefit from the new structured approximate posteriors. Through extensive regression experiments, we show that the proposed block -dickodia approximation consistently performs similarly to better than existing diagonal approximations while maintaining comparable computational costs. Furthermore, the new PEP framework with structured posteriors provides competitive performance across various power hyperparameter settings, offering practitioners flexible alternatives to standard variational approaches.
Introduction
Gaussian processes (GPs) provide principled framework for modelling functions that offer calibrated uncertainty and safeguard against overfitting, among many other benefits (see c.g., Raamussen & Williams 2006 However, their computational requirement, cubic in the number of training data N. is prohibitive for many practical applications. This bottleneck motivates the development of
plethora of scalable approximation methods (Quiñonero- Candell & Rasmussen, 20055 Liu et al.) 2020), with sparse variational methods using inducing points arguably the most popular (Titsias 2009
Hensman al. 2013).
The key idea behind sparse variational GPs (SVGPs) is to approximate the posterior process using small set of M < N inducing points, reducing the computational complexity to O(N or (MS) in the batch and stochastic settings, respectively. A key assumption in the standard SVGP approximation is the prior distribution of the non-inducing function values conditioned inducing points remains unchanged in the approximate posterior, that is, ((fzrll)) (u) = p(fzulu))
the
Titsias 2025); Bui et al. (2025) recently showed that relaxing this assumption yields provably tighter variational bounds. In particular, the key innovation slightly adjusting covariance of qf by diagonal scaling matrix M, leading to improved predictive performance while maintaining computational tractability. This approach has the original SVGP approach as a special case when achieve even better approximations by
M I. Such improvement begs the question: considering more expressive structures for while preserving efficient computation?
To this end, propose using block-diagonal structures for M and show that this choice leads to provably tighter variational bounds compared to existing diagonal approximations while maintaining the same computational complexity and ease of implementation. We then show that these structured approximations can also help with other inference schemes beyond variational inference. Specifically,
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

certain structural choices for M lead to tractable Power Expectation Propagation (PEP) updates and approximate log marginal likelihood. This greatly extends and improves over the unifying framework of Bui et al. (0077).
The remainder of this paper is organised as follows. Section reviews sparse variational GPs, recent advances in structured approximations, and the PEP framework for sparse GPs. Section 3 presents the proposed block-diagona variational approximation. Section Textends the existing PEP framework with various structured posteriors. Section evaluates the proposed methods on a suite of tasks. We then discuss related work in section and conclude with a discussion of future directions in section
Background
We first provide summary of inducing-poing sparse variational Gaussian processes (SVGP; Tit- stas, 2009; Hensman et al. [2013] 2015] Matthews et al., (016), a recently proposed tighter bound (Bui et al. 2025: Titsias (225). and a power-EP based approach (Bui et al., 2017). Consider the supervised learning setting with an unknown input-output mapping J. GP prior over this function p(fh)) = GP( ), and pointwise likelihood p(y) f. X, w) = II. p(yn (x),), w), where X E ERRD and y E RV are the training inputs and outputs, ky is the covariance function governed by hyperparameters y, and is the likelihood hyperparameters. In what follows, we will use 0 to denote these hyperparameters and, when clear, drop the dependence on 0 for brevity. Inference and learning in this model are computationally challenging for large-scale datasets due to the ((N5) complexity; thus, efficient approximations are required. Sparse variational methods parameterise an approximate posterior based on M inducing points, { RMD RM) J. with M < N, as follows,
q()) = = ul ())
(1)
where f = (f(x))... ... fxx))]. Note that the factorisation here mirrors that in the prior, p(f) = p(fzfuulf u)p(f (u)) where p(u) = N(wu; O, Kun), p(f(u)
Nff; Kr, u, D
Dæ = Kr Qu. Qr = Ku. K = k(X, X), Ktu k(X, z). and Kuu k(z, z). Note that we use I to denote the function and f to denote the function values at the training inputs. The resulting variational lower bound to the log marginal likelihood is
Fo -KL(g())ll(()) - / ((w))LLg((() (2) [P([/))] +
Σ/ alw)a( plgl(snl
When (ff()) p(flu)= = N(: Kra Kirkalu u, DRr,, the bound above becomes.
Fi()()),0) -K..((()))lp())
Σ/ ((u))((((n)))) log (f
commonly known the uncollapsed SVGP bound (Hensman al. 2015 Titsias 2009). This bound conveniently allows both (i) tractable computation (O(NM*) in the batch setting or O + M²)
where B is the batch size in the mini-batch setting] and (ii) tractably handling of non-Gaussian likelihoods using quadrature or Monte Carlo estimation for the expected log-likelihodd terms. For the Gaussian likelihood, the bound can be simplified to
(g.)()()).0) -KLg(u))lp(())) Σ
/ (()))ogg log K.uuu us 02)
dnn 20
(3)
where 0² is the observation noise and dnn = [Drlan. Furthermore, an optimal form for q(u) can be found, q(u) x pyu)N(y; KrrKalu, o² IN),, yielding the following analytic collapsed bound (Titsias 2009).
Fi,rc(0) = Mog:((y: 0, Qi I'IN)
20 Σ dnn-
The SVGP approach above has arguably been the most popular scalable GP approach in the literature. More recently, Bui et al. (2025): Titsias (2025) show that this approach can be improved by relaxing
the q(flu) = p(flu) (u) assumption. Specifically, when qffl (f: Kn,K,ll D MD
1/2 1/2
ft

---

where M = diag((m mx], , the uncollapsed and collapsed bounds in the regression case are:
(2(() (K) [p(e]] Σ
g(ve) log WYyni ky, uKiu u, u,) log
1
1 1+
dnn
]]
(5)
2
o²
(2) (00)) Qn I 1 = log 1+
dun 0²
(6)
Note that the optimal form for m. is m. = 0²2 + dnn) < 1; and eqs. eqs. (3) and (4) for fixed and q(u) since log (1 + dnn/d²) 0²) S dnn/
and 66 are tighter than
The posterior approximation in eq. ① can also be used in other deterministic inference strategies. For example, in the regression case, for q (f(lu) = p(flu). Bui al. (2017) showed that Power- Expectation Propagation (PEP) yields an analytic collapsed approximate marginal likelihood,
(o) = (3: O, Qn + aD + O'IN) - 2a Σ log 1
1-a a
dnn
1+ a
(7)
σ²
and a closed form q(u), ((u) x (y: K aD + INI,,, where a is the power
hyperparameter in PEP. This framework encompasses a multitude of approximations, such as the
SVGP approximation (as -> 0) and FITC (Snelson & Ghahramani 2005 et al. 2010) (a 1).
block block-diagonal structured variational approximation
We first consider the following posterior approximation:
{(f) f(f=f.ulf If. ())))))))(z))
q(flu) = N(f: Ka,Kmlu,
where we have not posited a form for the covariance C. Interestingly, this leads to the familiar optimal form for qu)). q(u) x p(u) KikKilu u. I'IN). The resulting collapsed bound is.
1 F(0)= = (g:: O, Qu + I'IN)= - trace((o (arcelloor I + D log |C-
+ N 2
2
Except for some special cases, the bound above is as expensive the original log marginal likelihood
to compute. Specifically, as shown in the background, C D MD 1/2 with M = I Tiisiaas
f 2
2009) or == admit tractability, and each move (from IN mI v. and from mI y to diag (m, n== () makes the bound tighter. It is thus natural to enquire what structure to encode in M to further improve the bound, retain tractable computation, and potentially improve predictive performance.
= mIv (Artemev et al. 2021 or M = diag( m,,, n= Titsias 2025
Bui et al. 2025)
We now consider one such structure, a block diagonal M, M = blkdiag( (m, } where is the number blocks and m, No xN, Substituting this into the bound above gives
(.) log (y: O, Qi + -o IN) Σ
trace m, Dr, + trace mm.]
log |m|
No
We obtain the optimal my, m, (Is + Dr..
leading to the following collapsed bound,
Farre ( log (y: O. Qu + d'IN) E log Ia -D,
(8)
Due to the Hadamard's inequality, [L, + - Dr. < [1, +
(.).) and thus log [
DDr, < E, log(1 + o (D,,,.)).. In other words, the bound eq. 8 [M is block-diagonal] is provably tighter than the bound in eq. (6) [M is diagonal].
Similar to the standard SVGP approach. for large datasets, it is more convenient to work with the following uncollapsed bound that supports stochastic optimisation.
B Σ
1 [ q(u) log Kr.ukuu o ITI)) - log [I, - (9)
If the blocks are of roughly equal size, computing the bound in eq. ( using the entire training set takes O(M3 2 + B[N/B(3). However, in practice, perform stochastic optimisation,
3

---

M
Amale
bildding (mw)
B N
B.C. N
B.C N
B
BCN
B
N
B N
a
PTCC
Eq 14
Eq Intractale
Ittactable intaaliinttttaaoo
Ouinomeroo
8
PEP Intractabi
PEP
PEP Intractable
PEP Introctable
Buin,
Eq
Eq
intractblle
Titstos, 2009
Titslas 2025, Buis, 2025; Artameve, 2021
Titias, 2025, Buis, 2025
Equ 8,
Figure I: Connections between the sparse GP regression methods from the Power-EP perspective. Green means computationally tractable, red means intractable, and blue represents the new methods presented in this paper. B < N means the training points are partitioned into B disjoint blocks. B N means having the same number of blocks as training points, i.c., block size equal to
where we unbiasedly approximate the sum blocks in eq. p using random block to obtain the stochastic bound
Fa.-()) = (LL[([))ll[[[ + / g((u) log (y K Ku u o I.) log [T, + D
(10)
based which we perform stochastic gradient updates by cycling over the B blocks. If judiciously choose the block size N B to be M (i.e., block size equals to the number of inducing points), the computational requirement per iteration is only O(MS). Therefore, eq. 10 has a small implementation overhead compared standard stochastic sparse GP objectives. The precise extra overhead involves taking the Cholesky decomposition of I, + DD.,, needed when computing the log- dogeeenermnnant regularisation term.
We will next consider a special case. When we let all m, matrices to be the same, m, = m, we arrive at the optimal m. (Ig + B- E, Dr and the resulting collapsed bound,
B Fo log Wy: 0. Qi + d'IN) log
Bo² D
(11)
2
Since the log-determinant a concave function on the cone of positive definite matrices, we can apply Jensen's 's inequality to show that the bound above is less tight compared to eq. (8) As the block size equals one, this becomes the spherical bound in Titsias 2025); Artemev et al. (0011)
A disadvantage of diagonal and block diagonal structures in M is the expensive predictive covariance, uf at test time. Bui al.
However, we can approximate it by reverting using q( ((flu) (2025) noted that this approximation does not degrade the performance compared to the expensive exact predictive distribution. In other words, practice, we only use the new structured posterior in training, and therefore, any improvement in predictive performance at test time will come from better q(u)) and hyperparameters.
4 A more general approximation based on Power Expectation Propagation
Although variational sparse GP approach has captured the spotlight in the sparse GP literature, Bui et al. (2017) showed various variants of PEP can be competitive better. We will now revisit the framework of Bui et al. (2017) and explore how can be improved by leveraging the recent innovation in structured posterior approximations (Titsias (2025); Bui et al. (2025) and section 3) originally developed in the variational inference setting. We first write down the joint density of the exact model and the approximate posterior,
p(f.y) = p(f=s w))(ff....ppfllppet II p(ynlfm)
(12)
n=1
B
ql) x p(54s,ulff 4)((( (
(w)
(13)
hm

---

where the N training points are partitioned into B disjoint blocks, and the factors (2) are assumed to be Gaussian. Instead of using {( [) = p(flu) as in Bui et al. (2017). we consider ((flu)
(f: Ku u: D of MD 1/22 where M = blkdiag( (m (mms | ). that is, the blocks in M match that of the likelihood partitions.
it
The PEP procedure (Minka, 2004) iteratively updates (s(u) by (i) first remove an a-fraction of (s(u) from ((f) to form the cavity distribution, o (f) = ((f)/t ig (u), (ii) incorporate an a-fraction of No the likelihood for the b-th block p fo) I (f)) to form the tilted distribution, i(f)
q g(()))) (2)).f)) (iii) find a new approximation q(f) that minimises (L(G(f))l((f) and (iv) adjust (a(u) based on the new posterior using (s(u)) (g())//((())4) or (s(u) < (1)"(()) (4)(((()))) (N)- These steps are repeated for all blocks until convergence. Readers might have noticed that step (iii) is daunting task involves moment matching for the entire Gaussian processes; however, due to
the structure of the approximate posterior ((.). it is sufficient to perform moment matching for the finite function values u (Bui al. 2017). In addition, this procedure returns an estimate of the log marginal likelihood that can be used for hyperparameter optimisation.
Mirroring the derivation in Bui et al. 2017), we can show the optimal form for to (u() has rank No- (foll, to (u) (K u: gb. vb)). In the regression case, go = y/ and v. a 2DD I The full derivation rather lengthy and is included in the appendix; however, one can verify that this is a stable fixed point of the procedure by noting that the a-fraction of (u) is identical to the contribution of pd (ybl f6) to the posterior at u. f dffo(foll.)) (o The optimal q(u) is thus q(u) ((w)N(B K ut ablkddiag( ((D))
+
1/ MD 1/2 Jos To ) + o²IN). Furthermore, for the regression case, we can further derive the approximate marginal likelihood,
Fo, M) log.N(;: O, Qa + ablkdiag (D
1/2
MD
1/2 If h6
( + ININ)
Σ
1
log
I, +
1D/3
MD 2²
bbb
log [T. + a(ms I.)] + log |ms]
2a
2a
2
We note that, for a general a and M, including diagonal and block-diagonal cases, the PEP proce- dure above as well the approximate marginal likelihood for the regression case is computationally intractable (i.e.. cubic in NN due to the need to find the Bloock)ddiagonal of D //D now discuss the tractable special cases.
1/2 We will
ff
Remark When M is diagonal or block diagonal, the approximate marginal likelihood and poste- rior approximation above are only tractable as -> O. Specifically, when M diag( {mm) 1). the objective becomes the variational bound of Titsias (2025); Bui et at. (2025). and when M = blkdiag( ( )=1)) the objective matches the variational bound in eq.
Remark 2 When M = mly, the approximate marginal likelihood and posterior approxima- tion are computationallyy tractable for all a In particular, the optimal q(u) is q(u) & p(u)N(y; Kr,Kuuu, mablkdiag(((Daa) 18 + +zIIN),, and the approximate marginal likelihood becomes,
Fo(0, m) = log Ngy, O, Qu + mablkdiag( (eriddiag(((DA,))) (E InIN)
1 - O
N
N
log
I,
log (1 a(m 1)) log(m). (14) 2
2a
x²
2a
In this special case, we note the following. First, as a sanity check, we can see that when m = 1, we recover the Power-EP approximate marginal likelihood of Bui et al. (2017):
1 -
Di.. σ²
Fo.m=1 (0) = log BoggN(rr O, QF + ablkdiag(DDp) IN)
log I, + a
2a
Second, when a = 1, the objective in eq. 14 becomes Fo.o=1 (0, m) mblkdiag(DPp + σ IN).. When m = 1. this becomes the PITC marginal likelihood and for block size equal to one it further reduces to FTTC (Quihonero-Candlla & Rasmussen 2005).
(By; 0, QF -
Third, as a + 0, we recover the spherical bound in (Bui et al.) 2025 2021 Only in this setting, we can derive the optimal m (1 N
Titsias
2025
Artemev et al.
En dn
5

---

Finally, inspired by the uncollapsed variational bound, we can optimise an uncollapsed version of eq. (1) that supports stochastic optimisation follows,
[ ((u) log Kr,ukuu maDre + =
(o.r( -KL(g()))l(()) + Σ
1 2cx Σ log I
am Dr.r
N
N 2 log(m).
log (1 + a(m 1)) 2a
That is, instead of running the PEP procedure, we can optimise the objective above to yield the same fixed point PEP. We attempt to visualise the connections between the methods, the special cases and the broader literature in fig. T
Experiments
Having described the block-diagonal structure in sparse variational GPs and revisited the unified work of Bui et al. (2017) light of the new approximate posteriors, we will detail the experiments to qualitatively investigate (i) if the proposed block-diagonal approximation in section yields better performance and, if yes, how, and (ii) whether having m * 1 benefits power expectation propagation in section the same way does variational inference. All experiments were done on either V V100 GPU or a MacBook laptop. We provide an implementation here https
wttp://c/./moagnb/ooiiiiiie sp
5.1 1-D regression and biases in hyperparameter estimation
We first illustrate the difference between the proposed and existing methods on simple ID regression
problem (Snelson & Ghahramani] 2005). In particular, we compare Titsias' collapsed bound in eq. (4)
[SGPR]. the bound of Titsias (2025); Bui et al. (2025) in eq. (6) [T-SGPRJ. the bound with block
diagonal M in eq. (8) with 10 and 20 blocks [20 and 10 data points per block, respectively, BT-SGPR),
the PEP approach of Bui et al. (211)) with = 0.5 [PEP], and the PEP approach in eq. (14) with B N and a 0.5 [T-PEP] We used 5 inducing points in this experiment. The key results are summarised in fig. 2 It can be observed that (i) the Block-diagonal approximation improves over the diagonal one this example, (ii) increasing the number of training points each block tightens the bound, (iii) the structured posterior approximation also helps in PEP, and (iv) hyperparameter optimisation using a more structured approximation tend to result a smaller noise variance and larger kernel variance. We note that the PEP approximate marginal likelihood is not guaranteed to be a lower bound and therefore optimising it can result pathological behaviours, for example, when 1, the noise variance can be severely underestimated (Bauer et al. 2016).
encore bleoo TAF wit
Ex30
05
1x00
GARR 6.500 STOOPI 201 NO 323,
6.11 1-0
1.22 OLS 133, (H0 423 100
-25
Deac: STw0.203 of 2000000
10/0 67
Figure Left: Variational bounds during training on the Snelson dataset. Middle and right: Predictive mean and intervals using various methods and the final hyperparameter values.
To further investigate point (iv), we picked a subset of the KIN40K dataset with 5,000 data points, and ran experiment to compare the sparse approximations with exact GP. For each method, we recorded in table I the exact or approximate log marginal likelihood, the predictive performance measured by root mean squared error (RMSE) and log likelihood (LLL), the noise standard deviation and the kernel lengthscales. Similar the observation in the Snelson dataset above, the noise estimate is smaller when moving from M Iy to increasingly more structured M, translating

---

to better predictions. The trend seems to be consistent across two numbers of inducing points. In addition, there is no notable difference in the lengthscales between PEP, T-PEP, and the structured variational approximations, however, these methods tend to leverage more dimensions than SGPR for M 256.
Table I: Exact/approximate marginal likelihoods, predictive performance, and lengthscales given by various methods on 5,000 samples from the KIN40K dataset.
256
M = 512
Method
Obj.
RMSE
LL
lengthscales
Obj.
RMSE
LL
0
lengthscales
Exact
-0.66
0.12
0.80
0.00
-0.66
0.12
0.80
0.00
SGPR
0.88
0.26
-0.14
0.30
0.66
0.22
0.02
0.25
T-SGPR
0.78
0.22
-0.06
0.26
0.51
0.18
0.11
0.21
BT-SGPR [50]
0.75
0.22
-0.05
0.25
0.50
0.18
0.12
0.20
BT-SGPR [10]
0.66
0.20
-0.03
0.23
0,44
0.17
0.13
0.19
PEP [0.5]
0.66
0.23
-0.02
0.22
0,42
0.20
0.14
0.19
T-PEP [0.5]
0.48
0.20
0.02
0.18
0.18
0.16
0.19
0.14
5.2 Block- structured variational approximation
We next ran an experiment to validate the utility of the proposed block-structured approximation in section on four real-world regression datasets] For each dataset and each inducing point configuration (M = 256 or M 512), we compare the uncollapsed variational bounds of Titsias
Hensmaan et al. (2015) leq. (3). SVGPJ, Titsias (2025); (Bui al. (2025) [eq. (5), T-SVGP], and
[BT.)VG)), corresponding to M IN. M diag([m, H and M = bikdiig [m,) bi=) ). respectively. We repeated the experiment 10 times, each using a random train/test split, batch size of 500 (also the block size), random partitioning of the training data into blocks, and 300 epochs for training. The average variational bound (ELBO) and test performance after training are shown in fig. 3 Similar to the earlier experiments, the benefit of the block-structured approximation is also clearly demonstrated here: it tightens the variational bound compared to that of the diagonal M and consistently yields comparable or better predictive performance. We note again that (i) the estimated observation noise tends to be smaller when employing the new bound (see the appendix), and (ii) there is minimal implementation overhead compared to Titsias (2009 2025); Bui et al. (2025) to result in these gains.
the proposed bound in eq.
94-250
M-532
M--55
K-512
H-:5b
R-5522
R--255
M--522
7as 335E
0340
0.10
5.514
00:44
8.04
8.19
9.00
5.90
0.28
0.54
6.35
1,00
1,000
5.35
1000
0.00
6.32
0400
5.00
6.30
-1000
0.000
-234
-2000
-0.0
-2000
8.00
5.10
-.3
-1.00
5400
584
C40
06.0
SVGP
RSVGP
B13VCP
SVSP
PSSGPP
Sr.ver
3V3P
ST-vaP
SVEP
RSVGP
on.syge
T.SKPP
0T:3VP
SYGP
7:5VPP
87 87.5vee
ar.sycP
SYGP
Tover
ersvor
SVOP
T-SYPP
STSWGP
M-256
M=512
Tigher is better
Titter Monee
Figure 3: Lower bounds (ELBO) and predictive performance of various variational methods with
M = Iy [SVGP], M diag( [m, M= [T-SVGP), and M = blkdiag( ((m,) [BT.SVGP]. (,
'ee used the splits available in this repository http::///tttaaboooooooooott//o/a
7

---

5.3 Power-EP with a structured approximate posterior [M = mIy1
As shown in section + the structured approximate posterior considered by Titsias (025) can be utilised in PEP and in the regression case, the approximate posterior and marginal likelihood are analytically available. To evaluate practical utility, we ran an experiment on five small regression datasets, comparing the PEP approach of Bui et al. (2017) [M = IN]] ] to the proposed approach in section 4 [M mly]. The typical performance across various inducing point configurations is shown in fig. with the full results included in the appendix. It is noticeable that the Power-EP scheme with m * 1 tends to outperform the corresponding setting when m = 1. To clucidate the trend, plot the difference between the performance of M In and M mIN in fig. We note that m * outperforms m = on all datasets in terms of RMSE, but log- tog-likelihood performance degrades when is closer to 1. These results suggest that for * 1, intermediate values such as 0.5 are most competitive in terms of both RMSE and LL, in line with recommendations from Bui let al. (2017) when m 1.
2-44
M--22
3m91
-44
940
22.00
0.20
Figure 4: Predictive performance of power expectation propagation with M = In and M mm on two UCI datasets. Results for other datasets are the appendix.
alifoll
concrete
energy
wine
housing
0.000
RMSE <D better
3.005
-0.023
-0.02
-0.001
-0 0050
-2.015
0.000
8,025
0..000
0.013
0125
0.025
-0.20
-0.04
0.550
3.099
0.0175
0.02
0.02
L.A.. 200 is better
0.03
0.20
0.01
0..2
-0.02
-0.01
-0.04
-0.05
-0.02
-0.05
-9.05
0.3
-0.03
2=0 5g
e-
gw0.2
0-0
o-
Zgn0 Spnd Rgv0.9
Figure 5: Difference in PEP performance between M = In and M = mIx on five UCI datasets.
6 Related work
The use of inducing points for sparse approximations in Gaussian processes has rich history, to name few approaches, sparse online GPs (Csató Opper 2002). DTC (Seeger al. 2003).
Carhonnro-Candder & Rasmussen 2005). The most notable was the variational approach Titsias (2009), who introduced a principled method for selecting inducing points by optimising a variational lower bound. Hensmmn et al.
FITC approximation (Snelson & Ghahramani] 2005), and PITC
(2013 2015) extended this approach to enable stochastic optimisation and non-Gaussian likelihoods, significantly broadening the applicability of sparse GPs to large datasets. Other work on inducing point methods have exploited Kronecker products (Wilson & Nickisch, (2055), nearest neighbour
structures (Tran et al. 2021
Wu al. 2022) and inter-domain inducing points (Lázaro-Gredilla &

---

Figueiras- Vidal 2009 Hensman et al. 2018). Also, recent theoretical work (Burt et al.
2020) studied
the approximation convergence with respect the number of inducing points.
Our work is most closely related to the recent advances by (Titsias) (2025); Bui et al. showed that relaxing the standard assumption with diagonal scaling matrices improves the variational bound. Our block-diagonal extension naturally builds upon this line of work, showing practical benefits. Similarly, our extension of the PEP framework builds directly on Bui et al. (2017), expanding their unifying perspective by incorporating structured posterior approximations. Note that our work distinct from the PITC approximation. PITC is derived from the prior modification perspective. where the prior is modified so that the blocks of function values are conditionally independent given the inducing points. [uii et al. (2017) showed that this is equivalent to EP when retaining the prior conditional in the approximate posterior, which differs from our proposed structured conditional distribution.
2025 who
key component in the sparse GP approximate posterior is q(u), and imposing additional structures for this object will likely lead improvement. For example. [Shi al. (2020) showed that q(u) can be parameterised by two sets of inducing points, orthogonal to each other, leading to better predictive performance at much lower compute cost compared to doubling up the inducing points in the standard SVGP approximation. This line of work is complementary to our work here, as it focuses on a different aspect of the posterior, and thus, the two approaches can be combined.
well-known pathology of variational sparse GP regression is the large estimated observation noise variance (Bauer al. 2016). It can be partially alleviated by changing the objective function (Jankowiak et al. (0119) or mixing separate schemes for learning and inference (Li et al. (023). Our work shows that principled structured variational approximations can also partly address this issue.
Summary
Approximation schemes using inducing points are the method of choice for scaling GP models to large datasets. We show that (i) these methods can be improved by introducing additional structures in the approximate posterior and (ii) these new structures can be applied to various inference strategies, including PEP and variational inference. The resulting methods show comparable or better predictive performance and smaller hyperparameter estimation biases in many standard regression tasks.
There are several potential future directions. First, we have assumed that the size of the data blocks in a dataset is the same and the data partitioning in the experiments was random, but these can be adjusted based the data characteristics, potentially tightening the variational objective further. Second, the power hyperparameter a in PEP can be made private per block; this will require an understanding of when variational EP might work best and how dynamically select Third, a full discussion for non-Gaussian likelihoods and models beyond GP regression (e.g., deep GPs, GP latent variable models) and how they benefit from structured approximations is promising exploratory direction.
Acknowledgments
We would like to thank the anonymous reviewers for the feedback. TDB would like to thank the National Computational Infrastructure (NCI Australia), an NCRIS enabled capability supported by the Australian Government, for the computing resources.
References
Artem Artemev. David R Burt, and Mark van der Wilk. Tighter bounds on the log marginal likelihood of Gaussian process regression using conjugate gradients. In International Conference on Machine Learning, pp. 362-372, 2021.
Matthias Bauer, Mark van der Wilk, and Carl Edward Rasmussen. Understanding probabilistic sparse Gaussian process approximations. Advances Neural Information Processing Systems, pp- 1533-1541, 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):172, 2017.
Thang D. Bui, Matthew Ashman, and Richard E. Turner. Tighter sparse variational Gaussian processes, 2025.
David R. Burt, Carl Edward Rasmussen, and Mark van der Wilk. Convergence of sparse variational inference in Gaussian processes regression. Journal of Machine Learning Research, 21(131):1-63, 2020.
Lehel Csató and Manfred Opper. Sparse on-line Gaussian processes. Neural Computation, 14(3): 641-668, 03 2002.
James Hensman, Nicolò Fusi, and Neil D. Lawrence. Gaussian processes for big data. In Conference Uncertainty in Artificial Intelligence, pp. 282-290, 2013.
James Hensman, Alexander Matthews, and Zoubin Ghahramani. Scalable variational Gaussian process classification. In International Conference on Artificial Intelligence and Statistics, pp. 351-360, 2015.
James Hensman, Nicolas Durrande, and Arno Solin. Variational fourier features for Gaussian processes. Journal of Machine Learning Research, 18(151):1-52, 2018.
Jose Miguel Hernández-Lobatte Yingzhen Li, Mark Rowland, Daniel Bui, Thang D. and Heráández- Lobato, and Richard Turner. Black-box alpha divergence minimization. In International Conference on Machine Learning, pp. 151 1511-,, 2016
Martin Jankowiak, Geoff Pleiss, and Jacob R Gardner. Sparse Gaussian process regression beyond variational inference. 2019.
Miguel
and Aníbal Figueiras- Viiuel Inter-domain Gaussian processes for sparse infer- ence using inducing features. Advances in Neural Information Processing Systems, volume 22, 2009.
Rui Li, ST John, and Arno Solin. Improving hyperparameter learning under approximate inference in Gaussian process models. In International Conference on Machine Learning. 2023.
Yingzhen Li, Jose Miguel Hemúndez-Lobato, and Richard E. Turner. Stochastic expectation propaga- tion. In Advances in Neural Information Processing Systems, pp. 2323-2331, 2015.
Haitao Liu, Yew Ong, Xiaobo Shen, and Jianfei Cai. When Gaussian process meets big data: A review scalable GPs. IEEE Transactions on Neural Networks and Learning Systems, pp. 1-19, 01 2020.
Alexander de Matthews, James Hensman, Richard Turner, and Zoubin Ghahramani. On sparse variational methods and the Kullback-Leibler divergence between stochastic processes, In International Conference on Artificial Intelligence and Statistics, pp. 231 231-239, 2016.
Thomas Minka. Power EP. Technical report, Microsoft Research, 2004.
Yuan Qi, Ahmed Abdel-Gawad, and Thomas P Minka. Sparse-posterior Gaussian processes for general likelihoods. In Conference on Uncertainty in Artificial Intelligence, pp. 450-457, 2010.
Joaquin Salihonero-Candee and Carl Edward Rasmussen. A unifying view of sparse approximate Gaussian process regression. Journal of Machine Learning Research, 6(65): 1939-1959, 2005.
Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian Processes for Machine Learning. The MIT Press, 2006.
Matthias W. Seeger, Christopher K. I. Williams, and Neil D. Lawrence. Fast forward selection to speed up sparse Gaussian process regression. In International Workshop on Artificial Intelligence and Statistics, pp. 254-261, 2003.
10

---

Jiaxin Shi, Michalis K. Titsias, and Andriy Mnih. Sparse orthogonal variational inference for Gaussian processes. In International Conference on Artificial Intelligence and Statistics, pp. 1932 1942. 2020.
Edward Snelson and Zoubin Ghahramani. Sparse Gaussian processes using pseudo-inputs. Advances in Neural Information Processing Systems, 18, 2005.
Michalis K. Titsias. Variational learning of inducing variables in sparse Gaussian processes. In International Conference Artificial Intelligence and Statistics, pp. 567-574, 2009.
Michalis K. Titsias. New bounds for sparse variational Gaussian processes, 2025.
Gia-Lac Tran, Dimitrios Milios, Pietro Michiardi, and Maurizio Filippone. Sparse within sparse
Gaussian processes using neighbor information. In International Conference on Machine Learning. pp. 10369-10378, 2021.
Andrew Wilson and Hannes Nickisch. Kernel interpolation for scalable structured Gaussian processes (K)SS-GP). In International Conference Machine Learning, pp. 1775-1784 2015.
Luhuan Wu, Geoff Pleiss, and John P Cunningham. Variational nearest neighbor Gaussian process. In International Conference on Machine Learning, pp. 24114-24130, 2022.

---

A
Full derivation of the block-diagonal variational bound
We start with a general posterior approximation of the form:
((f) = p(fff,ulf (f ((f()) - K C)
(15)
(16)
where we have not specified the form of the covariance matrix C. The variational lower bound to the log marginal likelihood
F(g, 0) = KL (T([(()))lp(( -
+ ((())(((()) logr(ulf) (17)
Setting the gradient wrt qrt (() to zeros gives, q(u) x p(u) expl. ((flu) log f)). In the regression case, p(y) f) = (w; IN)) and thus,
/ ((flu)logp(e)
/ N(f: C) (y: f. I'IN)
(18)
1
log Nyy: K Kulu, o I'I)
20 trace( (C).
(19)
The middle term in the bound can be simplified to,
f
1 ((u))Lio((Fw)))ll((f u)] = trace D C) + log D log
1
N
C|
(20)
2
Substituting this and the optimal g(u) back to the bound gives,
N
F= log.ggN (y:; O, Qu + +I - trace (D + IN) log D# +
log CI +
2
2
2
2
When C ff
MDD
1/2 the collapsed bound above becomes, if
F(0)= Noy: O, Qu +o'IN)) Σ
1 trace Srace|m,DDr Dr., + trace[ms] - log |ms - No.
We can find the gradient of the bound wrt m..
1
F
2 0 DE.. + + - m
(21)
0mB
Setting this to zero gives m,
(Is + DE,.) -1 and the resulting m-collapsed bound:
J = QF I'IN)= 1 Σ log [I -Danl
(22)
When the block size is 1, the above bounds become the bounds presented in Titsias 2025 (022)).
Bui et al.
We now consider a special case when we let all m, matrices to be the same, m = m. The gradient wrt m in this case is,
0 F
_
1
Σ
0 2 Deat + +I
m m1
(23)
8m
2
b
This leads to the optimal m bound,
(I, B 1-10-2 -2 Drar.) and the corresponding m-collapsed
=o O. QF + InIN)
B log I
1
D
(24)
Bo²
A special case is when the block size is only 1, we arrive at the spherical diagonal approximation M mI (Titsias) 2025 Artemev et al. 202T). Note that, since the log-determinant a concave function on the cone of positive definite matrices, we can apply Jensen's inequality to show that the bound above (when all m blocks are the same) is less tight compared to the bound when all blocks different.
12

---

B Power-EP posterior and approximate marginal likelihood
B.1 Power EP steps
Given data of N input-output pairs (I Un Ji we use M pseudo-points y locations z to approximate the exact posterior. Power-EP EP posits the following approximation to the joint:
v) p(Srs u)
II (h) = plSfs u (1)((u) II (a(u) s())
where we have partitioned the data into B disjoint blocks, b indexes blocks of data and tz(u) the approximate factors. Crucially, we employ structured conditional approximate posterior (ffl) = [f;; Krn u, D 3 The Power-EP procedure with power a iteratively updates the factors {10)6-1 as follows:
1. Deletion step: Compute the cavity distribution by removing a fraction a of one approximate factor:
(f) (1)
plSff u) (f g(u)
(f) x
(2)
=
(25)
where q(u) = p(u) II, (u() and q ( = q (1)
2. Projection step: First, compute the tilted distribution by incorporating a corresponding fraction of the true likelihood factor:
f() = q (f)p"" fi) = p(f=s ( (u)) (v.l.)
(26)
Second, project the tilted distribution onto the new approximate posterior using KL diver- gence:
q(f) + arg min KLLā())la(()) +(f)
(27)
Due to the structure of the approximate posterior, this minimisation is achieved when the moments at the ceeutosnppuss are matched: E (u)] (x(r(o())) where Φ(u) {u, ) are the sufficient statistics (Bui et al. 2017). In practice, this can be done by using the moment- matchngg shortcut involving the gradients of the log-normalising constant of the tilted distribution.
3. Update step: Compute the new fraction by dividing the new approximate posterior by the cavity:
g(f)
ne (u)
(28)
qvi(S)
The factor then is updated using (e(xx) = Li, (u) or with damping, (+(u) = 11, (u)
f.now (u))
B.2 Optimal factors
The factors parameterised follows,
to tu(u)= NN; sb. TiS. T2,) 2b exp(p Ti. Ti,s z² - T2
(29)
The posterior distribution over u is therefore q(u)= = NNw m. m. S), where
s s-1= = = + Σ T2.
(30)
S s- m = Σ Ti
(31)
13

---

Similarly, the cavity distribution over a q' - Nu;m (S)), where
s =Kuu
2, + (1 a2T2 S-1 aT2,
(32)
S/
m
1
Σ
T + (1-1))111) (1 = S - aT
(33)
The moments of the tilted distribution (and the new posterior) can be computed efficiently using the following shortcuts,
d log 2,
mm
V uf:
(34)
dm 8 V de log Z
V =
(35)
))
where 2, -
is the normaliser of the tilted distribution.
At convergence, the optimal form of T2, is rank-Nb, T2, = waV, ww where W; =
V/b-iv>' ub KuuKun, V, -d2 Vid Vib-ivVe bu
ub and d2
de log 2n. d(m "
In the regression case, at convergence, ls (+) (K u; yo, ar D /DD bb + I²I)). We can check this by computing the contribution of an fraction of the exact likelihood to the posterior q(u).
/ (ffoll))) NN N( (D 1/2
MD
1/2
IN'' fo. I
(36)
Ny(b: K u. D MD 1/2 lm + o I//a,
(37)
which is exactly a-fraction of the optimal factor listed above.
B.3 Power-EP approximate marginal likelihood
After convergence, Power EP provides approximate log marginal likelihood:
log ZPEP = log /
p p(f+suuf: u ()((f(u)p(u)
II (u)df
(38)
1 ((((u)) ((p(u))
Σ
log 2o + (()) -(((u))
1
log Zgo
(39)
a
a
where
log 2, log
q ((fol)gg ( (w.l))dd
(40)
log Z, log
(l-a())l)) (folu)dfadu
(41)
((q(w))
M
log(2m)
log V
マ m m
(42)
2
2
M
((p(w)) log(2m)
log Kuul
(43)
2
M (w)) log(2m)
1
log V |1
M. VTV V (1 'mm
(44)
2
In the regression case, following closely the steps in (Bui et al. approximate log marginal likelihood
2017), we can derive the closed-form
log ZPEP- = log. Ayy; 0, Qa + ablkdiag [
MD
1/2 bbl)))) + o I²IN)
f
if
bb
+ Σ
log
I, + a
D
f
1/2
MD 1/2
bb
log II + a(ms - 1b] + log m
1
2
2x
2
(45)
14

---

B.4 Extension to classification
Instead of working with individual factors, we can use the stochastic Power-EP parameterisation et al. 2015) ie., assuming contributions from all blocks to the posterior equal (u = t(u). In addition, instead of running stochastic Power-EP iteration, we can directly work with q(u) xx p(u)t (u) and optimise the Power-EP energy, also known the black-box a-divergence objective
Heernndde--oobar et al. 2016). We will explore this direction in future work.
C Additional experimental results
C.1 Experimental set-up
In addition the details in the main text, we provide additional information here. For all experiments involving the block-diagonal matrix M, we randomly partitioned the training data into blocks. In the Snelson, kin40k, and Power-EP experiments, we optimised the collapsed bound using the L-BFGS optimiser. In the block-diagonal experiments with medium-scale datasets, we used the Adam optimiser with learning rate of 0.005. To initialise the inducing point locations, we picked M random training inputs in the Snelson experiment, and employed k-means clustering for all other experiments. For the later datasets, we used the median distance between the data points initialise the lengthscales and set the initial observation noise variance to 0.1.
C.2 Snelson dataset
We compared several sparse variational GP variants, including SGPR, T-SGPR, and BT-SGPR, with M = 10 to exact GP regression, and the objective and hyperparameters collected during optimisation are included in fig. [6] We note that, by using structured approximations, (i) the variational bound that is provably tighter for fixed hyperparameters indeed is tighter in practice, and (ii) the observation noise variance (the kernel variance) is smaller (larger).
SCPS
10
ragitre log Moilhoos TOOM boans/ T4
124
T-SCPR
BTTISIPAA (20 sock)]
SPT 110 vollaa
Exaci GITI
2
2x
10
ltrratonn
dbserviiinn roke vaturne
6X
12-1
10
Itoration
Iteroon
Figure 6: Objectives and hyperparameters provided by sparse variational and exact methods.
15

---

C.3 KIN40K hyperparameters
We include the full results, including the standard errors, for the KIN40K experiment in table
C.4 Block- Block-diagonal structured variational approximation
In addition to the predictive performance metrics in the main text, we also recorded the estimated hyperparameters when using the new structured variational approximations. These results are included in fig. 7 and agree with observations in smaller datasets (Snelson and kin40k); kernel variance and observation noise variance tend to be larger and smaller, respectively, when using the improved bounds.
C.5 Power-EP
We include the full results for all five datasets considered in the main text in fig.
16

---

lengthscale les
Exact/appros Exxet/approximaai marginal likelihoods, predictive performance, and lengthscales given by various methods on 5,000 samples from the KIN40K dataset, including the standard errors across three repeats
1 0.000
0.252 + 0.001
0.206 0.001
0,201 + 0.001
0.186 +. 0.001
0.187 + 0.000
0.138 4: 0.001
0.001
0.796 1 0.004
0.022 + 0.002
0.1 1155100000 0.001
0.120 + 0.002
0.133 +. 0.002
0.140 + 0.002
0.190 +: 0.003
M=
L1
0.117 1 0.000
0.215 + 0.001
0.184 1. 0.000
0.181 + 0.000
0.173 +. 0.000
0.200 + 0.001
0.164 +: 0.000
RMSE
-0.666 1 0.00S
0.660 0.006
0.515 0.005
0.499 + 0.006
0.437 0.006
0.422 + 0.005
0.184 0.006
Obj.
lengthscales
1 0,000
0.299 + 0.001
0,259 1 0.001
0.250 + 0.001
0.227 + 0.001
0.225 + 0.001
0,182 4 0.001
0.001
0.79 1 0.004
-0.136 + 0.002
-0.057 1 0.002
-0.055 + 0,002
-0.032 +. 0,002
-0.015 + 0.002
0.024 + 0.002
= 2.6
LL
0.117 1 0.000
0.256 + 0,001
0.223 1 0.001
0.217 + 0.000
0.200 +. 0,001
0.235 1 0.001
0.200 +: 0,000
RMSE
0.656 1 0.005
0. 0.883 + 0.006
0.779 0.006
0.752 + 0.006
0.659 0.006
0.661 + 0.006
0.480 +: 0.005
Obj.
BT-SGPR B= 50
T-PEP 0.5
oEP = 0.5
Method
Exact
SGPR
T-SGPR
BT-SGPR
Table 2:
17

---

BrSS SESTOO
kegoundreeen
YSVGR GGP
SVOP
0.2
030
324
320
380
376
374
$30000
11.00
1000
11000
BT.SVO PI-SNOR
kegoundreccm
M-256
Tighter olves lower values
Figure 7: Estimated hypprrprraeteess using SVGP. T-SVGP and BT-SVGP on four UCI datasers.
SVOP
520
0.2
5.00
20.
200
104
T YSYOO P
SYOP
Tighter aivss higher values
6.40
6.3
6.30
8.33
15.22
14.75
14.25
13.72
1000
0009
vospelircoo
T-VGP
T.5
SVOP
0.99
0.00
$7.30
16.0
10.00
13.25
15.20
1000
0000
M=512
M-512
M=256
SVOP
0.4
040
0.30
L.30
1,09
1.90
1,00
L335
0442
1400
BT.SVO
BTSVVP
protest
M-2
S00P
SVSGP
0.30
1.09
2.00
3.39
40000
10405
SWGP
Mi5112
STOP
KSYOP
TS
SVOP
243
10:44
120
ET.SV6 STSYOP
kin40k
003233
02930
81

---

Avg. Rank
RE
PEP
UPE
4.0
35
3.0
200
PEP
Avrg Kank RMSSE
REP
PER
VF2
Figure 8: A comparison between M. = and M = ml for Power Expectation Propagation
4.5
3.0
2.5
3.0
2.
4.5
3.5
25
3.5
3.0
355
P3
PEP
M-128
LL T
PSP
PE
6
0.25
0.35
00.25
0.30
0.75
8000
00.30
40.00
0.00
0000
75
R6
M-128
ROYSI
19
VIE
0.37
0.38
0.25
0.34
40.200
1.000
0122
61000
000
0.000
0.000
0.445
30
0.30
0.34
Mi64
RED
LL
PER
0.55
0.55
0025
0.25
000
0.6
0.5
03
10.00
1000
000
1060
420
0.30
0.35
0.45
PEP
M-64
RASE
HEP
PSP
WE
0.40
0.30
0.30
0.30
0.30
0.22
0:20
0.22
455
450
455
450
-400
0.40
0.30
0.30
10P
PE
M=32
LL
VPE
50
50 40
30
0.5
0.
0.00
04
00
10 00
-0.00
0.90
M4 11
M--22
RED
PEP
VRE
0.45
0.2
GR2
033
0.28
628
0.26
0.20
0.15
0210
0.05
0.48
0.45
0.70
048
0450
0.2
0.34
0.2
19

---

NeurIPS Paper Checklist
The checklist is designed encourage best practices for responsible machine learning research, addressing issues of reproducibility, transparency, research ethics, and societal impact. Do not remove the checklist: The papers not including the checklist will be desk rejected. The checklist should follow the references and follow the (optional) supplemental material. The checklist does NOT count towards the page limit.
Please read the checklist guidelines carefully for information on how to answer these questions. For each question in the checklist:
You should answer [Yes], [No] or [NA].
[NA] means either that the question is Not Applicable for that particular paper or the
relevant information is Not Available.
Please provide short (1-2 sentence) justificatin right after your answer (even for NA).
The checklist answers are an integral part of your paper submission. They are visible to the reviewers, area chairs, senior area chairs, and ethics reviewers. You will be asked to also include (after eventual revisions) with the final version of your paper, and its final version will be published with the paper.
The reviewers of your paper will be asked to use the checklist as one of the factors in their evaluation.
While "[Yes] " is generally preferable to "[No]" ", is perfectly acceptable to answer "[No] "provided proper justification is given (c.g., "error bars are not reported because it would be too computationally expensive" or "we were unable to find the license for the dataset we used"). "[No] or "[NA] is not grounds for rejection. While the questions are phrased in binary way, we acknowledge that the true answer often more nuanced, so please just use your best judgment and write justification to elaborate. All supporting evidence can appear either in the main paper or the supplemental material, provided in appendix. If you answer [Yes] to a question, in the justification please point to the section(s) where related material for the question can be found.
general, answering
IMPORTANT, please:
Delete this instruction block, but keep the section heading "NeurIPS Paper Checklist",
Keep the checklist subsection headings, questions/answers and guidelines below.
Do not modify the questions and only the provided macros for your answers.
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The claims: a new biock-diagonal structured variational approximation and an improved Power-EP EP framework for sparse GPs. These are presented in sections 3 and and experiments are provided in section 5
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made in the paper and important assumptions and limitations. A No
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results be expected to generalize other settings.
is fine include aspirational goals as motivation long is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors? Answer: [Yes]
20

---

Justification:
We provided several limitations (and potential future directions) at the end of section
Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors encouraged to create separate "I "Limitations" section in their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, c.g., if the approach was only tested on few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, worse outcome might be that reviewers discover limitations that aren't 't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proof's
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification: The posterior approximation and bounds are provided in the main text, with the full derivations in the appendix.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but they appear in the supplemental material, the authors are encouraged provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless whether the code and data are provided or not)?
Question:
Answer: [Yes]
Justification: We provided all necessary information (number of epochs, batch size, block size, number of splits, source of data) to repeat the experiments. Code will be provided upon
acceptance.
21

---

Guidelines:
The answer NA means that the paper does not include experiments.
If the paper includes experiments, a No answer to this question will not perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) If the contribution is primarily a new algorithm, the paper should make it clear how to reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) the contribution is a new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g., with an open-sonnee open-source dataset or instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors are welcome deseribe the particular way they provide for reproducibility. In the case of closed-source models, it may be that access to the model is limited in some way (e.g., to registered users), but it should be possible for other researchers to have some path to reproducing verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: We used publicly available datasets, typically used in the sparse GP litera- ture. The implementation was built on GPflow, and released here https: //gi htt//ggtttttut thub. com/
Baangbbu/tighhtr__pprss.._
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurIPS code and data submission guidelines (https://nips.co cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (e.g., for new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurlPS code and data submission guidelines (https
/ //nips
ssionooiccy for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
22

---

At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code is permitted.
6. Experimental setting/details
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification:
We have provided sufficient information to replicate.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification:
We have reported all error bars in the main text and the appendix.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
should be clear whether the error bar is the standard deviation or the standard error of the mean.
It is OK to report 1-sigma error bars, but should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful to show in tables or figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures or tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information the puter resources (type of compute workers, memory, time of execution) needed to reproduce the experiments?
Question:
Answer: [Yes]
Justification:
We have added comment about the compute workers in section
Guidelines:
The answer NA means that the paper does not include experiments.
23

---

The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurlPS Code of Ethics stttps://neuripst http:/s....c..lll/ttcccc.uoceaa.
Answer: [Yes]
Justification: We, at all times, follow all general research ethics and the NeurIPS code of ethics.
Guidelines:
The answer NA means that the authors have reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e.g., if there is special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [NA]
Justification: The proposed methods will benefit probabilistic ML practitioness. We do not expect any near-term societal impacts of our work.
Guidelines:
The answer NA means that there no societal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and tied to particular applications, let alone deployments. However, if there is direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional unintentional) misuse of the technology.
there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML).
11. Safeguards
Question: Does the paper describe safeguards that have been put in place for responsible release of data models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?

---

Answer:
[NA]
Justification:
We did not use any models or datasets that require safeguards.
Guidelines:
The answer NA means that the paper poses no such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (c.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [NA]
Justification:
We used open-source data sets and software packages.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswi thcode. com datasets has curated licenses for some datasets. Their licensing guide can help determine the license of dataset.
For existing datasets that re-packaged, both the original license and the license of the derived asset (if has changed) should be provided.
If this information is not available online, the authors are encouraged to reach out to the asset's s creators.
New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [NA]
Justification:
No new assets.
Guidelines:
The answer NA means that the paper does not release new assets.
Researchers should communicate the details the dataset/codeemoodll as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create an anonymized URL or include an anonymized zip file.
14. Crowdsourcing and research with human subjects
25

---

Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given to participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
No crowdsourcing or human subjects involved.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Including this information in the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then as much detail as possible should be included in the main paper.
According to the NeurIPS Code of Ethies, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/review based the requirements of your country or institution) were obtained?
Question:
Answer: [NA]
Justification: No human subjects involved.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is an important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Answer: [No]
Justification:
We did not use LLMs for this submission.
Guidelines:
The answer NA means that the method development in this research does not involve LLMs as any important, original, or non-standard components.
Please refer to our L.MM policy futtps://neurd: ce. /c// for what should or should be described.
26


Reviews:
[{"review_id": "3P2uJcIAQB", "reviewer": "Reviewer_tnwE", "summary": "Sparse Gaussian processes are popular methods for making GP inference computationally tractable. These methods use a set of small number of function points to approximate GP functions. While they yield tractable inferences, this comes at a cost of expressiveness. The paper proposes an advanced technique to improve the expressiveness of the approximate method while maintaining comparable computational cost. It builds upon an earlier result Bui et al. 2025, a variational approximation method, by further relaxing its diagonal scaling matrix assumption to a block-diagonal matrix. The paper further shows how this approximation technique can also be applied to other approximate methods, such as PEP. Through 1D regression problems, the paper highlights the improvement over standard approximate methods.", "questions": "* Eqns 2,3 5,9 show integration but omit the integrand variables. \n\n* Line 255: “however, these methods tend to leverage more dimensions than” any insights why this is the case? Is this related to the expressiveness of the methods? It might be important to see how the block sizes play a role in this as well.", "limitations": "yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 4, "significance": 3, "originality": 3, "strengths_and_weaknesses": "**Strength**\n\n**S1**: The paper deals with improving the accuracy of a popular GP approximation method while still maintaining comparable computation cost. This has a broader utility for both GP community and other ML sub-areas. \n\n**S2**: The paper presents a recipe that generalizes and recovers prior improvements on sparse GP methods (Tistias 2025, Bui 2025, Artemov 2021) from the block-diagonal assumption, giving a clear picture of the proposed method. \n\n**S3**: The experiments show the improvement of using the proposed block diagonal assumption in terms of improved ELBO, and more faithful noise variances (relative to the exact method). \n\n\n**Weakness**\n\n**W1**: The experiments, focused on 1D regression tasks, make it limited to judging how well this could improve >=2D experiments. \n\n\n**W2**: While the authors claimed the method “performs similarly or better to diagonal approximation while maintaining comparable computational costs”, there is a lack of experiments to support it. Such experiments could show we can trade accuracy for computational cost (as this work and related works are all striking a good balance between the two). \n\n**W3**: While authors appropriately cited inspiration from Bui et al (2017) framework, I think there could be more discussion on the challenges that come with relaxing the assumption to a block diagonal setting.", "was_revised": true, "final_justification": "I raised a concern regarding the application of the method beyond 1D datasets. The authors pointed out the experiments they did on >1D datasets. Therefore, I withdraw my concern on this point.   \n\n\nAn additional concern I raised was on runtime comparison to see how the computational complexity translates to CPU time.  The authors responded with an example to support the computational advantage of their method.\n\nI also noted the issue on the overall accessibility of the paper to a broader audience, which is also raised by other reviewers. The authors have acknowledged this limitation and noted they will revise the related work section to improve readability.", "rebuttal": "We thank the reviewer for the constructive feedback. We would like to answer your questions and clarify a few misunderstandings.\n\n*Weakness 1*: Please note that the inputs in our regression tasks (except the Snelson dataset in section 5.1) have >1 dimensions. Specifically, the input dimensions are 8 for kin40k, 9 for protein, 20 for keggdirected, 27 for keggundirect, 5 for airfoil, 8 for concrete, 8 for energy, 11 for wine, and 13 for housing.\n\n*Weakness 2*: in lines 109-112, page 4, we noted the precise difference between the methods in terms of compute costs. M is typically small, so the overhead is small compared to the overall runtime. For example, on a 2021 Mac laptop, for the kin40k dataset, when M=256, SVGP takes ~15ms for one mini-batch of 500 data points, and the overhead for BT-SVGP is 1ms. We will report the full runtimes in the next iteration. \n\n*Weakness 3*: Could you please expand and point us to the specific section in the text?\n\n*Questions*:\n* Equations 2, 3, 5, 9: Thanks! We will fix these.\n* Line 255: Yes, this is related to the tightness of the bound and expressiveness of the variational approximation. We observe that optimising a tighter variational bound yields smaller observation noise, larger kernel variance (please see Figure 7 in the appendix) and more input dimensions being utilised. And thank you for the suggestion on the block-size investigation. We will include this in the next iteration.\n\nPlease let us know if you have further questions."}, {"review_id": "BJacet3fMy", "reviewer": "Reviewer_kA7g", "summary": "This paper focuses on developing a more structured sparse approximation for Gaussian processes by leveraging recent advances in tighter variational bounds (Bui et al., 2025; Titsias, 2025). Specifically, it replaces the diagonal scaling in the posterior covariance with a block-diagonal structure, resulting in a even tighter variational bound. Moreover, this paper extends power expectation propagation (PEP) under this improved bound, or more exactly new posterior approximation. Experiments show that the proposed structured approximation achieves a tighter bound and improved regression performance compared to existing approaches.", "questions": "The two methodological foundations for this paper are: \n1. **Bui et al., 2025:** tighter variational bound by relax the assumption of $q(f)=p(f_{\\neq \\boldsymbol{f}, \\boldsymbol{u}})p(\\boldsymbol{f}|\\boldsymbol{u})q(\\boldsymbol{u})$ to $q(f)=p(f_{\\neq \\boldsymbol{f}, \\boldsymbol{u}})q(\\boldsymbol{f}|\\boldsymbol{u})q(\\boldsymbol{u})$ where we set $q(\\boldsymbol{f}|\\boldsymbol{u})=\\mathcal{N}(\\boldsymbol{f}; \\boldsymbol{K_{fu}} \\boldsymbol{K_{uu}}^{-1}\\boldsymbol{u}, \\boldsymbol{D_{ff}}^{1/2}\\boldsymbol{M}\\boldsymbol{D_{ff}}^{\\top/2})$, $\\boldsymbol{M}$ is diagonal.\n\n2. **Bui et al., 2017:** PEP updating framework for $q^*(f)=p(f_{\\neq \\boldsymbol{f},\\boldsymbol{u}})p(\\boldsymbol{f}|\\boldsymbol{u})p(\\boldsymbol{}u)\\prod_{b}t_b(\\boldsymbol{u})$. \n\nIf my understanding is correct, the contribution of the current paper lies in extending the above by using a block-diagonal structure for $\\boldsymbol{M}$, instead of a diagonal one, and then computing both the corresponding variational lower bound and the PEP objective accordingly. \n\nHence my concerns are:\n\n1. The proposed use of a block-diagonal structure for the matrix $\\boldsymbol{M}$ closely resembles the PITC approximation. In PITC, the conditional distribution is modeled as \n$q(\\boldsymbol{f}|\\boldsymbol{u}) = \\prod_b \\mathcal{N}(\\boldsymbol{f}_b;\\boldsymbol{K}\\_{\\boldsymbol{f}_b}\\boldsymbol{K}\\_{\\boldsymbol{uu}}^{-1} \\boldsymbol{u}, \\boldsymbol{D}\\_{\\boldsymbol{f}_b,\\boldsymbol{f}_b})$, which leads to a block-diagonal structure in the covariance of $q(f)$. Hence, the current method seems more like an extension of PITC under the tighter bound/posterior approximation, which may weaken the originality of the proposed approach.\n\n2. The current framework appears more like an incremental extension, by replacing the old bound/posterior approximation in PEP with a new one, without offering deeper justification for why PEP is particularly suitable or necessary in this new setting, such as improved robustness, analytical tractability, or computational efficiency.\n\n3. The experimental section lacks some important details. For example, in the 1-D regression task, what kernel was used? How were the training and test sets split? The conclusions drawn from the middle and right plots of Figure 2 are not clearly explained. In Figure 3, what does “tighter is worse” refer to, and under what conditions? The conclusions could be refined to better support the claims.", "limitations": "yes, in future work", "rating": 4, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Although the proposed method is technically complex and requires substantial background knowledge, the overall presentation is smooth and coherent. Still, it may be difficult to follow for readers less familiar with this field. \n\nThis paper builds upon recent work on tighter variational bounds (Bui et al., 2025; Titsias, 2025), particularly the inference framework proposed in Bui et al., 2025. The proposed block-diagonal structured approximation closely resembles the Partially Independent Training Conditional (PITC) approximation. While the paper theoretically and  empirically verifies that combining this structured approximation with the tighter bound leads to improved variational bound and performance, the overall contribution appears somewhat incremental, as it primarily extends existing ideas with limited novel theoretical or methodological insights.", "was_revised": true, "final_justification": "I had initially questioned whether introducing a block-diagonal structure under the tighter bound framework of Bui et al. (2025) and Titsias (2025) brings a sufficient contribution, given that similar ideas have already appeared in PITC. However, since this work applies the structure at the posterior level, it can be seen as a new direction. However, as I am not fully in the scaling GPs community, I cannot fully assess the impact of this work on the field, and therefore I will lower my confidence.", "rebuttal": "We thank the reviewer for the questions and feedback. We will now answer the reviewer’s questions.\n\n*Concern 1*: We would like to clarify that our motivation/approach is very different from that of PITC. The PITC approximation is derived from the model/prior modification perspective, whilst our derivation is based on the posterior approximation view (see (Bui et al, 2017) for more discussion). Specifically, PITC modifies the prior so that in the prior, the blocks of function values are conditionally independence given $\\mathbf{u}$. Bui et al (2017) showed that this is equivalent to posterior approximation using EP when using multiple function values in a factor. This approach is distinct from what our proposal, we make the approximate posterior richer by utilising $q(\\mathbf{f}|\\mathbf{u})$ and a structured $\\mathbf{M}$ instead of using $p(\\mathbf{f}|\\mathbf{u})$ (see lines 132-134). \n\n*Concern 2*: We provided supporting evidence in section 5.3 that the new structured posterior approximation improves PEP, in the same way as it does to VFE, in a range of datasets and alphas. We also noted that intermediate alpha values tend to work best. These observations are empirical and, in our view, strong enough to justify the utility of the new approximation. And, generally, we can’t claim that PEP or VFE, or a flexible posterior is guaranteed to lead to a better performance, since the performance of each inference method is dataset and model dependent.\n\n*Concern 3*: Thanks for pointing these out - we will clarify them here and also in the text. \n* kernel: we used the Matern-3/2 kernel. \n* train/test splits: we used the same random train/test spits used in previous benchmark studies such as Yang et al (2015, A la Carte - Learning Fast Kernels) or Wilson et al (2016, deep kernel learning). The link to the splits is provided in the footnote on page 7.\n* middle and right plots in figure 2: The key message we want to get crossed is “hyperparameter optimisation using a more structured approximation tend to result in a smaller noise variance and a larger kernel variance”, as in lines 192-194.\n* figure 3: “tighter” means a method that gives a provably tighter variational bound (for the same hyperparameters), “better” means better lower test RMSE or higher test log-likelihood or lower train lower bound. In this figure, we want to highlight whether moving to using a tighter bound, from SVGP to T-SVGP or from T-SVGP to BT-SVGP, yields a better predictive performance. \n\nPlease let us know if you have further questions."}, {"review_id": "uM0JzbM06w", "reviewer": "Reviewer_sm4Y", "summary": "The applicability of Gaussian process models is gated by their costly inference, and the go-to solution to GP inference is inducing point sparse variational Gaussian processes (SVGP). SVGP replaces the standard Gaussian process posterior \n$$ p(f\\vert\\mathbf{y}) \\propto \\prod_{n=1}^N p(y_n \\vert f_n) p(f_{\\ne \\mathbf{f}, \\mathbf{u}}\\vert \\mathbf{f}, \\mathbf{u})p(\\mathbf{f}\\vert \\mathbf{u})p(\\mathbf{u})$$\nwith the variational posterior \n$$q(f) = p(f_{\\ne \\mathbf{f}, \\mathbf{u}}\\vert \\mathbf{f}, \\mathbf{u}) q(\\mathbf{f}\\vert \\mathbf{u}) q(\\mathbf{u}),$$\nand minimizes the KL-divergence. \nThe common wisdom of SVGP replaces $q(\\mathbf{f}\\vert\\mathbf{u})$ with $p(\\mathbf{f}\\vert\\mathbf{u})$, but Bui et al. (2025) and Titsias (2025) suggests that the variational lower bound can be tightened by an additional parameter $\\mathbf{M} \\in \\mathbb{R}^{N\\times N}_{\\succeq 0}$. Instead of the diagonal scaling matrix $\\mathbf{M}$, this paper proposes a block-diagonal structure that tightens the variational lower bound (VLB) without introducing a significant computational overhead, and illustrates theoretically and empirically that the block diagonal structure performs similarly or better to their diagonal matrix alternatives. The paper also draws a connection between SVGP and Power Expectation Propagation, and sees a similar boost in flexibility when a block-diagonal scaling matrix is introduced.", "questions": "- I have some trouble understanding the big picture, especially the part about replacing $q(\\mathbf{f}\\vert \\mathbf{u})$ with $p(\\mathbf{f}\\vert \\mathbf{u})$ at prediction. Does the improvement in performance for BT-SGPR mainly comes from better hyperparameter optimization, as the gap between log-likelihood and variational lower bound decreases with the better flexibility of the variational approximation? If we use SGPR with the same hyperparameter value, would the performance be the same?  \n- I don't understand the statement in line 87-88. How does the flexible covariance matrix $\\mathbf{C}$ yield an optimal $q(\\mathbf{u})$?\n- I understand the variational lower bound formula (e.g., eqs 1-4) are functions related to the variational posterior $q$ and the hyperparameter $\\boldsymbol{\\theta}$ , but it would be good to state explicitly the arguments to the lower bound functions. \n- The block-diagonal structure pre-supposes a partition of data in batches. How would the data be partitioned?", "limitations": "Yes. I think the authors have adequately addressed the limitation of this paper, as in, the paper's methodological nature has little negative societal impact.", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "## Strengths\n- The paper is well-written and theoretically sound, as the block-diagonal matrix structure generalizes the diagonal scaling matrix, and its added flexibility creates a tighter variational bound. The paper also has a good narrative pace that links a collection of variational lower bounds by emphasizing their connection. \n- I am not familiar with the work on power expectation propagation, but the paper's discussion on this topic seems novel and might be of interest to further development in the approximate inference of GPs. \n- The paper showcases convincingly that the added flexibility is helpful in obtaining better approximate posteriors and parameter inference via empirical assessments. \n\n## Weaknesses\n1. The methodology of introducing an added scaling matrix provides an appealing way to tighten the VLB is an interesting approach, but I am not sure of its applicability: much effort is devoted to learning a flexible parametric form for $q(\\mathbf{f}\\vert\\mathbf{u})$, but it is hard to make use of that as part of the variational predictive distribution (lines 118-123). I understand this drawback applies to all methods that investigates the scaling matrix parametrization, but this paper does not explain further on this point. For example, the second panel of Figure 2 illustrates different approximate posteriors, yet it is hard to see meaningful differences between them. \n2. The methodology seems incremental. Replacing the diagonal matrix with a block-diagonal matrix is a natural next step, and it naturally leads to better flexibility. The paper does not discuss the difference in the dimensionality of variational parameters, or its effect on computational complexity. Therefore, if increasing the expressivity only leads to a marginal improvement at the cost of higher computational expense, it would seem that the method has more limited applicability. \n3. I think the paper and supplementary materials do not provide adequate details (see the _formatting concerns_ section).", "was_revised": true, "final_justification": "I have interacted with the authors regarding my review, and it has been helpful for me to better understand the paper. \n\nThe main concern from my previous assessment is that the paper's contribution seems incremental, especially due to the fact that the random partition of training data introduces a limited performance improvement, and a random partition seems like a suboptimal way to exploit the block-diagonal structure. I am mostly satisfied with the author's promise to look further into this issue and amend the paper accordingly. \n\nMy 2nd point of concern about non-conjugate likelihood is not correct, and it does seem that the optimal q(u) with non-conjugate likelihood does not follow the same form, so I think this requires a separate investigation. Therefore, I withdraw my concern in this respect.", "rebuttal": "We thank the reviewer for constructive feedback. We attempt to answer the reviewer’s questions below.\n\n*Weakness 1*: even when we don’t use $q(\\mathbf{f} | \\mathbf{u})$ at test time, the test performance can be improved in two ways: more accurate $q(\\mathbf{u})$ and better inducing point $\\mathbf{z}$, and/or better model hyperparameters. Bui et al (2025) did some ablation studies to isolate the impact, showing that (i) the improved T-SVGP bound improves both $q(\\mathbf{u})$ and the hyperparameters, and (ii) using $p(\\mathbf{f} | \\mathbf{u})$ at test time does not degrade the performance whilst significantly reducing prediction time. We expect the same observation for BT-SVGP and can repeat their studies in the next version. In Figure 2, the key difference between the methods is the hyperparameters, noted in lines 192-194.\n\n*Weakness 2*: We disagree on your comments on applicability. We would like to clarify that the collapsed bound is available (equation 8), and a version that supports data minibatching is also available (equation 10). Thus, using these introduces *no* additional variational parameters compared to SGPR/SVGP and T-SGPR/T-SVGP. The implementation and computational overhead over SVGP and T-SVGP is minimal, as explained in lines 108-112.\n\n*Weakness 3*: Thanks for pointing this out. We will provide code and scripts to reproduce the experiments as promised.\n\n*Question 1*: The performance improvement comes from better $q(\\mathbf{u})$, inducing inputs, and hyperparameters. If we use SGPR with the same hyperparameters, the performance of T/BT-SGPR will be different since the $q(\\mathbf{u})$ and $\\mathbf{z}$ obtained by these methods are different. Please also see our response to Weakness 1.\n\n*Question 2*: in lines 87-88, we assume a general matrix $\\mathbf{C}$ and subsequently write down a variational bound that depends on $q(\\mathbf{u})$, $\\mathbf{z}$ and the hyperparameters. We can optimise this bound wrt $q(\\mathbf{u})$ to give the $q(\\mathbf{u})$ in line 88, resulting in the collapsed bound on the next line. We hope this is clear.\n\n*Question 3*: Thanks for the suggestion. We will include these arguments.\n\n*Question 4*: As noted in lines 268-270, we used random partitioning for all experiments. \n\nThanks again, and please let us know if you have further comments."}, {"review_id": "MxdHm18Lkm", "reviewer": "Reviewer_XNNT", "summary": "The authors propose a larger variational family for inference in inducing-point-based Gaussian process regression. They show that the optimal posterior in this family leads to a tighter variational lower bound (when using variational inference with a reverse KL-divergence), as compared to several recent works (Titsias, 2025; Bui et al 2025). They then show that this larger variational family can be used within the power expectation propagation framework (Bui et al 2017) to lead to new inference schemes for scalable Gaussian process regression. The authors show the performance, primarily in terms of root mean square error and log likelihood of held-out data.", "questions": "- Equations 5 and 6: Why does $\\bf{M}$ not appear? Was an optimal form already plugged in given the values of all other parameters? Similarly, in equation 9, have the optimal $\\{\\mathbf{m}_b\\}$ been plugged in?\n\n- In Line 104: Is the square bracket indicating something in the big-O notation? Otherwise simplifying would improve exposition.\n\n- Why would I use equation 5 instead of selecting different values for each block? Would you ever recommend using it? If not, it should be moved to a supplement as essentially an interesting observation.\n\n- `A disadvantage of diagonal and block diagonal structures in M is the expensive predictive covariance.` How expensive would using this predictive covariance be? This should be stated. Could the diagonal elements be calculated efficiently? It looks to me like your experiments only use predictive variances anyway.\n\n- What are the dots above the curves in figure 2? Are these data or the $Z$? This should be stated in the caption.\n\n- What are the solid versus dashed lines in figure 2? This should be stated in the caption or legend.\n\n- What is the (empirical) run time of each method? It seems like (on paper) there is a small amount of additional overhead of the proposed method. But experiments should provide an indication of either run time for each method, or preferably, performance of the methods over time, to give a more detailed sense of the trade-off between performance and run-time for the existing methods versus for the proposed method. This is particularly important because the gains from the proposed method seem to be reasonably small relative to the gains from just increasing the number of inducing points (figure 3).\n\n- Why use the uncollapsed bounds versus the collapsed bounds in 5.2? This requires selection of many more training parameters. See discussion in Section 4.2 of Ober et al 2024.\n\nOber, Artemev, Wagenländer,  Grobins, van der Wilk. \"Recommendations for Baselines and Benchmarking Approximate Gaussian Processes\". 2024.", "limitations": "Yes", "rating": 4, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 3, "originality": 3, "strengths_and_weaknesses": "## Strengths:\n\n**Connections to existing literature.**\nThe authors describe how they build on existing work in a clear way. In particular, the authors repeatedly show how their contribution generalizes prior work in this area.\n\n**Some improved performance in regression.**\nThe authors show some improvements over prior work on several regression datasets.\n\n## Weaknesses:\nOverall, my largest criticism of the paper is that it seems largely written for those who are already both experts in, and interested in, inducing-point-based inference in Gaussian process regression. I'm not convinced that the paper is more widely accessible because of notation not being fully introduced and a large number of equations that are not central to the contribution of the paper. I'm also not convinced the authors have fully motivated the contribution for those outside of this community.\n\n**Motivation.**\nThe motivation for the work is very much presented for someone who already cares a lot about scalable Gaussian process regression. For example the authors propose as a motivating question: ` \"Such improvement begs the question: can we achieve even better approximations by considering more expressive structures for M while preserving efficient computation?\" ` The paper would be stronger if the authors made and supported a claim about improving a prediction task or data analysis, that is more directly of interest to those who are not already familiar with closely related literature.\n\n**Clarity.**\nAt the end of section 4, I still don't have a concrete idea of what I am supposed to practically use as a bound. The other have introduced many. Some of these are tractable, some are recovering existing bounds. It should be more clear which bound or bounds the authors think are candidates for obtaining the best posterior/hyperparameter inference.  I would strongly suggest moving some of the bounds (especially intractable bounds) to the supplement, and focusing on bounds you would advocate for the use of in the main text.\nFigures: Figure 2 has far too many lines that are not visually distinguishable. I can't get a clear takeaway from looking at it and reading the caption. \n\n\n**Introducing Notation.**\nIn many places notation is not clearly introduced.\n\n- Line 27. Notation of u, f, q and p is all introduced at once without sufficient explanation. This is only readable to someone who is already familiar with related work using similar notation.\n- Line 51. The authors do not introduce the notation for the covariance function $k_\\gamma$.\n- Line 57. The authors should state that $\\mathbf{u}$ is a random variable when introducing it. Currently it is a bit hard to unpack how $\\mathbf{z}$ and $\\mathbf{u}$ relate unless someone is already familiar with this notation and related work using it.\n- Equation 1. The notation $f_{\\neq \\mathbf{f}, \\mathbf{u}}$ should be introduced. I think it's technically an abuse of notation (two random variables could be equal without being indexed by the same point in input space, and the stochastic process $f$ can't be equal to its evaluation on a finite collection of points since they live on different spaces.) I don't think the notation is unacceptable, but do think it needs an explanation.\n- General. The authors introduce lots of lower bounds. It is hard to keep track of them. There is a pattern to the notation using subscript $r, c$. The authors should state what this pattern is explicitly. I assume $r$ is meant to stand for regression (although I wonder if $G$ for Gaussian is more appropriate) and $c$ for collapsed, meaning the optimal posterior is calculated for a fixed set of inducing points and \"plugged in\".\n- Line 114: Differentiate between the optimal $\\mathbf{m}$ and a generic $\\mathbf{m}$, e.g. use $\\mathbf{m}^{\\star}$ for the optimal $\\mathbf{m}$.\n- Equation below line 151: Consider noting the dependence on $\\alpha$ in notation on the left hand side.\n\n**Minor:**\nLine 77. Remove \"Note that\". It isn't clear how the author would observe the optimal value of $m_n$ directly from equation 5 and 6.\nLine 147. ` The full derivation is rather lengthy and will be included in the appendix`. This should be a cross-reference to a specific section.\n\n**Question 12 in check list is not answered adequately.**\nJust because a dataset is publicly available does not mean it isn't licensed. The same applies to code. **The authors should cite datasets and software that are integral to this work.**", "was_revised": true, "final_justification": "My primary concern in my initial review was that the paper was primarily targeted to an audience very familiar with the sparse Gaussian process literature, and while containing interesting results the paper did not clarify how these could be practically applied to real analyses. In light of the authors proposed addition of a short section clarifying these practical takeaways, as well as proposed changes to improve clarity and ensure notation is defined, I lean towards accepting the paper.", "rebuttal": "We thank the reviewer for the detailed feedback and will now address their concerns and questions\n\n*Motivation*: We appreciate the suggestion to improve the motivation. We would like to emphasise that efficient methods for Gaussian processes (GPs) have a notable history of strong presence at NeurIPS and are generally of great interest to machine learners and statisticians, as GPs are widely used in many domains – machine learning, robotics, spatial statistics, and astrophysics, to name a few. \n\n*Clarity*: We kept the key bounds in the main text as we thought it would be useful to describe how the bounds are obtained and related to each other. To improve the clarity, we will summarise them in a subsection in the next version.\n\n*Notation*: We attempt to follow standard definitions and notation in the GP and sparse GP literature (Rasmussen and Williams, 2006; Mathews et al., 2016; Bui et al., 2017), and we are happy to clarify them here as well as in the next version.\n\n* $q$, $p$, $f$, $u$: $q$ and $p$ denote the approximate posterior and prior, respectively. $f$ denotes the entire function, i.e., collection of all function values, $\\mathbf{f}$ is the training function values, as explained in line 59, and $\\mathbf{u}$ is the inducing outputs, defined in line 57. The difference between $f$ and $\\mathbf{f}$ is explained in line 62.\n* $k_\\gamma$ is explained in lines 52-53.\n* $\\mathbf{u}$: thanks for the suggestion. We will also mention that $\\mathbf{u}_m = f(\\mathbf{z}_m)$.\n* equation 1: $f_{\\neq \\mathbf{f}, \\mathbf{u}}$ should be read as the function values that are not $\\mathbf{f}$ and $\\mathbf{u}$. \n* lines 114 and 151: we will fix these.\n\n*Questions*\n\n* Question 1 - equations 5, 6 and 9: Yes, the optimal $\\mathbf{M}$ was plugged in.\n* Question 2 - square bracket in line 104: $[N/B]$ simply means the block size.\n* Question 3 - equation 5: Could you please clarify? Did you mean equation 11?\n* Question 4 - expensive predictive covariance: the exact predictive variances are equally expensive, unfortunately, due to the need to invert an NxN matrix, $\\mathbf{D}_{\\mathbf{ff}}$. Bui et al (2025) showed that dropping the expensive term does not hurt the predictive performance.\n* Question 5 - dots in figure 2: these indicate the inducing inputs $\\mathbf{z}$.\n* Question 6 - dotted and solid lines: the solid lines are the predictive mean, and the dashed lines show the 95% credible intervals.\n* Question 7 - runtime: in lines 109-112, page 4, we noted the precise difference between the methods in terms of compute costs. $M$ is typically small, so the overhead is small compared to the overall runtime. For example, on a 2021 Mac laptop, for the kin40k dataset, when M=256, SVGP takes ~15ms for one mini-batch of 500 data points and the overhead when using BT-SVGP is ~1ms, \n* Question 8 - uncollapsed vs collapsed: because of the size of the datasets, the uncollapsed bound is used to enable stochastic optimisation with data minibatching. Using the collapsed bound would be too expensive. Specifically, the sizes of datasets in section 5.2 are 40000 [kin40k], 45730 [protein], 48827 [keggdirected], and 63608 [keggundirected], and thus are out of reach of collapsed bounds.\n\nThanks again, and please let us know if you have further feedback."}]

Author rebuttal:

Meta-review: The authors propose a larger variational family for inducing-point-based Gaussian process regression by replacing a diagonal matrix with a block diagonal matrix.  They show that optimization within the family admits tighter variational lower bounds compared to recent works, and describe how to use the power expectation propagation framework for along with this family for a new inference scheme for scalable GPs.  The authors illustrate the performance of the method in terms of RMSE and log likelihood of held-out data.

During the discussion, tmwE had asked for some performance numbers, which were provided to their satisfaction.  Reviewers kA7g and sm4y expressed some skepticism about the contribution of replacing a diagonal with a block diagonal and the connection to the PITC method; kA7g seemed satisfied with the reviewer response, while sm4y remained of the opinion that going from diagonal to block diagonal seemed a bit incremental.  Otherwise, XNNT commented that this paper seems rather inaccessible to those not already steeped in the lore of variational GPs and maybe familiar already with the prior work on PEP; this sentiment was echoed by reviewer TnWE.  The authors responded that they would add a short clarifying section with context, which seemed to satisfy both reviewers.

The main case for accepting this paper is that the change does seem to get good empirical results with a natural generalization of previous work (replacing a diagonal matrix with a block diagonal); I do not necessarily see this as incremental on its own.  The main case against is the presentation, which seems geared primarily toward specialists.
Decision: Accept (poster)