Adaptive Cholesky Gaussian Processes
Authors: Simon Bartels
Kristoffer Stensbo-Smidt
Pablo Moreno-Muñoz
Wouter Boomsma
Jes Frellsen
Søren Hauberg
Venue: NeurIPS 2022 Submitted (2022)

Abstract: We present a method to fit exact Gaussian process models to large datasets by considering only a subset of the data. Our approach is novel in that the size of the subset is selected on the fly during exact inference with little computational overhead. From an empirical observation that the log-marginal likelihood often exhibits a linear trend once a sufficient subset of a dataset has been observed, we conclude that many large datasets contain redundant information that only slightly affects the posterior. Based on this, we provide probabilistic bounds on the full model evidence that can identify such subsets. Remarkably, these bounds are largely composed of terms that appear in intermediate steps of the standard Cholesky decomposition, allowing us to modify the algorithm to adaptively stop the decomposition once enough data have been observed. Empirically, we show that our method can be directly plugged into well-known inference schemes to fit exact Gaussian process models to large datasets. 

A Adaptive Cholesky Gaussian Processes
Anonymous Author(s) Affiliation Address email
Abstract
We present a method to fit exact Gaussian process models to large datasets by considering only subset of the data. Our approach is novel in that the size of the subset is selected on the fly during exact inference with little computational overhead. From empirical observation that the mor-inagi likelihood often exhibits a linear trend once a sufficient subset a dataset has been observed, we conclude that many large datasets contain redundant information that only slightly affects the posterior. Based on this, we provide probabilistic bounds on the full model evidence that can identify such subsets. Remarkably, these bounds are largely composed of terms that appear in intermediate steps the standard Cholesky decomposition, allowing to modify the algorithm to adaptively stop the decomposition once enough data have been observed. Empirically, we show that our method be directly plugged into well-known inference schemes to fit exact Gaussian process models to large datasets.
3
5
7
10
11
12
13
Introduction
14
15 has been observed (Chalupka al., 2013) that the random- subset-of-data approximation can
16 be a hard-to-beat baseline for approximate Gaussian process inference. However, the question of
17 how to choose the size of the subset is non-trivial to answer. Here we make an attempt. The key
18 computational challenge in Gaussian process regression is to evaluate the log-magginal likelihood of
19 the N observed data points, which known to have cubic complexity (Rasmussen & Williams, 2006).
20 In order to arrive at computationally less expensive approximation of this log- marginal likelihood,
21 we first empirically study its behavior as we increase the number of observations. Figure show this progression for variety of models. We elaborate on this figure in Section 3.1, but for now note that after a certain number of observations, determined by the model and the dataset, the log-marginal likelihood starts to progress with linear trend. This suggest that we may leverage this near-linearity 25 to estimate the log marginal likelihood the full dataset after having seen only a subset of the data. However, the point-of- linearity differs between models, this must be estimated on-the-fly to keep computations tractable.
23
24
27
28 In this paper, approach the problem from (frequentist) probabilistic numerics perspective
29 (Hennig al., 2015). By treating the dataset as collection of independent and identically distributed random variables, common assumption in the frequentist literature, we provide expected upper and
30
31 lower bounds on the log-marginal likelihood, which become tight when the above-mentioned linear trend arises. We provide a particularly efficient algorithm for computing the bounds that leverage
33 the intermediate computations performed by the Cholesky decomposition that is commonly used for evaluating the log- marginal likelihood. The bounds therefore practically free to evaluate. We 35 further show that these bounds allow us to predict when the linear trend determines the full-data log-nn marginal likelihood, such that can phrase optimal stopping problem to determine suitable 37 subset of the data for a particular model. We refer to our method as Adaptive Cholesky Gaussian Process (ACGP). Our approach has complexity of ((M3), where M is the processed subset-size,
Submitted to 36th Conference on Neural Information Processing Systems (NeurIPS 2022). Do not distribute.

---

x 10
Marginal logliiklldood
x 100
Zoom-in
0.0
Cumnlative logply)
0.2
00.4
00.6
-4
0.8
-6
1.0
2-01,a7 0.5
=
-12
Z- 100,0 -20
8-
5m 10.0
200
400
600
800
1000
20 80 100
Observations
Figure 1: The figure shows total log-marginal likelihood function of the size of the training set of different permutations of a simple, synthetic dataset. The different colors correspond to different Gaussian process models, using the squared exponential kernel with length scale l and amplitude a². It can be seen that the log-marginal likelihood exhibits a linear trend after sufficiently many observations have been processed.
39 inducing an overhead of ((M) to the Cholesky. The main difference previous work is that our algorithm does not necessarily look at the whole dataset, which makes it particularly useful in settings where even ( operations are intractable. When a dataset contains large amount of redundant data, ACGP allows the inference procedure stop early, saving precious compute-especiie when the kernel function expensive to evaluate.
41
43
2
Background
44
We PYTHON- PYTHON-inspired inspired index notation, abbreviating for example [yı,...,, 46 that the indexing starts at 1. With Diag we define the operator that sets all off-diagonal entries of matrix to O.
you observe
48
2.1
Gaussian Process Regression
19
We start by briefly reviewing Gaussian process (GP) models and how they are trained (see Rasmussen
& Williams (2006, Chapter 2)). We consider the training dataset D = {Ix, Vn hizz with inputs
50
51 In E IRP and outputs Yn E IR. The inputs are collected in the matrix X [ f NN E NNx, A GP f (P((m(ar) kannaa))((a,e..
is collection of random variables defined in terms of a mean function, (x), and a covariance function or kernel, k(x,
53
cov(s(2) [ f(x')) such 54 that any finite amount of random variables has a Gaussian distribution. Hence, the prior over
55 f == f(X) is [(f: m(X). K#), where we have used the shorthand notation K = k( X, (XX). We
56
will consider the observations y as being noise-corrupted versions of the function values f. and we
57 shall parameterize this corruption through the likelihood function p(3) f), which for regression tasks
58 typically assumed to Gaussian, p(y if) Nff. I). For such a model, the posterior over test ~m,.. . S.), where
59
inputs X. can be computed in closed-form p(f.
m. = k(x X)K-
and
S. = kv X x.) k(X X)K X)
80 with K Kr III By marginalizing over the function values of the likelihood distribution,
61 we obtain the marginal likelihood, p(y) p(y f) the de facto metric for comparing
82 the performance of models the Bayesian framework. While this integral not tractable in
83 general, it does have a closed-form solution for Gaussian process regression. Given the GP prior,
64 p(f) (0, K and the Gaussian likelihood, the oogmaagginal likelihood distribution can be
85 found to be
1 log p(y) log det [K]
N K log 2n.
2
(1)

---

66 2.2 Background on the Cholesky decomposition
Inverting covariance matrices such K is slow and numerically unstable procedure. Therefore, in practice, one typically leverages the Cholesky decomposition of the covariance matrices to compute the inverses. The Cholesky decomposition of symmetric and positive definite matrix K is the unique, lower' triangular matrix L such that K LL (Golub & Van Loan, 2013, Theorem 4.2.7). The advantage of having such a decomposition is that inversion with triangular matrices amounts to Gaussian climination. There are different options to compute L. The Cholesky of 1 1 matrix is the square root of the scalar. For larger matrices,
68
71
73
[chol [ 10,5a
chol
T
chol
Ka.,:: TT
74
where T := ,,, , s.:ee] and s is any integer between and the size of K, Hence, extending given Cholesky to larger matrix requires three steps:
76
1. solve the linear equation system T,
2. apply the downdate K s:,s: TT and
78
3. compute the Cholesky of the down-dated matrix.
An important observation is that Kat,s: - TT is the posterior covariance matrix S. + o²1 when considering X as test points. We will make use of this observation in Section 3.5. The log- determinant of K can be obtained from the Cholesky using log det [K] "n log Lnn- A -1 similar recursive relationship exists between the quadratic form y K y and L y (see appendix, Equation (22)).
79
80
81
82
84
2.3
Related work
Much work has gone into tractable approximations to the log- oog-magginal likelihood. Arguably, the most popular approximation methods for GPs are inducing point methods (Quiñonneo--aan Caadella & Rasmussen, 2005; Snelson & Ghahramani, 2006; Titsias, 2009; Hensman et al., 2013, 2017; Shi et al., 2020; Artemev et al., 2021), where the dataset is approximated through set of pscudo-data points (inducing points), summarizing information from nearby data. Other approaches involve building approximations to K (Fine & Scheinberg, 2001; Harbrecht et al., 2012; Wilson & Nickisch, 2015; Rudi et al., 2017; Wang al., 2019) or aggregating of distributed local approximations (Gal et al., 2014; Deisenroth & Ng. 2015). One may also consider separately the approximation of the quadratic form via linear solvers such conjugate gradients (Hestenes & Stiefel, 1952; Cutajar et al., 2016) and the approximation of the log-determinant (Fitzsimons et al., 201 7a,b; Dong et al., 2017). Another line of research is scaling the hardware (Nguyen al., 2019).
86
86
90
82
93
g4
95
96
All above referenced approaches have computational complexity at least O(N) (with the exception of Hensman et al. (2013) since it uses mini- mini-batching). However, the size of dataset is seldom particularly chosen value but rather the ad-hoc end of the sampling procedure. The dependence on the dataset size implies that more data requires more computational budget even though more data might not be helpful. This is the main motivation for our work: to derive an approximation algorithm where computational complexity does not depend redundant data.
98
100
102
The work closest in spirit to the present paper is by Artemev et al. (2021), who also propose lower and upper bounds quadratic form and log-determinant. There are number of differences, however. Their bound relies on the method of conjugate gradients where we work directly with the Cholesky decomposition. Furthermore, while their bounds deterministic, ours probabilistic, which can make them tighter in certain cases, as they do not need to hold for all worst-case scenarios. This is also the main difference to the work of Hensman et al. (2013). Their bounds allow for mini-batching but these are inherently deterministic when applied with full batch size.
104
106
108
"qquivalently, define L to be upper triangular such that K L L.
3

---

109 3 Methodology
110 In the following, we will sketch our method. Our main goal is convey the idea and intuition. To this end, we use suggestive notation. We refer the reader to the appendix for a more thorough and 112 formal treatment.
118
3.1
Intuition on the linear extrapolation
114
The marginal likelihood is typically presented as a joint distribution, but, one can also view it from cumulative perspective as the sum of cog-conditionas::
115
log p(y) =
Σ log p(yn y:n)
(3)
116
With this equation in hand, the phenomena Figure/ becomes much clearer. The figure shows the value of Equation (3) for increasing number of observations When the plot exhibits linear trend, it is because the summands log p(yn ymm) become approximately constant, implying that the model is not gaining additional knowledge. In other words, new outputs are conditionally independent given the output observations seen so far.
118
120
The key problem addressed in this paper is how to estimate the full marginal likelihood, p(y). from only subset of M observations. The cumulative view of the log-marginal likelihood in Equation (3) is our starting point. In particular, we will provide bounds, which are functions of seen observations,
122
124 on the estimate of the full marginal likelihood. These bounds will allow us decide, on the fly, when
125 we have seen enough observations to accurately estimate the full marginal likelihood.
3.2 Stopping strategy
127 Suppose that we have processed M data points with N M data points yet to be seen. We can then decompose Equation (3) into sum of terms which have already been computed and a remaining sum
log ply) log p(yn v..)
p(yn V.n)
M
A: processed
8: remaining
129
Recall that we consider the Fz. V/s as independent and identically distributed random variables. Hence, we could estimate B as (N - A/M. (1 WM++)....... NI yN) interact non-linearly through the kernel function. Instead, we will derive unbiased lower and upper bounds, L and U. To obtain unbiased estimates, we use the last-m processed points, such that conditioned on the points up to == M m, log p(y) can be bounded from above and below:
130
Yet this is estimator is biased, since
131
132
133
E|C X .. yi.l 5 E X. y.al < EM X.,, V.ol.
and the observations from s to M can be used to estimate L and U. We can then detect when the upper and lower bounds are sufficiently near each other, and stop computations early when the 137 approximation is sufficiently good. More precisely, given a desired relative error r, we stop when
135
U-C
1min(144114)) <r. and sign((/) sign
138 If the bounds hold, then the estimator (4 + 21)2 achieves the desired relative error (Lemma in appendix). This is in contrast to other approximations, where one specifies computational budget, rather than a desired accuracy.
140
141
3.3 Bounds
From Equation (1), we see that the log-marginal likelihood decomposes into the log-determinant of the kernel matrix, quadratic term, and a constant term. the following we present upper and lower bounds for both the log-determinant (Upp and LD. respectively) and the quadratic term (Ug and Lq).
143

---

145 We will need the posterior equations for the observations, that is p(yn V.n) and we will need them as functions of test inputs x, and x²_ To this end, define
148
m (1) (z.) k(z,, X K mn
and
(") (x,, a'.) == k(a,,, x²,) + Sz
k(x,, x
K in,
k(X in x!,)
such that plyn vnn) N(un (1) (Ix). {11) (I kn)), which allows to rewrite Equation (3)
148
as
(Un m. (n (n-1) (In)))
log log (n (1x ann +
(5)
5"' --- 1) (xn 2.n)
149
This reveals that the log-determinant can be written as sum of posterior variances and the quadratic form has an expression as normalized square errors. Other key ingredients for our bounds are estimates for average posterior variance and average covariance. Therefore define the shorthands
150
151
V : Diag ( (X s:M, X x
and
C=
E (6)
(0,+20 (2
2,+2i- 1))
1=1
152
where e; E I is the j-th standard basis vector. The matrix V is simply the diagonal of the posterior covariance matrix E. The matrix C' consists of every second entry of the first off-diagonal of Σ.
153
154
These elements are placed on the diagonal with every second element being 0. The reason for taking every second element is of theoretical nature, see Remark the appendix. In practice we use the full off-diagonal.
155
158
157
3.3.1
Bounds on the log-determinant
158
Both bounds, lower and upper, use that log det [K] = log del [Kis,ss] + log del which follows from the matrix-determinant lemma. The first term is available from the already processed datapoints. is the second addend that needs to be estimated, which we approach from the perspective of Equation (5). It is well-established that, for a fixed input, more observations decrease the posterior variance, and this decrease cannot cross the threshold o² (Rasmussen & Williams, 2006, Question 2.9.4), This remains true when taking the expectation
ا (X x.))
160
162
the input. Hence, the average of the posterior variances for inputs s: is with high probability an overestimate of the average posterior variance for inputs with higher index. This motivates our upper bound the log-determinant:
184
186
Up log det [K se.ssl
+ (N s)μD
(6)
1 712
PD
log (V:) average log posterior variance
167
To arrive at the lower bound on the log-determinant, we need an argument about how fast the average posterior variance could decrease which governed by the covariance between inputs. The variable pn measures the average covariance, and we show in Theorem 6 in the appendix that this overestimates the decrease per step with high probability. Since the decrease cannot exceed x², we introduce vi to denote the step which would cross this threshold.
188
170
VD $
LD = log det [ [...s>] + (ØD s) UD
PD) + (N - - logoz.
2
2 m
PD
2,25 I/ average square covariance
m
DD == max(N, 1 (ND - 1l0--ogg"" o²)
!/ steps can decrease by
(8)
172 Both bounds collapse to the exact solution when 8 = N. The bounds are close when the average covariance between inputs, pn, is small. This occurs for example when the average variance is close
174 2 since the variance an upper bound to the covariance. Another case where pp is small when points are not correlated to begin with.
5

---

178
3.3.2
Bounds on the quadratic term
Denote with r. := yı: m. ! (xx) the prediction errors (the residuals), when considering the first points as training set and the remaining inputs test set. Analogous the bounds on the
177
178
log- deg-determinant, one can show with the matrix inversion lemma that y K K-y y V. :.:
179
x (X)) rr.. Again, the first term will turn out to be already computed. With slight abuse [Σ(
180
of notation let r. ys:M m!!) [x xs:M), ) that is, we consider only the first m entries. Our lower bound arises from another well-known lower bound: a 2 2a
(*)
181
b² Ab for all b (see for example Kim & Teh (2018); Artemev et al. (2021)). We choose b == al where a is chosen
182
183
184
to maximize the bound. The result, after some cancellations, the following lower bound on the
quadratic term: Lo Vi K.KK J.. + (N - 8)0: (2nq apq)
1QQ == r rx
|/ average square ermor
m
N PQ == T.
22
186
The a maximizing above bound is "a/p2. which is the value we chose in our implementation.
187
However, note that Lo is an expected lower bound only if a depends on variables with index smaller
188
than s.
189
Our upper bound arises from the element-wise perspective of Equation (5). We assume that the
expected mean square error -y. (n- "mmmm"""(" (a,n""(xx)) (zx,)) decreases with more observations. However, though mean square error and variance decrease, their expected ratio may increase decrease depending on the choice of kernel, dataset and number of processed points. Using the average error calibration with correction for the decreasing variance, we arrive at our upper bound on the quadratic term:
190
192
Uq V.s K U:s + (N - s) (v' + Pl2)
(10)
1/Q = r. V r. m
I average error calibration
N
PQ = m o4 Cr. Cr. // uvverage increase in error calibration
In our implementation use slightly different upper bound. The estimate of the possible decrease of the variance uses the same technique as the lower bound for the determmnmnn. Therefore we can define analogue to Equation (8) determining the step when the variance estimate falls below o², In our implementation, addends of the quadratic after this step are estimated by the more conservative
195
196
197
2 uQ- Again, the bounds collapse to the true quantity when s = N. The bounds will give good 198
199
estimates when the average covariance between inputs, represented by the matrix C, is low or when the model can predict new data well, that is, when r. close to O.
200
3.4
Validity of bounds and stopping condition
202
For the upper bound on the quadratic form, we need to make a (technical) assumption. It expresses the intuition that the (expected) mean square error should not increase with more data model should not become worse as its training set increases.?
204
205
Assumption
1. Assume that
E [((,, x' (') -m. () (x:))) X Vio < E [( x') m (+) (x)))
X Vio
206
for all s E (1,...,N) for all < S N, where xxxxx') is either
(.,)) or Yا'! (m,,))
Theorem 2. Assume that (x1, vI).... (æN UN) are independent and identically distributed, assume 208 that Assumption holds, and assume that a in the definition of La depends only on I1, y1,... .. E. Un-
2 Empirically, we confirmed this assumption for all experiments considered in Section 4.1 and in Appen- dices B.3.1 to B.3.3.

---

209
For any [l.....N] the bounds defined in Equations (6), (7). (9) and (10) hold in expectation:
[ X os Vial S Ellog det [K] xs: M.wl - BuD X: Y.sl and
ElLq] :x, V:sl
K y
Xx,, y:s)
E( X U:s)
Proof and proof sketch can be found in the appendix. 210
211
Theorem 3. Let > 0 be a desired relative error and set U := Up + Uq and L := CD + Lq. If the stopping conditions hold, that is, sign(U) sign( () and Equation (4) is true, then log p(y) can be estimated from (L + / such that, under the condition LD S log(dct [() s Up and LQ S
213
y' "ky S Uq. the relative error is smaller than v, formally:
214
llog p(y) (u + /2 S <rllog(y))
(11)
215
The proof follows from Lemma 16 in the appendix.
218
Theorem 2 is a first step to obtain a probabilistic statement for Equation (11), that is, a statement of
217
form P (1
log p(ar) +u st
log p(y)
>)) S Theoretically, we can obtaín such a statement
using standard concentration inequalities and union bound over In practice, the error guarding constants would render the result trivial. A union bound can be avoided using Hoeffding's 's inequality for martingales (Fan et al., 2012). However, this requires to replace == M by stopping time independent of M, which we regard as future work.
219
221
222
3.5
Practical implementation
The proposed bounds turn out to be surprisingly cheap to compute. If we set the block-size of the Cholesky decomposition be m, the matrix Σ(() is exactly the downdated matrix in step 2 of the algorithm outlined in Section 2.2. Similarly, the expressions for the bounds on the quadratic form appear while solving the linear equation system L y. A slight modification to the Cholesky algorithm is enough to compute these bounds on the fly during the decomposition with little overhead.
225
226
228
The stopping conditions can be checked before or after Step of the Cholesky decomposition (Section 2.2). Here, we explore the former option since Step 3 is the bottleneck due to being less parallelizable than the other steps.
229
230
231
Note that the definition of the bounds does not involve variables æc y which have not been processed. This allows an on-the-fly construction of the kernel matrix, avoiding potentially expensive kernel function evaluations. Furthermore, it is not necessary to allocate ((N2) memory in advance; a user can specify a maximal amount of processed datapoints, hoping that stopping occurs before hitting that limit. We provide the pseudo-code for this modified algorithm, our key algorithmic contribution, in the appendix. Additionally, we provide a PYTHON implementation of our modified Cholesky decomposition and scripts to replicate the experiments of this paper.
233
234
235
236
237
Experiments
239 We now examine the bounds and stopping strategy for ACGP. When running experiments without GPU support, all linear algebra operations are substituted for direct calls to the OPENBLAS library (Wang et al., 2013), for efficient realization of in-place operations. To still benefit from automatic differentiation, we used PYTORCH (Paszke et al., 2019) with a custom backward function for log p(y) which wraps OPENBLAS. The details of our experimental setup can be found in Appendix A.
241
243
244
4.1
Bound quality
245
In this section we examine the bounds presented in Section and compare them to those proposed by Artemev et al. (2021, Lemma 2 and Lemma 3) (CGLB). Specifically, for the determinant we compare to their ((N) upper bound (Artemev et al., 2021, Eq. I) and their log(det [Q]) as lower bound.
246
247
We set the number of inducing inputs M for CGLB to 512, 1024, 2048, and 4096. For ACGP, define m == 40- 256 = 10 240 which the number of cores times the default OPENBLAS block
249
TTe code is available the following repository: anonymized
7

---

10400
512
ACGPP
100000
ACGP
S13r
1000
30733
43360
30229
OCILB
OGURR
-2.5-
150000
107
1004
-1.
7'K%
42300
500000
30233
42303
200000
-2:
107
-00000
ACTPP
OCILD
30
1141
150
2000
1500
150
35
1041
150
2000
200
Erton socortcle
socords
(a) Upper and lower bounds the (b) Upper and lower bounds on the quadratic term.
(c) Upper and lower bounds the full log-marginal likelihood.
log-determinant term.
Figure 2: Comparison of the upper and lower bounds for ACGP and CGLB on the metro dataset using the OU kernel with a length scale of log l and the time it takes to compute them. The black line indicates the result obtained using exact GP regression with points above and below it marking the upper and lower bounds, respectively. The experiment was repeated five times with different seeds to illustrate the variability in the computation time, shown here as multiple points of the same color. For ACGP the number near the points shows M, the size of the used subset; for CGLB it is the number of inducing inputs. The color of the points reflects these numbers to help discern the size of the subset or number of inducing inputs from the repeated experiments.
size for machines. We compare both methods using squared exponential kernel (SE) and the Ornstein-Uhlenbee kernel (OU),
251
(s, 2) == exp (- | x112
x z
kou = exxp
212
252
where we fix o² = 10-3 and 0 =- 1, and we vary l as log l E 1-1, As benchmarking datasets we use the two datasets consisting of more than 20000 instances used by Artemev et al. (2021): kin40k and protein. We further consider two additional datasets from the UCI repository (Dua & Graff, 2019): metro and pm25 (Liang et al., 2015). We chose these datasets in addition as they are of similar size, they are marked as regression tasks and no data points are missing.
253
254
255
256
257
Empirically, CGLB seems to better estimate the quadratic term, whereas ACGP is faster to identify the log-determinant. Figure 2 shows typical example. Note that, for the quadratic form, the upper bounds tend to be less tight than the lower bounds. Generally, there is no clear winner; sometimes ACGP estimates both quantities faster and sometimes CGLB, For other results, see the appendix.
259
260
261
The reason why CGLB has more difficulties to approximate the log-determinant is that the bound involves trace K Q] where Q is low rank approximation to K. If K2 is of high rank, the gap in the trace can be large. For CGLB the time to compute the bounds is dominated by the pivoted Cholesky decomposition to select the inducing inputs. This overhead becomes irrelevant for the following hyper-parameter tuning experiments, since the selection is computed only once in the beginning. One conclusion from these experiments is to keep in mind that when high precision is required, simply computing the exact solution can be a hard-to-beat baseline.
263
265
266
267
288
4.2
Application in
tuning
We repeat the hyper-parameter tuning experiments performed by Artemev et al. (2021) using the same set-up, see Appendix A for details. We use the same kernel function, a Matérn and the same optimizer: L-BFGS-B (Liu & Nocedal, 1989) with ScIPY (Virtanen et al., 2020) default parameters. Artemev et al. (2021) report their best results using A = 2048 inducing inputs. For reference, we also compare against Sparse Variational Gaussian process regression (SGPR) by Titsias (2009) initialized with the same 512, 1024 and 2048 inducing inputs as CGLB. We use root mean square error (RMSE), negative log predictive density (NLPD) and exact, marginal log-likelihood on the training set, log p()), as performance metrics. The results for all experiments discussed in this section can be found in Appendix B. Here, we will focus the behavior of each method during training.
270
272
274
276
277
278
possible application of ACGP is that an optimizer can decide how precise function evaluations need to be. To explore this possibility, we successively decrease the "relative change in function
8

---

Metrie Evolution over GP training protein
Metric Evolution over GP training netro
OGLB
ACOP
300001
Boad
Foad
XCGP
OCLB
25000
-Iepply)
1000
2000)
2000
2000
0-
1P
socorale
socorade
(a) protein dataset. The iteratively increasing preci- sion may allow ACGP to reach better solutions faster than exact inference the price of later convergence.
(b) metro dataset. Function evaluations with CGLB are generally the fastest at the cost of plateauing at. higher objective function values.
Figure 3: Typical examples of the evolution of the exact log marginal likelihood p(y) while optimizing hyper-paramettrs. See Appendix B.1 for additional plots for all datasets, as well for SVGP runs.
value" " (ftol) convergence criterion of L-BFGS-B as (2/3) and set this as value for r. With this choice, ACGP does not have any more free parameters than a standard optimizer. The block size 282 is a problem independent parameter and it is set to the same value as in Section 4.1.
280
281
We explore two different computing environments. For datasets smaller than 20 000 data points, ran our experiments on a single GPU. The results can be summarized in one paragraph: all methods converge the latest after two minutes. The time difference between methods is less than twenty seconds. Exact Gaussian process regression is fastest, more often than not. The results can be found in Appendix B.I. We conclude that in an environment with significantly more processing resources than memory, approximation may just cause overhead.
284
285
286
287
288
For datasets larger than 20 000 datapoints, our setup differs from Artemev et al. (2021) in that we use only CPUs on machines where the kernel matrix still fits fully into memory. On all datasets, ACGP is essentially exhibiting the same optimization behavior as the exact Gaussian process regressor, just stretched out. ACGP can provide results faster than exact optimization but may be slower in convergence Figure 3a shows for the protein dataset. This observation is expected. However, approximation can also hinder fast convergence as Figure 3b reveals on for the metro dataset. CGLB
290
292
283
294
benefits from caching the chosen inducing inputs and reusing the solution from the last solved linear equation system. The algorithm is faster, though it often plateaus at higher objective function values. The results for kin40k are similar to protein and the results for pm25 are similar metro. These results and the evolution of the root mean square error over time can be found in the appendix. Again, when the available memory permits, the exact computation a hard-to-beat baseline. However, the Cholesky as a standard numerical routine has been engineered over decades, whereas for the implementations of CGLB and ACGP there opportunity for improvement.
296
297
298
299
300
301
Conclusions
302
In this paper we have revisited the use of Cholesky decompositions in Gaussian process regression. We have shown that the Cholesky decomposition almost computes expected lower and upper bounds 305 on the marginal log-likelihood associated with GP regression. With only small modifications to this classic matrix decomposition, we can use these bounds to stop the decompasition before all observations have been processed, This has the practical benefit that the kernel matrix K does not have to computed prior to performing the decomposition, but can rather be computed on-the-fly.
303
307
309
Empirical results indicate that the approach carries significant promise, but no clear winner can be crowned from our experiments. In general, we find that exact GP inference leads to better behaved optimization than approximations such as CGLB and inducing point methods, and that a well-optimized Cholesky implementation is surprisingly competitive in terms of performance. An advantage of our approach is that is essentially parameter-free. The user has to specify a requested numerical accuracy and the computational demands will be scaled accordingly. Finally, we note that
310
311
312
313
315
ACGP complementary to much existing work, and should be seen as an addition to the GP toolbox,
rather than a substitute for all existing tools.

---

References
317
318
Artemev, A., Burt, D. R., and van der Wilk, M. Tighter bounds on the log marginal likelihood of gaussian process regression using conjugate gradients. In Meila, M. and Zhang, T. (eds.), Proceed- ings of the 38th International Conference on Machine Learning, volume 139 of Proceedings of Machine Learning Research, pp. 362-372, 2021.
320
322
Chalupka, K., Williams, C. K. I., and Murray, I. framework for evaluating approximation methods for Gaussian process regression. Journal of Machine Learning Research, 14(1):333 2013.
323
Cutajar, K., Osbore,, M., Cunningham, J., and Filippone, M. Preconditioning kernel matrices. In
325
Balcan, M. F, and Weinberger, K. Q. (eds.), Proceedings of The 33rd International Conference Machine Learning, volume 48 of Proceedings of Machine Learning Research, pp. 2016.
327
328
Deisenroth, M. and Ng. J. W. Distributed Gaussian processes. In International Conference on
329
Machine Learning (ICML). pp. 1481-1490, 2015.
330
Dong, K., Eriksson, D., Nickisch, H.,, Bindel, D., and Wilson, A. G. Scalable log determinants for gaussian process kernel learning. In Advances in Neural Information Processing Systems, pp. 6330-6340, 2017.
332
333
Dua, D. and Graff, C. UCI machine learning repository, 2019. URL httaaapp/// uci. edu/ml.
334
335
Fan, X., Grama, I., and Liu, Q. Hoeffding's inequality for supermartingales. Stochastic Processes and their Applications, 22(100)3555-3550 2012.
337
Fine, S. and Scheinberg, K. Efficient svm training using low-rank kernel representations. Journal of Machine Learning Research, 2:243-264, 2001.
338
339
Fitzsimons, J., Cutajar, K., Osborne, M., Roberts, S., and Filippone, M. Bayesian inference of log
340
determinants. Elidan, G., Kersting, K., and Ihler, A. T. (eds.), Thirty-Thire Conference on Uncertainty in Artificial Intelligence, UAI 2017, August 11-15, 2017, Sydney, Australia, 2017a.
Fitzsimons, J., Granziol, D., Cutajar, K., Oshorne, M., Filippone, M., and Roberts, S. Entropic trace estimates for log determinants. In Ceci, M., Hollmén, J., Todorovski, L.,, Vens, C., and Džeroski, S. (eds.), Machine Learning and Knowledge Discovery in Databases, pp. 323-338, 2017b.
343
344
345
Gal, Y., Van Der Wilk, M., and Rasmussen, C. F. Distributed variational inference in sparse gaussian process regression and latent variable models. arXiv preprint arXiv: 140.. 2014.
Golub, G. and Van Loan, C. Matrix computations. Johns Hopkins Univ Pr, 4 edition, 2013.
348
Harbrecht, H., Peters, M., and Schneider, R. On the low-rank approximation by the pivoted Cholesky decomposition. Applied Numerical Mathematies, 62(4):428 440, 2012.
349
Hennig, P., Osbore,, M., and Girolami, M. Probabilistic numerics and uncertainty in computations.
351
Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 421(2179), 2015.
353
Hensman, J., Fusi, N., and Lawrence, N. Gaussian processes for big data. In Uncertainty in
354
Artificial Intelligence (UAI), pp. 282-290, 2013.
Hensman, J., Durrande, N., Solin, A., et al. Variational fourier features for gaussian processes. J. Mach. Learn. Res., 18(1):5537-5588, 2017.
356
357
Hestenes, M. and Stiefel, F. Methods of conjugate gradients for solving linear systems. Journal of Research of the National Bureau of Standards, 4966)44994336, 1952.
358
359
Kim, H. and Teh, Y. W. Sealing up the automatic statistician: Scalable structure discovery using gaussian processes. In Storkey, A. and Perez-Cruz, F. (eds.), Proceedings of the Twenty-First International Conference on Artificial Intelligence and Statistics, volume 84 of Proceedings of Machine Learning Research, pp. 575-584, 2018.
361
10

---

363
Liang, X., Zou, T., Guo, B., Li, S.. Zhang, H., Zhang, S., Huang, H., and Chen, S. X. Assessing beijing's 's pm2.5 pollution: severity, weather impact, apec and winter heating, Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 471(2182):201502577 2015.
384
Liu, D. C. and Nocedal, J. On the limited memory bfgs method for large scale optimization. Mathematical Programming, 45(1):503-528; 1989.
367
368
Nguyen, D.-T., Filippone, M., and Michiardi, P. Exact gaussian process regression with distributed computations. Proceedings of the 34th ACM/SIGAPP Symposium on Applied Computing, pp. 1286-1295, 2019.
369
Paszke, A., Gross, S., Massa, F., Lerer, A., Bradbury, J., Chanan, G., Killeen, T., Lin, Z., Gimelshein, N., Antiga, L., Desmaison, A., Kopf, A., Yang. E., Vetito, Z., Raison, M., Tejani, A., Chilamkurthy, S., Steiner, B., Fang, L., Bai, J., and Chintala, S. Pytorch: An imperative style, high-performance deep learning library. In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc F., Fox, E., and Garnett, R. (eds.), Advances in Neural Information Processing Systems 32, pp. 8024-8035. 2019.
372
373
374
375
378
Oundele J. and Rasmussen, C. A unifying view of sparse approximate Gaussian process regression. of Machine Learning Research, 6:1939-1955 2005.
378
379
Rasmussen, C. and Williams, C. Gaussian Processes for Machine Learning. MIT, 2006.
Rudi, A., Carratino, L., and Rosasco, L. Falkon:
An optimal large scale kernel method. In Guyon, I., Luxburg, U. V., Bengio, S., Wallach, H., Fergus, R., Vishwanathan, S., and Garnett, R. (eds.), Advances in Neural Information Processing Systems, volume 30, 2017.
381
383
Shi, J., Titsias M., and Mnih, A. Sparse orthogonal variational inference for gaussian processes, In International Conference on Artificial Intelligence and Statistics, pp. 1932-1945 PMLR, 2020.
Snelson, E. and Ghahramani, Z. Sparse gaussian processes using pseudo-inputs Advances in neural
386
information processing systems, 18:1251 127,, 2006.
387
Titsias, M. Variational learning of inducing variables in sparse gaussian processes. In Artificial intelligence and statistics, pp. 567-574. PMLR, 2009.
388
389
Virtanen, P., Gommers, R., Oliphant, T. E., Haberland, M., Reddy, T., Cournapeau, D., Burovski, E.,
Peterson, P., Weckesser, W., Bright, J., der Walt, S. J., Brett, M., Wilson, J., Millman, K. J., Mayorov, N., Nelson, A. R. J., Jones, E., Kern, R., Larson, F., Carey, C. J., Polat, i., Feng, Y., Moore, E. W., VanderPlas, J., Laxalde, D., Perktold, J., Cimrman, R., Henriksen, I., Quintero, E. A., Harris, C. R., Archibald, A. M., Ribeiro, A. H., Pedregosa, F., van Mulbregt, P., and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing Python. Nature Methods, 17:261-272, 2020.
381
399
394
385
396
Wang, K., Pleiss, G., Gardner, J., Tyree, S., Weinberger, K. Q., and Wilson, A. G. Exact gaussian processes on million data points. Advances in Neural Information Processing Systems, 32: 14648-14659, 2019.
397
399
Wang, Q., Zhang, X., Zhang, Y., and Yi, Q. AUGEM: Automatically generate high performance Dense Linear Algebra kernels on x86 CPUs. In SC '1: Proceedings of the International Conference on High Performance Computing, Networking, Storage and Analysis. pp. 1-12, 2013.
400
402
Wilson, A. and Nickisch, H. Kernel interpolation for scalable structured gaussian processes (kiss-gp). In International Conference on Machine Learning, pp. 1775-1784, PMLR, 2015.
Checklist
405
1. For all authors...
408
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope? [Yes]
407
(b) Did you describe the limitations of your work? [Yes]
11

---

409
(c) Did you discuss any potential negative societal impacts of your work? [N/A]
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
411
412
2. If you are including theoretical results...
(a) Did you state the full set of assumptions of all theoretical results? [Yes]
414
(b) Did you include complete proofs of all theoretical results? [Yes]
3. If you experiments..
(a) Did you include the code, data, and instructions needed to reproduce the main exper- imental results (either in the supplemental material or as a URL)? [Yes] The 50MB file size limit of OpenReview, may prevent us from providing everything during the reviewing phase. In any case, all data, code and results will be available from a github repository after acceptance.
417
419
421
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes]
422
(c) Did you report error bars (e.g., with respect to the random seed after running experi- ments multiple times)? [Yes]
424
425
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, cloud provider)? [Yes]
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
428
(a) If your work uses existing assets, did you cite the creators? [Yes]
(b) Did you mention the license of the assets? [Yes]
(c) Did you include any assets either in the supplemental material or as URL? [Yes]
431
The URL to our github repository is currently anonymized but will be made available after acceptance.
433
(d) Did you discuss whether and how consent was obtained from people whose data you're using/curating? [No] This questions seems not applicable.
434
(e) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [No] We do not process sensitive information.
436
437
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given to participants and screenshots, if applicable? [N/A]
439
440
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, applicable? [N/A]
441
442
(c) Did you include the estimated hourly wage paid to participants and the total amount spent participant compensation? [N/A]
12


Reviews:
[{"review_id": "l186Q1BbZ7m", "reviewer": "Reviewer_shrn", "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. \n\nThe 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$.", "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.", "rating": 3, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "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. \n\n\n## Conceptual Limitations\n\n**The first $M$ points might not be representative of the input space or even the full dataset**\n\nThe 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. \n\nThis is clearly not the case!\n\n**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.\n\nClearly, 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! \n\nFor 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})$. \n\nYet, 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$. \n\nIs 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)$.\n\nIn 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.\n\n**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.**\n\n\n## Theoretical Limitations\n\nThe paper also makes a few theoretical mistakes. \n\n\n**Data Points Are Not i.i.d.; Theorem 2 is Questionable.**\n\nFirst, 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. \n\n$(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.\n\nIt would appear that the main result, Theorem 2, relies on this assumption.\n\n\n**The Upper Bound $\\mathcal{U}_D$ is questionable**\n\nNote that the log determinant of a positive definite matrix is always smaller or equal to the log determinant of its diagonal. \n\nHint: 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$).\n\nHence, 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)$.\n\n\n**The Lower Bound $\\mathcal{L}_D$ is questionable**\n\nFundamentally, 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. \n\nIf 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!", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "L11pV5T9H1v", "reviewer": "Reviewer_xgeQ", "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.\nThese bounds are shown to hold on expectation, which I understand to be taken over the randomness of the subsampling order.\nImportantly, 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.\nIn 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.\nIn 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.\nAn 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.\n", "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?\n", "limitations": "The authors strive to honestly discuss the limitations of the work. Two things that should probably be considered in more detail:\n1. 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.\n2. Assumption 1 is merely claimed to “appear empirically true”. What is the evidence for that? More discussion or experiments would be useful.\n", "rating": 6, "confidence": 3, "soundness": 2, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n* 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.\n* Also the paper is very well written and presents technically challenging ideas in an accessible manner in a small amount of space.\n \nWeaknesses:\n* 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.\n* 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.\nIn 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.\n* 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.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "-xrhce6ozAp", "reviewer": "Reviewer_rfLj", "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. \nthe 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.\n\nThe paper presented the bound and proof and numerical simulation study shows that a tight bound can be achieved using some subset data, s.\n", "questions": "Please see Strengths and Weakness. ", "limitations": "yes. ", "rating": 4, "confidence": 3, "soundness": 2, "presentation": 2, "contribution": 1, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "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. \n\nThere are some uncertainties and unclearness about the method.\n\na)\tBounds  (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.\nb)\tThe 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.\nc)\tI 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.\nd)\tCould authors comment on the behavior of bounds with s/M for different model dimensions? \ne)\tFrom 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? \n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Reject