New Bounds for Sparse Variational Gaussian Processes

Thank you for the insightful comments.

I would say that the evidence for 2 and 3 is convincing, while the evidence for 1 is somewhat weak. If this is to be a central claim of the submission, it should be supported with more empirical evidence, in particular in the simulated data setting where "bias" has a particularly straightforward meaning.

The authors note the tendency of SVGP to overestimate the observation noise and provide hints that the new objective lessens this bias. It would be great to see a more extensive study of this particular question, especially on simulated data, and reporting more detailed results for other kernel hyperparameters (e.g. the learned kernel scale).

We provide here the learned hyperparameters for the 1-D Snelson dataset:

        | Noise Variance  | Scale/amplitude  | Lengthscale|

Exact GP | 0.0715 | 0.712 | 0.597 |

SVGP-new | 0.087 | 0.485 | 0.615|

SVGP | 0.108 | 0.331 | 0.617 |

which show that SVGP-new has less bias than SVGP in this example. We plan to add further simulated regression examples in the appendix by following Reviewer's suggestion.

I would love to get more details on the nature of the learned vs. Can the authors report a histogram and/or summary statistics

Yes, we can add histograms for the learned $v_i$s for some of the GP regression datasets (e.g. Snelson, Pol, Bike, Elevators) in the appendix. We give here some summary statistics (min, median and max values) of the learned $v_i$s for the 1-D Snelson example (from Figure 1):

|min | median | max |

|0.172 | 0.952 | 0.9998 |

This indicates that there are few $v_i$s having quite small values but most of them are close to one in this example. We will update the Appendix to add histograms that visualize the final values of $v_i$s.

Regarding the scalar values of $v$ for the non-Gaussian Poisson regression runs, we report here the learned values. For the toy Poisson regression the learned value was around $v=0.675$. Interestingly, for the real Poisson example in NYBikes, the value gets very small below $0.01$. Perhaps this suggests that in this example the expected Poisson log-likelihood reconstruction term has a very strong influence in ELBO optimization, so that the optimization prefers to set $v$ very small (which allows to reduce the posterior variance in $q(f_i)$ since the term $v( k_{ii} - q_{ii})$ becomes small). If the number of inducing points becomes sufficiently large and each $k_{ii} - q_{ii}$ becomes small, then of course the learned $v$ will become close to 1.

| Have you considered adding a new $N$-dimensional "mean shift"

This is a great suggestion. Yes, we tried to do this but we were unable to find to way to add a drift vector $\Delta$ that gives a computational efficient ELBO of $O(N M^2)$ cost. Here, is a short derivation of this. In the KL divergence $KL[q(f|u) || p(f|u)]$ and since the $q(f|u)$ and $p(f|u)$ do not have the same mean anymore the term $\frac{1}{2} \Delta^\top (K_{f f} - Q_{f f })^{-1} \Delta$ appears which has $O(N^3)$ cost. If we try to add some preconditioning to $\Delta$, e.g., $(K_{ff } - Q_{ff})^{1/2} \Delta$ then the $O(N^3)$ remains (although this time will appear in the expected log-likelihood term). We can add an appendix about these derivations, since they could be useful for future research.

We would like to thank the reviewer. We respond to the main comments below.

Yes, mainly. The proposed methods and evaluation criteria are appropriate for the problem at hand, and the experiments effectively demonstrate the benefits of the tighter variational bound. However, the datasets used in the evaluation seem somewhat simple, and the paper could have further strengthened its empirical validation by including a more complex benchmark, such as MNIST or a real-world dataset. This would provide additional evidence of the method’s scalability and effectiveness in practical applications.

Thank you for this comment. Here, we have focused mainly on the GP regression, large GP regression and we provide also an experiment with Poisson regression which is a non-Gaussian likelihood example.

The main weakness of the paper is that the results demonstrate only marginal improvements with the new bound across all evaluated tasks, without a clear example where the standard approach would fail without it. This makes it difficult to assess the practical necessity of the proposed refinement. Additionally, in the non-Gaussian likelihood setting, the paper opts for a simpler approximation rather than leveraging the more expressive variational form, seemingly due to computational constraints.

Please note that we do observe noticeable improvements in predictive performance in some experiments (Pol, Bike, Kin40k, Protein, Buzz). For some datasets, like Pol and Kin40k, the improvement is significant. We believe that a very practical feature of the new ELBO is that in GP regression it requires a minor modification to existing code. So a practitioner can easily train a GP with the new ELBO while keeping the computational cost the same as in the previous SVGP bound.

Tighter sparse variational Gaussian processes, Bui et al (2025), Under review TMLR.

