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: This paper proposes speeding up Gaussian Process inference by estimating the log marginal likelihood with a subset of the data inexpensively chosen on the fly. 

The core idea is to bound the full log marginal likelihood by two quantities $L$ and $U$ that may be estimated inexpensively with a small subset of the data grown incrementally, to stop growing the subset when the error $|U-L|$ is small enough, and to use as estimate for the log marginal likelihood $(U+L)/2$.
strengths_and_weaknesses: While the core idea above is interesting, the implementation of this idea by this paper has a few fundamental limitations, both conceptual and theoretical. 


## Conceptual Limitations

**The first $M$ points might not be representative of the input space or even the full dataset**

The fundamental idea underpinning the stopping strategy proposed in this paper is that, if we estimate the log marginal likelihood of $M$ data points selected sequentially, and we realize that the last $s$ points we added were somehow redundant in the estimation of the log marginal likelihood of the $M$ points, then we can conclude that the first $(M-s)$ points selected are good enough to compute the log marginal likelihood of ***any*** $n \gg M$ points. 

This is clearly not the case!

**Counter-Example:** As a simple counter-example, consider the following toy 1D regression problem: $x_i = i\frac{2\pi}{n}, ~ i \in [1, n]$ with $n=1000000$, and $y_i = \sin(x_i)+\epsilon_i$ where $\epsilon_i$ are i.i.d. Gaussian noise terms with standard deviation $1/100$. Basically, a noisy sine function.

Clearly, inputs $x_i$ are sorted and, no matter the interval $[0, \alpha]$ you consider, you do not need nearly as much as $1000000\frac{\alpha}{2\pi}$ points to represent the sine function! You don't even need $100\frac{\alpha}{2\pi}$ points! 

For most $M$ you'll take, you'll find that points $(x_{M-s}, \dots, x_M)$ are redundant (too many points relative to how fast the function is changing) given $(x_1, \dots, x_{M-s})$. 

Yet, no matter $x_M \in [0, 2\pi)$, knowing $(y_1, x_1), \dots, (y_M, x_M)$ is clearly not sufficient to estimate the full log marginal likelihood for the simple reason that you have no clue what the latent function can possibly look like beyond $x_M$. 

Is the latent function smooth, is it piecewise smooth, does it undergo changepoints after $x_M$, etc. All these scenarios would yield drastically different full marginal likelihood and you simply cannot tell them apart simply from observing $(y_1, x_1), \dots, (y_M, x_M)$.

In general, because points $(y_{M-s}, x_{M-s}), \dots, (y_M, x_M)$ are redundant relative to points $(y_1, x_1), \dots, (y_{M-s}, x_{M-s})$ does not mean that $(y_1, x_1), \dots, (y_{M-s}, x_{M-s})$ are representative of the latent function on the **whole input domain**, which is what you need.

**The main takeaway here is that the order in which you consider your inputs $x_i$ matters! At the very least, your dataset needs to be shuffled, you need to revise your algorithm to make sure it is not exposed to the peculiarities of a specific permutation, and you need to provide some analysis on the effect of this random shuffling on overall performance.**


## Theoretical Limitations

The paper also makes a few theoretical mistakes. 


**Data Points Are Not i.i.d.; Theorem 2 is Questionable.**

First, it is said a couple of times that data $(y_i, x_i)$ are assumed **i.i.d.** (e.g. Page 4, Line 125, and Theorem 2). That's incorrect and (obviously) never the case in Gaussian Process modeling. If that was the case, the Gram matrix $K$ would be diagonal and you would not have any issue estimating the marginal likelihood. 

$(y_i, x_i)$ are usually assumed **i.i.d. conditional on the latent function**! In other words, it is the noise terms affecting observations that are assumed to be i.i.d., not the observations.

It would appear that the main result, Theorem 2, relies on this assumption.


**The Upper Bound $\mathcal{U}_D$ is questionable**

Note that the log determinant of a positive definite matrix is always smaller or equal to the log determinant of its diagonal. 

Hint: Relate the log-det of a psd matrix to the entropy of a multivariate Gaussian, note that the entropy of a multivariate distribution is the sum of the entropies of its marginals and the entropy of its copula, and note that the entropy of a copula is always non-positive (the independence copula is maximum-entropy and has entropy $0$).

Hence, in Equation (6), the term $(N-s)\mu_D$ should really be the sum of the posterior variances of $y_{s:N}$. No need for the incorrect i.i.d. assumption on $(y_i, x_i)$.


**The Lower Bound $\mathcal{L}_D$ is questionable**

Fundamentally, is there really any hope to lower-bound the log-det of the conditional covariance matrix? Depending on how correlated any two unobserved outputs are unconditionally, this conditional covariance matrix can have as low a determinant as possible, irrespective of any observations we have made. 

