Venue: NeurIPS 2022 Conference

Reviews:

final_score: rating: 6: Weak Accept: Technically solid, moderate-to-high impact paper, with no major concerns with respect to evaluation, resources, reproducibility, ethical considerations.
review_content: summary: The paper contributes a specific inference method for sparse variational Gaussian process regression in a Bayesian setup allowing inference of  hyperparameters which doesn't require sampling the inducing variables u. The new algorithm is derived by exploiting the exact bound derived in Titsias, 2009 to "collapse" the full ELBO (with hyperarameters) leading to a more effective algorithm in the Bayesain setting.

The algorithm is evaluated on a relevant set of synthetic and benchmark regression problems which demonstrates that the proposed method outperforms existing sparse inducing points methods in terms of RMSE, predictive log-likelihood, and wall clock time. It is not clear that the approach consistently outperforms a basic SGPR+MLII approach.
strengths_and_weaknesses: 

- Novelty and signiificance: The approach appears to be leverage a new insights about the ELBO and contributes a sensible addition to the sparse GP regression toolbox. The paper provides a better understanding of GP inference in the sparse regime which is valuable in its own right, although it is arguably not very far reaching and significant from a practical point of view in part because many people simply resort to MLII which also seems to perform on par with the proposed method on most of the UCI benchmarks.
- Quality: The basic idea appears sensible and the core algorithmic work well-executed. I find that the experiments could have been presented and discussed better (see specific suggestions/questions below).
- Clarity and presentation: I find that there are several questions (see below) that I would like to see addressed in order for it to be easily accessible and comphrensive.
questions: 
Comments/suggestions/questions:
- l105 mentions the classification case, although only regression with the Normal likelihood is considered going forward. I'd suggest removing or clarifying.

- Table 2: I am not sure I follow what's going on in the table; e.g. how is the "quality" determined? Why does MAP not support a non-Gaussian likelihood?... I'd suggest adding a reference per entry where these facts are extracted unless clearly explained e.g. in the appendix. Also, perhaps explicitly indicate which approach is yours.

- Algorithm 1 (page 6): I am slightly confused by the comments and notes in the algorithm box (line 2-3 and 11); are the kernel hyperparameters part of the optimization or not? Additionally, I think it would be beneficial to include reference to the relevant equations (specifically for the \mathcak{L} elements )

- Figure 2: Please clarify the caption; I (initially) found it difficult to work out exactly what was happening in the various panels.

- Figure 3: There are several SGPR+HMC variants; I assume the authors prefer the reader to focus on "Adapt 2" but it would be helpful with some guidance.

- Table 4: $\delta$ is not explained/defined

- Sec 5.2: I believe the discussion fails to comment on the performance of SGPR (ML-II) vs. the proposed SGPR+HMC (Adapt 2) as observed in table 4 and figure 4? Would this not suggest that the benefit of a fully Bayesian approach is not evident in most cases?

- l250 (figure on the right): It would be informative with the GPR+HMC and SGPR+HMC performance shown as horizontal lines.

Very minor:
- Conditioning on X seems a bit arbitrary and briefly introduced in 101-108, then disappears in sec 3.1, but back in l130... I'd suggest more consistency/clarity.  

- l101 vs introduction: the text switches from "kernel" to "covariance function" and it is used interchangeably throughout; perhaps consider if a consistent use would be beneficial

- Eq (1) missing comma before f

- l151: "...sample the u..." -> "...sample u...."

- l190: check missing punctuation prior to  "...the outer..."

- l234: I am not sure where this common practice is used; could you provide a reference?

- l248, l266: "don't" -> "do not" I would suggest avoiding contractions.

- l261: not sure what "this" refers to




limitations: See points made elsewhere regarding the discussion of results etc.



------------------------
Following the authors' response, I have updated my score to 6.
review_title: null

final_score: rating: 6: Weak Accept: Technically solid, moderate-to-high impact paper, with no major concerns with respect to evaluation, resources, reproducibility, ethical considerations.
review_content: summary: The paper proposes a novel framework for tackling Gaussian Process (GP) inference by combining sparse GP (based on inducing variables) and Bayesian hyperparameter inference. The proposed method offers great potential to avoid overfitting compared to previous methods that rely on maximizing the marginal likelihood. Experiments include both synthetic data (one-dimensional) and five small/medium-sized UCI benchmarks. Empirical results show that the proposed methods improve existing sparse GP models.
strengths_and_weaknesses: Strengths:
- Considering hyperparameter uncertainty in GP inference is a meaningful and valid direction. The paper thoroughly examined existing methods and proposed a practical method to advance this direction, which brings value to GP's applications.
- All contributions are well motivated and related works are well cited and documented. The entire paper is self-contained and the presentation is logically ideal.
- The proposed method is technically correct, with Algorithm 1 helping understand how the model works.
- Empirical results indicate that the proposed model obtains strong performance on a diverse set of tasks, in particular with reasonable running time.

Limitation: In my opinion, a limitation of the method is that the compared methods are relatively old (SGPR in 2009 [1] and SVGP in 2015 [2]). So it would be interesting to see how this work compares to recent GP inference works (or at least discuss the opportunities of adopting the proposed Bayesian hyperparameter inference to recent works).

Reference
1. Titsias, Michalis. "Variational learning of inducing variables in sparse Gaussian processes." Artificial intelligence and statistics. PMLR, 2009.
2. Hensman, James, Alexander Matthews, and Zoubin Ghahramani. "Scalable variational Gaussian process classification." Artificial Intelligence and Statistics. PMLR, 2015.
questions: N/A
limitations: Yes, the authors have adequately addressed the limitations and potential negative societal impact of their work
review_title: null

