Venue: NeurIPS 2022 Conference

Reviews:

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 paper presents a Gaussian process method focusing on data in which the number of samples and the number of dimensions are both high. The GP prior is built upon Bezier surfaces and is formulated as a combination of Gaussian random variables and coefficients taken from Bernstein polynomials. The paper suggested adjusted priors and scalable construction via Bezier buttress. Experiments are performed on various data sets. 
strengths_and_weaknesses: Strengths:

- The problem of handling high-dim GP is interesting , the approach is novel to tackle the problem in GP research.
- The paper gives a good analysis for the method and provides a scalable solution for GP.

Weakness:
- I have one concern of how well Bezier d-surfaces can approximate data. What is the error as the number of dimensions increases? If there is any reference for this, it would be nice to add the clarification.
- It might be out of scope in this paper, there is no theoretical justification for that the proposed approach can address effectively high-dim data .This may be related to the previous point and the way the paper suggests constructing Bezier buttress. The empirical results (first plot of Fig 5.) may. The main concern is that we do not know whether the ability to handle high-dim data comes from deep networks or from the proposed approach.

questions: - How to perform the matrix inverse of $k(Z,Z)$? Is it done analytically by using inverse squared Bernstein.  If so, I would suggest to put more explanation or reference for inverse squared Bernstein and Bernstein polynomials which I feel is missing in this paper. (I’m not familiar with the literature here).
 
- The construction of $f(x)$ somewhat resembles deep structures (of course with no activations), the comparison can be done extensively with other models like deep Gaussian process or deep kernel learning. However, the paper reports just one result with DGP for the “year” data set. Are the results for the remaining data sets (buzz, house electric, slice) available?

Minor points:
	-The optimization is somehow pragmatic. Since the paper suggests KL terms may be big at initial steps, I think KL annealing techniques can be used here. 

limitations: Please see above.
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 manuscript proposes a methodology to extend the Gaussian process framework towards datasets that are both high-dimensional (wide) and have a large sample size (tall) using multivariate Bézier curves. In particular, the authors propose to perform variational inference over a multivariate linear model over Bernstein polynomials having exponentially many control points as parameters. To render the model scalable, the authors propose two approximations: They use a factorisation assumption to reduce the number of control point parameters and they neglect the order of the summation over the dimensions. Experiments on several UCI benchmarks conclude the paper.
strengths_and_weaknesses: 1) Clarity
The paper is mostly well written and the technical content is rather well accessible.

2) Originality
The proposed polynomial covariance structure using Bernstein polynomials seems to not have been explored before.

3) Significance
The method is (currently) specific to large scale regression and could be useful for GP regression. However, in the current state, the underlying approximation assumption need to be better analysed before the model can be used as a black box.

4) Reproducibility
Unfortunately, only some vague outline of the code is contained in the appendix suggesting that there might be an implementation that might be shared at some point. The relevant classes i.e. BezierButtress are not contained hence this cannot be judged. Assuming that the code will be shared, the experiments should be well reproducible. If not, this should be rather difficult.

5) Empirical analysis
The manuscript contains experiments on several standard regression datasets regarding marginal likelihood and prediction error.
a large number of numerical results gathered in a big table regarding marginal likelihood, prediction error and runtime. Some more results on more subtle aspects of the covariance function itsel would have been much appreciated.
questions: * Computations with Bernstein polynomials usually require special care regarding numerical stability i.e. De Casteljau's algorithm. Do you encounter issues in that respect?
* The "adjustment" of the Bernstein polynomials might be numerically challenging as shown by Figure 2 right Panel. Is this correct? What is the condition of the linear system in Equation (10)?
* How would experimental results change without that correction?
* I'm missing an analysis and more detailed evaluation on the effects of the two approximations made to render the model scalable: the factorisation structure and the order of summation.
* In particular, how is the marginal likelihood bound affected? Is it still a bound or a mere approximation that cannot be trusted as a comparative measure in the experiments? Columns 2,4,6 of Table 2 suggests exactly that behavior. Then, one could only compare RMSE and not probabilistic model properties. More specifically, the comparisons in Figure 5, in the left panels and Table upper part would be questionnable in their meaning.
* What would need to change if the likelihood was non-Gaussian? I'm assuming that the same construction can be used.
limitations: Yes.
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 work introduces a type of kernel for Gaussian process regression that requires only a linear computational cost in both number of observations and input features. The scaling is achieved by means of the introduction of the Bézier buttress structure, which allows approximate inference without computing matrix inverses or determinants.
strengths_and_weaknesses: Originality: The work seems original and broadly different to the classical sparse Gaussian process approaches.

Quality: The experiments provided in the paper can be useful for practitioners in the field of Gaussian process models and can open an alternative path for researching about the scalability of Gaussian process models. 

Clarity: The work is generally well written and not difficult to follow, though it presents some problems in the mathematical notation.

Significance (importance): Nowadays the application of Gaussian process models in the context of large datasets with a wide number of features in the input space is still limited due to scalability problem. Therefore, alternative GP approaches that aim to enhance the applicability of GP in such contexts of large and wide data are of great value for the research community.
questions: ---Specific comments---

