Venue: NeurIPS 2022 Conference

Reviews:

final_score: rating: 3: Reject: For instance, a paper with technical flaws, weak evaluation, inadequate reproducibility and incompletely addressed ethical considerations.
review_content: summary: The paper proposes a scalable GP based on Karhunen-Loeve expansions of kernels. The paper takes a fully Bayesian approach where Gibbs sampling is used to infer the kernel parameters, therefore, leading to forward variable selections.

strengths_and_weaknesses: Strengths:

- The variable selection using GP is an interesting problem.

I think the weaknesses outweigh the strengths, include

- The paper seems to be not well-structured, lacking a proper writing of what the problem that the paper tries to address here is.  It would be better to have an introduction section that would clarify those for us as readers.
- The paper does not provide sufficient descriptions or analyses for the method, for example, lacking literature reviews and no comparison with alternative methods. 

questions: How to scale $\mu$ and $\Sigma$ in Eq. (10-11) which involve matrix inversion? If they are evaluated on 500 data points, such estimations are just local. 

limitations: Besides the presentation issues, I think the paper does not provide an extensive comparison with existing models. For example,an obvious baseline is sparse inducing GPs. 

I believe the paper needs an additional revision to be ready for publication.
review_title: null

final_score: rating: 4: Borderline reject: Technically solid paper where reasons to reject, e.g., limited evaluation, outweigh reasons to accept, e.g., good evaluation. Please use sparingly.
review_content: summary: Regarding the scalable Karhunen-Loeve decomposed kernel BSS-ANOVA for Gaussian process modeling, this paper proposed a BIC/AIC assisted forward variable selection method to have fast and accurate predictions for large tabular datasets. The performance of the proposed model was tested on two datasets.
strengths_and_weaknesses: Pros:
1) A forward variable selection method via BIC/AIC is developed in the ANOVA regime for GP modeling;
2) The new method has been applied to time series regression.

Cons:
1) Benefits brought by forward variable selection are unclear;
2) Differences to existing works have not been well discussed;
3) Insufficient numerical experiments to showcase the superiority of proposed model.

questions: 1) Additive Gaussian processes within the ANOVA formulation for high-dimensional inputs have been studied for a long history, see a recent paper
[1] Lu X, Boukouvalas A, Hensman J. Additive Gaussian Processes Revisited. arXiv preprint arXiv:2206.09861, 2022.
These however have not been discussed and compared in this paper. 
2) Besides, the usage of BIC/AIC to perform truncation is not a new idea, what is the difference to previous works? Please highlight it.
3) What benefits can we have when using the ANOVA type decomposition for GP? Reducing computational budget? Figuring out import variables? The authors should make it clear rather than just showcasing the accuracy of time series regression.
4) The proposed model is suggested to be compared to the ANOVA-type GP competitors on more benchmarks to offer solid conclusions.

limitations: yes
review_title: null

final_score: rating: 5: Borderline accept: Technically solid paper where reasons to accept outweigh reasons to reject, e.g., limited evaluation. Please use sparingly.
review_content: summary: The authors describe a method for variable selection in a GP type model for large datasets called BSS-ANOVA. This (existing) approach uses a kernel that has a basis decomposition where the basis functions can be precomputed and the training complexity becomes O(NP) where P is the number of terms from the kernel decomposition used. The papers’ contribution is on selecting the components in this decomposition with non-zero coefficients using the natural order of complexity and a forward selection paradigm. The coefficients and an information criterion is calculated via Gibbs sampling.
strengths_and_weaknesses: The paper describes a novel approach to fit a GP type model for large datasets and I think as such has a place in the literature.

I would like the authors to be more clear about the impacts of the proposed Gibbs sampler, as described below, i.e. both with regards to the novelty of it and whether it changes the fitted model. My assumption is that this method is new, but if it already exists in the literature, I do not think the use of an information criterion for variable selection alone is contribution enough to warrant acceptance.
questions: What is the worst-case complexity of this algorithm? I.e. BSS-ANOVA has complexity O(NP) if I pre-specify the number of terms P, but your algorithm appears to determine that number dynamically. Would it naturally stop at N^2, could it be unbounded or is guaranteed to be better than N^2?