Thank you for pointing to concurrent work. We will cite this work in the next version of our paper.

Table 2 is confusing. The top methods are baselines but have nothing to do with your method I assume. The bottom methods both have ours but isn't your method only SVGP-new?

We agree with the reviewer. We will modify the table to keep "ours" only for SVGP-new as suggested by the reviewer. The upper part of table gives two strong baselines from the literature (discussed in the first paragraph in the Related Work) that are based on different-type of extensions of the SVGP that allow to increase the number of inducing points. Indeed, they are unrelated to our method which replaces $p(f|u)$ by $q(f|u)$, but we have included them for comparison reasons.

When moving to the minibatch setting, the paper is essentially minibatching a KL term for the $q(f|u)$ . How does this work exactly and how do you scale all the terms.

Starting from the following expression for the ELBO (Equation 18 in the paper)

we observe that each term in the sum (i.e., the full term inside the big brackets) depends only on a single data point $(x_i, y_i)$, which is what is needed for minibatch training. Based on this, we can obtain an unbiased ELBO (and gradient) using a minibatch of $b$ data as

where note also that $E_{q(u)} [\log \mathcal{N}(y_i | k_{f_i u} K_{u u}^{-1} u, \sigma^2 )]$ is analytic. We can include an appendix to describe the above details.

It’s unclear how the hyperparameters interact with the new approximation. Do they lead to better-calibrated predictive variances? Providing more details or empirical results on this aspect would help clarify the impact of the proposed method.

After training with the new ELBO, we still use the previous standard SVGP predictive density as discussed in Section 3.1. So the only difference is that the new ELBO can provide different hyperparameters (and inducing inputs $Z$). Figure 1 and Table 1, 2 indicate that this can lead to better predictions in terms of test log-likelihoods. We can try to add an ablation to further study the effect on the predictive variances.

| In the non-Gaussian case, the paper resorts to the standard SVGP prediction, as using the proposed approach directly would be computationally too expensive. This suggests that the main benefit of the method in this setting is primarily improved hyperparameter estimation rather than a fundamentally better approximation of the posterior. If so you should also mention some related work in this ara.

Just to clarify here that for both Gaussian and non-Gaussian likelihoods we do predictions with previous SVGP predictive equations. We can clarify this further in the paper.

Thank you for your comments. Below we provide some responses.

The key insight, replacing $p(f|u)$ with a diagonal-covariance $q(f|u)$, is sound. However, the gap between diagonal V and spherical V is not thoroughly assessed.

Please note that in the medium-size regression experiments reported in Table 1 and Figure 2 we do compare with the method that assumes spherical $V = v^* I$ which uses the optimal scalar value $v^* = \left( 1 + \frac{\text{tr}(K_{f f} - Q_{f f})}{N \sigma^2} \right)^{-1}$. This is precisely, the method denoted in the experiments as "SGPR-artemev". We will clarify this further in the paper. Note, that currently we briefly explain the spherical case $V = v I$ for GP regression and the connection with Artemev et al's bound in Related Work (Section 4) and also in Appendix B.4. From Table 1 and Figure 2 we can observe that diagonal V does work better than the scalar $v$.

The insight presented in this paper is promising, but the details may require further discussion. The key question is: to what extent does assuming a spherical V influence the results? The potential risks associated with this assumption are not explored, nor is there a theoretical or experimental analysis addressing its impact.

We agree with the reviewer that the fact that the non-Gaussian likelihood case requires spherical $V = v I$ is not ideal. In seems that for such non-Gaussian likelihoods the only option that works is to use a spherical $V = v I$ and, as explained in Section 3.3, the diagonal $V$ (that works for GP regression) has cubic cost when trying to obtain each marginal $q(f_i)$. Notice also that if someone heuristically tries to use diagonal $V$ for the non-Gaussian likelihood case and approximate the marginal $q(f_i)$ by

(which is not the correct marginal under $q(f|u) q(u)$ since the variance term $v_i (k_{ii} - q_{ii})$ is wrong) then this creates an inconsistency in the variational distribution in the ELBO, since the $q(f_i)$ used to compute the expected log-likelihood term will be inconsistent with the $q(f|u) q(u)$ in the KL divergence term $KL[q(f|u) q(u) || p(f|u) p(u)]$), and the objective is not a rigorous ELBO anymore.

We can add further discussion to clarify that the spherical $V= v I$ is needed for non-Gaussian likelihoods in order to obtain a rigorous ELBO.

Thank you for very accurately describing the contribution of the paper and for pointing to concurrent work. As also mentioned in the response to Reviewer Pyzb below, we plan to discuss the concurrent work in the next version of our paper.