lines 13 to 16:
It might be useful to include references for Gaussian process model and the approximate methods for their scalability.

line 33:
Include references for Bayesian optimisation, active learning and reinforcement learning.

line 63:
As per Titsias' work (2009), I would suggest to use inducing variables when referring to $\boldsymbol{u}$, and inducing inputs/locations or points when referring to $\mathbf{Z}$.

line 65:
bold for $\boldsymbol{\Sigma_u}$?

line 72:
Why was the lower case variable $\mathbf{p}_{i_1,...,i_d}$  from Eq. (4)

changed to upper case variable $\mathbf{P}_{i_1,...,i_d}$?

is it a matrix now? What is the dimensionality of $\mathbf{P}_{i_1,...,i_d}$? Is it the same of Equation (1) where it was defined $\mathbf{p}_i \in \mathbb{R}^D$?
Is it just a scalar? if so, then most of the bold notations are problematic.

line 73:
It is necessary to introduce the mean $\boldsymbol{\vartheta}_{i_1,...,i_d}$ 

and covariance $\boldsymbol{\Sigma}_{i_1,...,i_d}$ of the Gaussian variables. 

Since they are defined in bold it means they are not scalars, doesn't it?

line 81:
Typo? shouldn't it be $\mathbf{0}$ since the mean is a vector?
Also why bold "$i$" in "for all $\mathbf{i}$"? What is $\mathbf{i}$?

line 82:
why is the covariance matrix $\boldsymbol{\Sigma}_{i_1,...,i_d}$ equal to a scalar? 

I mean in $\boldsymbol{\Sigma}_{i_1,...,i_d}=1$.

line 90: What is $1_{\nu_{\gamma}+1}$? The term was not introduced. 

Should it be $(A_{\gamma})_{i,j}$? It is not clear how both equations in (10) are related. 

line 91:
How is $\zeta_{\gamma}(i_{\gamma})$ related to equation (10)? It is not clear.

line 95:
The statement: "This adjustment works up to $\nu_{\gamma}= 25$, after which negative values occur." is an empirical fact? or is there any mathematical explanation for such a fact? Does it happen regardless of the dimensionality?

line 104 (Figure 2):
Typo, get rid of "is" in "the variance is of $f$ is more..."

---The comments below will appear separated a bit strange, but it is due to a problem in the reviewing platform---

line 108 :
It is a bit confusing to write, $\mathbf{P}_{i_1,...,i_d}$ following a Gaussian with

the mean $\boldsymbol{\hat{\vartheta}}_{i_1,...,i_d}$ 

and covariance $\boldsymbol{\hat{\Sigma}}_{i_1,...,i_d}$

since it was already defined for the prior, $\mathbf{P}_{i_1,...,i_d}$ (after Eq. (7)) to follow a Gaussian with

the mean $\boldsymbol{\vartheta}_{i_1,...,i_d}$ 

and covariance $\boldsymbol{\Sigma}_{i_1,...,i_d}$.

Actually for the "variational control points" it should be something like $\mathbf{P}_{i_1,...,i_d}|\mathbf{y}$ following a Gaussian with

the mean $\boldsymbol{\hat{\vartheta}}_{i_1,...,i_d}$ 

and covariance $\boldsymbol{\hat{\Sigma}}_{i_1,...,i_d}$,

given that it is the (approximate) posterior. Nonetheless, I suggest to define in terms of each variational distribution $q(\mathbf{P}_{i_1,...,i_d})$ equal to  a Gaussian with

the mean $\boldsymbol{\hat{\vartheta}}_{i_1,...,i_d}$ 

and covariance $\boldsymbol{\hat{\Sigma}}_{i_1,...,i_d}$,

this in order to ease the understanding of the following equations (15) and (16).  

Specially to differentiate the prior $p(\mathbf{P}_{i_1,...,i_d})$ 

from the approximate posterior $q(\mathbf{P}_{i_1,...,i_d})$.

line 111:
How is $q(\mathbf{P})$ related to each $q(\mathbf{P}_{i_1,...,i_d})$? It need to be defined.

lines 110, 113 and 120:
the notation for the variable $y$ should be bold, I mean a vector or set that stacks all the $y_j$, say $\mathbf{y}$={$y_j$}$^n_{j=0}$. 

line 116:
Why is it mentioned a variable $\sigma^2$ if the likelihood has not even been introduced?

line 120 Eq. (15):
There should be big brackets containing all the terms of the summatory.

line 127:
It might be better to make Figures appear after they are mentioned in the text.  

line 151:
please introduce the variables $1_{\nu_{1}+1}^{\top}$ and $1_{\nu_{d}+1}^{\top}$. 

line 225:
correct wording in "..., to gives insight into when..."

line 243:
Typo? include "be" in "..., $f$ should in the region..."

line 256:
introduce acronym SGPR.

lines 267 to 269:
It might be useful to know what are the results for keggdirected and elevators, but guaranteeing that the test points are not outside of the domain for the splits.