If any two unobserved outputs are perfectly correlated, which depends solely on the kernel, then the log-det of the conditional covariance matrix will be $-\infty$, no matter what outputs $y_{:s}$ were observed!
questions: Please address the concerns in bold above. Additionally, what are the fundamental differences between what you are trying to achieve and inducing point selection?
limitations: Yes.
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: This paper proposed a generic marginal likelihood approximation in Gaussian process models by using a subset of the data. A logged marginal likelihood is written as  a summation of recursive logged predictive densities. Although a predictive density is analytically calculated, the computation associated with a covariance matrix becomes expensive easily with a number of data points. 
the data is divided into two folders, a processed A (size M) and a remaining B (size N-M).  A logged marginal likelihood has two terms, a determinant and a quadrative term. Instead of taking a recursive sum, a logged marginal likelihood is bounded using a subset of the data (say, s-samples in A folder) and the logged marginal likelihood is approximated by the median of a lower and an upper bounds.

The paper presented the bound and proof and numerical simulation study shows that a tight bound can be achieved using some subset data, s.

strengths_and_weaknesses: The marginal likelihood is a model evidence that it is an important estimation, particularly for model comparison. Any reliable evidence approximation will receive lots of attention. 

There are some uncertainties and unclearness about the method.

a)	Bounds  (both  a determinant and a quadratic term) were formed by for by a posterior for s data in the A folder and a weighted/adjusted prediction for the remining data in A (M-s). Such argument makes sense for one prediction, p(x_{s+1}|x_{s}) and I do not get why this is a reasonable formation for a summation of recursive predictions.
b)	The approximation will be sensitive to M (size of the A folder) and s (subset used for a posterior). It is clear that  as s approaches N, we will get the exact evidence. I feel the paper did not investigate enough to present the approximation quality with s and M both theoretically and numerically.
c)	I don’t think authors did not provide any guide on how to choose s and M. In the simulation study, only M values are shown, no s values.
d)	Could authors comment on the behavior of bounds with s/M for different model dimensions? 
e)	From the practical point of view, s is likely to be chosen with a tight bound using some threshold. When we compare models, it is crucial that we get the same comparison result using evidence approximations. Could authors comment on whether your method will be reliable in this context? 

questions: Please see Strengths and Weakness. 
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: The paper provides a stopping criterion for random data subsampling (without replacement) for Gaussian process inference by computing upper and lower bounds to the final log marginal likelihood.
These bounds are shown to hold on expectation, which I understand to be taken over the randomness of the subsampling order.
Importantly, the bounds can be computed as a by-product from computing the Cholesky decomposition of the kernel matrix that would need to be computed anyway by the standard exact inference algorithm. That means that one can determine the random subset without computational overhead (at least in terms of asymptotic order of growth) and on stopping one already ends up with the Cholesky decomposition of the kernel matrix with respect to the data subset.
In particular, the authors use the standard decomposition of the log marginal likelihood into a quadratic and log determinant term, and develop bounds for both that can be computed efficiently from the partial Cholesky decompositions.
In contrast to previous methods for large scale GP inference, e.g., those that construct representative subsets of inputs or pseudo-inputs, the presented methods can be computed without seeing all the input data. This makes it appealing for very large scale datasets where small samples can be sufficient to approximately determine the posterior.
An experimental evaluation of the idea is provided that shows that the two individual bound pairs (log determinant and quadratic) can perform rather differently compared to a related approach that uses conjugate gradients. In particular the bound on the log determinant seems to be more efficient.

strengths_and_weaknesses: Strengths:
* The idea is to my knowledge novel and appealing because of its principal independence of the overall sample size and its mathematical aesthetics as a function of the partial Cholesky decompositions. One could imagine these bounds to be useful in various contexts.
* Also the paper is very well written and presents technically challenging ideas in an accessible manner in a small amount of space.
 
Weaknesses:
* The bounds are only shown to hold on expectation and, while the derivation of a full probabilistic guarantee might be out-of-scope for this work, there seem to also be no experiments that test in a representative range of settings the probability with which the bounds hold in practice.
* Generally, the tested settings seem rather limited with only one dataset considered in the main text and two more in the appendix. Given that the theory is still in early stages a little more empirical testing would be useful.
In particular, since the idea seems so appealing, one would expect to see some synthetic settings that demonstrate that the proposed method can outperform existing approximate and exact inference methods by a wide margin. However, such settings are not developed.
* Also for the datasets that are considered (at least in the main text), per dataset we only see the performance of the bounds of one component of the log marginal likelihood. While the analysis on this level of detail is certainly interesting, the overall performance would depend on the interplay of both bounds per dataset.

questions: What is the reason for the peak in the log marginal likelihood in the second hyperparameter optimisation experiment? I might not understand the experimental setup here. Should these curves not be monotonically decreasing?

limitations: The authors strive to honestly discuss the limitations of the work. Two things that should probably be considered in more detail:
1. While the computation of the stopping criterion does not affect the asymptotic order of growth of the computational complexity there seems to be some notable overhead that leads the method to be computationally outperformed by exact inference in some cases. This overhead should be discussed in more detail.
2. Assumption 1 is merely claimed to “appear empirically true”. What is the evidence for that? More discussion or experiments would be useful.
review_title: null

Meta-review: metareview: This paper proposes some nice ideas on speeding up Gaussian process inference based on approximating the marginal using subsamples. However, several reviewers noted gaps and potentially flaws in the technical details. The reviews as well as detailed replies during the rebuttal period will help the authors prepare a stronger revision, but the work is not airtight and is not ready for publication in its current form
confidence: Certain
Decision: Reject