final_score: rating: 6: Weak Accept: Technically solid, moderate-to-high impact paper, with no major concerns with respect to evaluation, resources, reproducibility, ethical considerations.
review_content: summary: This paper proposes a new extension of the variational sparse GP framework (Hensman et al., 2015b) to bypass the need of sampling inducing points. This is achieved by a clever analysis of the optimal variational distribution of hyper-parameter $\theta$, which then reveals a way to collapse the inducing points $u$ from the ELBO. This doubly collapsed ELBO is intractable, but its gradient can be estimated via sampling $\theta$ from the optimal variational distribution $q^{\ast}(\theta)$ using HMC. This is a significant contribution for inference with high-dimensional data, since it improves the sampling efficiency. The proposed method achieves promising empirical results on several benchmarks.
strengths_and_weaknesses: The paper is clearly written on a high level, but some details are not quite clear to me. For instance, I do not quite follow what happens at Eq. (14), which was written with extreme brevity. Where does the addition come from? If I understand correctly, it should be an equality there as we substitute $\mathcal{L}^{\ast\ast}$ with $\mathcal{L}^{\ast}$ evaluated at $q=q^{\ast}$. Then, we apply chain rule to get the inner product. Still, I'm quite confused by the final approximation. How do we get to $\frac{1}{J}\sum_{j=1}^{J} \mathcal{L}_{\Phi,\theta_j}$ ? The rest of the math seems to check out, so maybe the authors can help to clarify this part. 
Minor issue: In line 120, it should be plugging back to the Eq. above Eq. (5), which is unlabeled.

As far as I am aware, this work is novel. I really like the clever idea of doubly collapsing the ELBO to avoid sampling inducing inputs. I can see that it would be a problem for inference with high-dimensional data, as we have to draw samples from a very large state space. This was explicitly pointed out by the authors: "While this approach is quite general, in the case of Gaussian regression, it vastly increases the dimensionality of the state space over which HMC must be run relative to HMC in GPR since the u are sampled in addition to the θ". However, I find it quite disappointing that the authors only demonstrate on fairly low-dimensional datasets (e.g., the largest was 13-dimensional). I do not see how these benchmarks justify the use of this method.

----------
I have read the author's rebuttal, which clarified some of my concerns. My overall score will be adjusted to 6. However, I would encourage the authors to continue to revise the empirical study by expanding the range of benchmarks and adding more ablation study to confirm the expected advantage over other methods.

questions: - The warm start strategy was introduced in algorithm 1, but was never described anywhere else. Can the authors elaborate?
- Why do SGPR + HMC and SVGP use different sampler implementations? Could it be that the runtime/RMSE differences are partly due to this choice?
- The latest benchmarks included in the empirical study were SVGP (Hensman et al., 2015a) and JointHMC (Hensman et al., 2015b). Since SGPR+HMC extends these methods, I believe it is appropriate to include them as benchmarks. However, since the introduction of SVGP there have been plenty of other sparse GPs. Table 1 even listed a few, such as Rossi et al., 2021 and Bui et al., 2018. I think this paper should additionally compare to these benchmarks to fully showcase the relevance of the proposed method in light of more recent developments.
- Since I'm not entirely clear about Eq. (14), I'm tentatively putting a score of 5. I will upgrade my score once the authors address my confusion.
limitations: This is a purely theoretical work. I do not foresee any potential negative societal impacts.
review_title: null

final_score: rating: 7: Accept: Technically solid paper, with high impact on at least one sub-area, or moderate-to-high impact on more than one areas, with good-to-excellent evaluation, resources, reproducibility, and no unaddressed ethical considerations.
review_content: summary: The paper addresses the problem of hyperparameter optimization for sparse GP kernels. The proposed solution utilises MCMC to keep the hyperparameter posterior, accounting for uncertainty in the hyperparameter space. The paper introduces a series of decisions, which improve computational efficiency (separating sampling of Zs and \thetas) and help with handling high-dimensional inputs (by using a doubly collapsed ELBO). The experiments show improvement over alternative techniques.
strengths_and_weaknesses: Strengths: The proposed fully Bayesian sparse GP is suitable for larger datasets and has more efficient sampling scheme that the alternatives. The proposed approach is theoretically sound.

The problem is well motivated and presented. The authors nicely demonstrate the problem of kernel hyperparameter optimization with Figure 1 and further compare relevant works and alternative approaches with Tables 1 and 2.

Weaknesses: It would be interesting to see more comparison methods in the experiments section. The benchmark datasets could be better motivated, e.g. by explaining their main differences, adding their descriptions.
questions: Minor issues:
Equation 2 is a little bit weird in a sense that \theta is missing after the first equal sign. I think you should include the dependency on \theta there for p(y|f) and p(f) too.
limitations: Yes.
review_title: null

Meta-review: metareview: The paper presents a new method for scalable inference in Gaussian regression based on the inducing-variable formalism. It shows that it is possible to sample covariance hyper-parameters while avoiding sampling the inducing variables, a consequence of their doubly collapsed bound. The reviewers agree that it is a technically solid paper and the authors have addressed their concerns satisfactorily, with one of the reviewers raising their scores and the authors providing additional results. I believe this work is worth presenting at NeurIPS and, therefore, recommend its acceptance. 
confidence: Certain
Decision: Accept