line 274:
The authors explain that for the concrete dataset the variance overfits yielding a non-optimal likelihood, and conjecture this might be due to $n/d$-radio being low. Nonetheless, the $n/d$-radio for energy is not far from concrete's, and we see from the Table 1 that the energy results were much better than SGPR. How did the authors come up with such a conjecture?  

line 325: please upper case names in the references like: Gaussian, Bayesian, Bernstein,

---Other Questions---

I do not remember having seeing it in the paper, how are the predictions of the model performed? How does the computational complexity of the predictions differ from the classical SGPR?

An interesting property of a Sparse GP is that when the number of inducing points approaches to the number of data observations, we would expect to obtain almost the same results of a full (non-sparse) GP. Is this fact tested in the Bézier GP approach, I mean that when the number of control points approaches to the number of observations we would almost end up with the original Non-sparse GP?

Also, if we set the inducing points as exactly the same to the original input observations, we would expect to recover the original full (non-sparse) GP, is it possible to think about the Bézier GP approach as able to recover the full GP where the input data observations are the same control points?

Would it be feasible to think about a Deep Bézier Gaussian Process? What could it be the challenges and limitations?

The construction of the GP allow the use of different orders $\nu_d$ per dimension, when might that case of having different orders per dimension be useful or applicable?
limitations: The authors were critic about the limitations 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: The paper proposes a new Bézier buttress-based kernel that allows for scaling GP regression both in number of inputs (tall data) and in number of input features (wide data). The scalability is shown both theoretically and empirically through evaluation on popular datasets.

strengths_and_weaknesses: *Originality*. The work is original. While there are other papers leveraging Bézier curves for Gaussian processes, the paper provides a way to parametrize Bézier control points to allow for a scalable training and inference (termed Bézier buttress).

*Quality*. The research is sound. The statements the authors make are supported by evidence. The authors make a detailed analysis of empirical results, addressing both advantages and limitations of the suggested method. They provide a review of the related work, and highlight areas of possible further research.

That said, a more detailed overview of performance against order of the Bézier would be nice. 
Another thing beneficial for the paper would be reports on the wall time consumed and computing power used in the experiments. 

*Clarity*. The work is written in a clear language, the presentation is on-point and laid out in a way that facilitates reading. The notation might be cumbersome sometimes.

One thing to pay attention to is in the supplementary material, line 5: the definition of the variance buttress weights. I think there is a forgotten dependence of $v$ on $\gamma$ and also $\zeta_{\gamma_i}$ doesn't make sense, I guess it should be $\zeta_\gamma(i_\gamma)$.
There is also a typo on line 142: "can be parametrized".

*Significance*. Since applying Gaussian processes to high-dimensional inputs is known to be hard, the proposed method that deals with both "wide" (high-dimensional inputs) and "tall" (many inputs) data is a significant contribution to the world of Gaussian process machine learning. 

questions: I have a several questions to the authors.

1. Am I correct in understanding that computing full covariance would be just rendered as another type of "non-linearity" in the Bézier buttress, that is, multiplying be a cross-product of Bernstein polynomials instead of the squares of them?
2. Have you tried to make even larger orders $\nu$ work? Does it make sense to try larger orders of the Bézier?
3. Have you tried to faithfully compute the true posterior (on a small example perhaps) and see how good the proposed variational approximation is?
4. How easy is it to introduce non-Gaussian likelihood?
5. In your opinion, does it make sense to vary the order of each dimension (that is, non-equal $\nu_\gamma$)? It would make sense as not all dimensions are equal. That would just render a Bézier buttress with different layer widths.
6. A major thing is an assumption of the line 142 which makes Bézier buttress possible. This introduces inter-relations between variational parameters. Do you think this is not very restrictive of an assumption? Again, comparing to the true posterior would probably be good enough of an answer.
7. What happens when jointly optimizing likelihood noise and variational parameters? Have you observed a worsening of the performance?

limitations: The authors address the limitations and constraint of their proposed method. The distribution of the data must be supported on an a-priori known compact cube. Another constraint is a tendency to overfit when $n/d$ ratio is small. Orders $\nu$ larger than 25 do not work.
review_title: null

Meta-review: metareview: The submission considers Bernstein polynomial-based kernels for large-scale Gaussian process regression. To deal with high dimensional inputs, two approximations were proposed: a parameterisation to amortise/share parameters between the random control points of the polynomials and ensembling over randomly permuted orders. In addition, a prior adjustment to the polynomial coefficients was introduced to retain the desirable properties of the variance in the prior. 

The reviewers raise some concerns about these approximations in terms of approximation quality and impact on predictive performance, and the overfitting issue. The authors acknowledged some of these concerns in the text and during the rebuttal.  The AC also notes that the "behaviour foundational to practitioners" [Line 250] is not illustrated well, and perhaps a comparison to non-Bernstein polynomial kernels could be provided to demonstrate the differences/benefits.

Overall, all reviewers think the proposed method is interesting and offers a fairly novel approach to scalable GP regression with high dimensional inputs (with some limitations that require further experiments/analyses), and thus the paper should be accepted. We hope the comments/discussions here are useful for the next iteration.
confidence: Certain
Decision: Accept