Sparse Gaussian Processes: Structured Approximations and Power-EP Revisited

We thank the reviewer for the detailed feedback and will now address their concerns and questions

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.

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.

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.

• $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.
• $k_\gamma$ is explained in lines 52-53.
• $\mathbf{u}$: thanks for the suggestion. We will also mention that $\mathbf{u}_m = f(\mathbf{z}_m)$.
• equation 1: $f_{\neq \mathbf{f}, \mathbf{u}}$ should be read as the function values that are not $\mathbf{f}$ and $\mathbf{u}$.
• lines 114 and 151: we will fix these.

Questions

• Question 1 - equations 5, 6 and 9: Yes, the optimal $\mathbf{M}$ was plugged in.
• Question 2 - square bracket in line 104: $[N/B]$ simply means the block size.
• Question 3 - equation 5: Could you please clarify? Did you mean equation 11?
• 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.
• Question 5 - dots in figure 2: these indicate the inducing inputs $\mathbf{z}$.
• Question 6 - dotted and solid lines: the solid lines are the predictive mean, and the dashed lines show the 95% credible intervals.
• 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,
• 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.

Thanks again, and please let us know if you have further feedback.

We thank the reviewer for constructive feedback. We attempt to answer the reviewer’s questions below.

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.

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.

Weakness 3: Thanks for pointing this out. We will provide code and scripts to reproduce the experiments as promised.

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.

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.

Question 3: Thanks for the suggestion. We will include these arguments.

Question 4: As noted in lines 268-270, we used random partitioning for all experiments.

Thanks again, and please let us know if you have further comments.

Thank you for your response.

I think my main question relates to my question 1, or more specifically about the author's claim about $q(\mathbf{u})$ being more preferable when all other hyperparameters are the same: at prediction, this paper and SVGP/SGPR share the same predictive distribution conditioned on $q$. The collapsed bound (eq. 4) suggests that when the variational posterior is $p(f\vert \mathbf{f}, \mathbf{u})p(\mathbf{f}\vert \mathbf{u}) q(\mathbf{u})$ with a Gaussian likelihood, the optimal $q(\mathbf{u})$ is the variational distribution on line 71. If our main objective is to just come up with a variational distribution $q$ and not to optimize eq. 3 with respect to $\theta$, then this goal has already been achieved: the optimal $q$ (line 71) minimizes the KL-divergence. On the other hand, the collapsed bound in e.g., eq. 8 implies the existence of a different $q(\mathbf{u})$, but it is only optimal with respect to a different $q(\mathbf{f}\vert \mathbf{u})$, which is not used in prediction. Assuming these 2 variational approximations of $q(\mathbf{u})$ are different, this would suggest that the author's claim about "a better q" is not accurate, at least when viewed from the perspective of the minimization of KL divergence. Does this make sense?

I understand that the optimization of GP hyperparameters can be improved with a tighter bound, and the experiments do seem to support this claim. My confusion comes from the observation that, after we are done with hyperparameter optimization, could we just use SGPR/SVGP to find the $q(\mathbf{u})$, with all other parameters unchanged? Would this result be significantly different?

We thank the reviewer for engaging in the discussion.

We would like to note that, loosely speaking, there are three types of variables in the bound: (1) $q(\mathbf{u})$ [or $q(\mathbf{u})$'s parameters, i.e., mean and covariance, if not collapsed out], (2) the variational parameters $\mathbf{z}$, the locations of the inducing points, and (3) the model hyperparameters such as kernel hypers and noise variance.

One important point the reviewer might have missed is that the optimal form for $q(\mathbf{u})$ depends on $\mathbf{z}$ (in addition to the model hypers), so even if $q(\mathbf{u})$ takes the same analytical form across methods, a tighter bound will give a better $\mathbf{z}$ that gives a different $q(\mathbf{u})$ numerically, and ultimately makes $q(f)$ closer to the true posterior.

We will include an experiment in the next iteration to make this clear.

Please let us know if the above makes sense and if you have any other questions.

We thank the reviewer for the constructive feedback and the good-faith discussion throughout the review process.

On improvement over random partitioning: Your suggestion is very reasonable. We agree that this correlation-based partitioning heuristic is a relatively small adjustment, so we will implement and evaluate this in the next iteration.

On non-Gaussian likelihoods: As noted in Titsias (2025) [Section 3.3], for non-Gaussian likelihoods, moving beyond the spherical case ($\mathbf{M} = m\mathbf{I}$) is computationally challenging. We agree that finding a tractable way to combine structured posteriors and non-Gaussian likelihoods is a natural and important next step (see line 273 in the conclusion), but believe it merits a dedicated investigation beyond the scope of the current work.

We hope we have addressed your concerns. Please let us know if you still have any further questions.

We thank the reviewer for the questions and feedback. We will now answer the reviewer’s questions.

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).

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.

Concern 3: Thanks for pointing these out - we will clarify them here and also in the text.

• kernel: we used the Matern-3/2 kernel.
• 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.
• 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.
• 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.

Please let us know if you have further questions.

We thank the reviewer for the constructive feedback. We would like to answer your questions and clarify a few misunderstandings.

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.

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.

Weakness 3: Could you please expand and point us to the specific section in the text?

Questions:

• Equations 2, 3, 5, 9: Thanks! We will fix these.
• 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.

Please let us know if you have further questions.

I thank the reviewers for their clarifying responses.

For weakness 3: to make the paper more complete and more accessible to broader audience, a general note more how the techniques used are similar or different compared to Bui et al (2017), and other related settings may be necessary.

This is inline with Reviewer XNNT's comment on readability,

Overall, 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.

We would like to thank all the reviewers for the questions and discussions so far. If you have any remaining questions/comments/requests before tomorrow, please let us know.