Is the Gibbs sampling algorithm described in Section 1.3 an original contribution of this paper or is it based on previous work? Does the crucial assumption (8) that enables the use of Gibbs sampling reduce the expressiveness of this model in a meaningful way or can it be shown that this assumption is wlog?

Please change the formula references to (N) to avoid confusion and clean up the citations ([Reich et al., 2009] vs Reich et al. [2009]).

Are there any other variable selection methods for BSS-ANOVA in the literature? Is there a way to use e.g. L1 penalties to find solutions faster?
limitations: I can not see any negative societal impact of this work. The authors do address some limitations in that the paper underperforms on some datasets.
review_title: null

final_score: rating: 4: Borderline reject: Technically solid paper where reasons to reject, e.g., limited evaluation, outweigh reasons to accept, e.g., good evaluation. Please use sparingly.
review_content: summary: The authors develop a Bayesian approach to  the Karhunen-Loéve (KL) decomposed kernel BSS-ANOVA
strengths_and_weaknesses: Strengths and Weaknesses:
Key weakness.  Most people in the GP community will not be familiar with the methods used.  So we need a _much_ clearer background, and perhaps more standard notation:
* Why are we using curly $\theta$ rather than $x$ for the input locations?
* Define kth Bernoulli polynomial.
* Tell us the shape of $\mathcal{B}_1$ etc.
* "the full kernel for a system with n inputs is written $\delta \sim MVN(0, \Gamma)$" is a very strange statement.  The kernel is written $\Gamma$, not $\delta \sim MVN(0, \Gamma)$?  Where did $\delta$ come from? Why are we using MVN rather than $\mathcal{N}$            (especially given that it can cause confusion with a matrix-variate Normal distribution.
* $\Gamma = \sigma_0^2 \tau_0^2 + ...$ makes no sense.  $\Gamma$ is a covariance matrix, while $\sigma_0^2 \tau_0^2$ are scalars?
* $\Gamma_{1, i}$ and $\Gamma_{2, ij}$ aren't defined.
* Where are the eigenvalues and eigenfunctions coming from?

Without these clarifications I do not believe most of the target audience will be able to understand the paper.
questions: See strengths and weaknesses
limitations: See strengths and weaknesses
review_title: null

Meta-review: metareview: Reading the reviews, I think there are ultimately two challenges for the authors to address in this work. First I think ends up being a somewhat simple "background for the community" problem, as both several reviewers and the authors in their general comments point out: significantly more background on KL decomposed kernels may be warranted, and this lack of background (along with perhaps some simple notational differences) led to what I felt like were some challenges understanding the full paper.

With that being said, I don't think the above alone should be sufficient by itself to result in rejection, despite this lack of context likely being a primary contributor to final review scores. However, I do agree with reviewers' concerns that there are some concrete comparisons to existing scalable GP literature missing. I think the inclusion of OAK + inducing points is a start, but e.g. the use of m=40 inducing points in section 3.4 is surprising to me, partly perhaps because it's not clear which task this was a problem for (you state early that you use m=200 for the cascaded tanks task), and partly because none of the dataset sizes involved appear to me to be even remotely beyond the capability of these existing scalable GP approximations in the literature -- even up to m=512 or m=1024 inducing points is fairly standard practice. Given the reasonably good performance of OAK in some of the new results even with limited inducing point sets, I think that clearly further investigation is warranted there. Beyond inducing point methods, there are also NNGP / Vecchia models (which have also been recently made variational via Wu et al., 2022), and even exact GPs seem readily applicable to some of the tasks considered here with access to even a single moderately powerful GPU.
confidence: Certain
Decision: Reject