A Scalable and Exact Gaussian Process Sampler via Kernel Packets
Authors: Haoyuan Chen
Rui Tuo
Venue: Submitted to ICLR 2023 (2023)

Abstract: In view of the widespread use of Gaussian processes (GPs) in machine learning models, generating random sample paths of GPs is crucial for many machine learning applications. Sampling from a GP essentially requires generating high-dimensional Gaussian random vectors, which is computationally challenging if a direct method, such as the one based on Cholesky decomposition, is implemented. We develop a scalable algorithm to sample random realizations of the prior and the posterior of GP models with Matérn correlation functions. Unlike existing scalable sampling algorithms, the proposed approach draws samples from the theoretical distributions exactly. The algorithm exploits a novel structure called the kernel packets (KP), which gives an exact sparse representation of the dense covariance matrices. The proposed method is applicable for one-dimensional GPs, and multi-dimensional GPs under some conditions such as separable kernels with full grid designs. Via a series of experiments and comparisons with other recent works, we demonstrate the efficiency and accuracy of the proposed method.

Under review as a conference paper at ICLR 2023
A SCALABLE AND EXACT GAUSSIAN PROCESS SAM- PLER VIA KERNEL PACKETS
Anonymous authors Paper under double-blind review
ABSTRACT
In view of the widespread of Gaussian processes (GPs) in machine learn- ing models, generating random sample paths of GPs is crucial for many machine learning applications. Sampling from GP essentially requires generating high- dimensional Gaussian random vectors, which is computationally challenging if a direct method, such as the one based on Cholesky decomposition, implemented. We develop a scalable algorithm to sample random realizations of the prior and the posterior of GP models with Matérn correlation functions. Unlike existing scalable sampling algorithms, the proposed approach draws samples from the the- oretical distributions exactly. The algorithm exploits a novel structure called the kernel packets (KP), which gives exact sparse representation of the dense co- variance matrices. The proposed method is applicable for onediiensioonal GPs, and multi-dimensione GPs under some conditions such as separable kernels with full grid designs. Via a series of experiments and comparisons with other recent works, demonstrate the efficiency and accuracy of the proposed method.
INTRODUCTION
Gaussian processes (GPs) have been widely used in statistical and machine learning applications
[Rasmussen] 2003 Cressie 2013 Santner et 2003). The relevant areas and topics include re-
gression (O Hagan 1978 Bishop et al. 1995 Rasmussen | MacKay et al. 2003). classification
(Kuss et al., 2005 Nickisch Rasmussen | Hensman et al.. 2015),. Bayesian networks (Neal, 2012). optimization (Srinivas et al., 2009), and so on. GP modeling proceeds by imposing a GP as the prior of an underlying continuous function, which provides a flexible nonparametrie framework for prediction and inference problems. When the sample size is large, the basic framework for GP regression suffers from the computational challenge of inverting large covariance matrices. A lot of work has been done to address this issue. Recent advances in scalable GP regression include Nyström approximation (Quinonero- Cainonero-Candela Candela & Rasmussen 20055 Titsias 2009] Hensman et al., 2013), random Fourier features (Rahimi & Recht] 2007)) local approximation (Gramacy & Apley 2015), structured kernel interpolation (Wilson & Nickisch) 2015), state-space formulation (Grig
Recht|[000)).
orievskiy al. 2017 Nickisch et al. 2018). Vecchia approximation (Katzfuss & Guinness 202
sparse representation Chen et al. 2022 Ding et al. 2021), etc.
In this article, we focus on the sampling of random GP realizations. Such GPs can be either the prior stochastic processes, or the posterior processes in GP regression. Generating random sample paths of the GP prior or the posterior of the GP regression is crucial in machine learning areas such as Bayesian Optimization (Snoek et al. [012) Frazier 201 8a b), reinforcement learning (Kuss) & Rasmussen| 2003 Engel et al. [005]] Grande et al. 2014), 4), and inverse problems uncertainty [
quantification MMuray-Smith & Pearlmutter] 2004 Marzouk & Najm 2009 Teckentrup 2020).
To generate the function of random GP sample, a common practice is to discretize the input space, and the problem becomes the sampling of high-dimensional multivariate normal vector. Sampling high-dimensional multivariate normal vectors, however, is computationally challenging as well, as need to factorize the large covariance matrices. Despite the vast literature of the scalable GP regression, the sampling methodologies are still underdeveloped.
Existing scalable sampling algorithms for GPs scarce. A recent prominent work is done by Wilson et al.) (2020). They proposed an efficient sampling approach called decoupled sampling by exploiting Matheron's 's rule and combining Nyström approximation and random Fourier feature.

---

Under review as a conference paper at ICLR 2023
They also generalized it to pathwise conditioning (Wilson al. 2021 based on Matheron's 's update, which only needs the sampling from GP priors and a powerful tool for both reasoning about and working with GPs. Motivated by those work, [Madoxx al. (2021) extended Matheron's 's rule for multi-task GPs and applied it to Bayesian optimization; Nguyen et al.
(021)) proposed the first use of such bounds improve Gaussian process posterior sampling. It worth noting that each of the above methods enforces certain approximation scheme to facilitate rapid computation, i.c., none of these methods draws random samples from the theoretical Gaussian distributions exactly.
In this paper, we propose algorithms of sampling from GP priors and GP regression posteriors ex- actly for one-dimension Matérn kernel with half integer smoothness v, then we extend it to noiseless data in dilli-dieensoon. We introduce the kernel packet (KP) (Chen et al. (0022) as a major tool for sampling and reduce the time complexity to O(v)). This produces a linear-time exact sampler if is not too large. Specifically, our work makes the following contributions:
We propose an exact sampling method for Gaussian processes with (product) Matérn cor- relations on one-dimensional or multi-diiensional grid points. The computational time grows linearly in the size of the grid points.
We propose exact and scalable sampler for the posterior Gaussian processes based on one-dimen mee-immnsiionll data or multi-dimennsoonl data on full grid designs.
We demonstrate the value of the proposed algorithm in the Bayesian optimization and dy- namical system problems.
2 BACKGROUND
This section covers the related background of the proposed method. Sections 2.1 and 2.2 introduce GPs, GP regression, and the basic method of GP sampling. In section 2.3 we review newly introduced covariance matrix representation called the Kernel Packet (KP) Chen et al. 2022). which will help expedite the GP sampling.
2.1 GPS AND GP REGRESSION
A GP is a stochastic process whose finite-dimensional distributions are multivariate normal. The probability law of GP is uniquely determined by its mean function nu(-) and covariance function ((s)), and we denote this GP as GP(((-), ((:))))
Let f X + R be unknown function, X' Rd. Suppose the training set consists of Gaussian observations yi (xx; + with noise ~ NN In GP regression, we impose the prior f~ GP((())),) K(.
Suppose that we have observed = (yI, ... on distinct points X = [= The posterior of given the data is also a GP. Specifically, the posterior evaluation at m untried inputs X {x] 1 follows fyly ~ NN n. K.,,, ..) with WWilliams & Rasmussen (2066):
Pajn = pb + K.,n(K.. [nn, + [
(y -nn);
(1)
K. K,,//n = =² Kr.. =2((K... K.,, Knn 0² In] "kn)).
(2)
where a2 is the variance of GP, In is a n n identity matrix, K,n K(x", XX) (K(X, x*) x²² KT (K(x*, (K(x*,x1);"* x1). ((xxxxxn) Kn,n (x)x,,x))]" 2m K.. [K(X] x;)) " Un = (u(xI), ... ,(xx))) and μ.. = (u(xi), ,((,)))))
In this work, we focus on Matérn correlation functions. One-dimensional Matérn correlation func- tions (teinn 1999) are defined
x') K(x,x')=
K(x,
(2 - T(v)
21-V V2
( Kv 2
x' w
w
for any x, x' E R. where x² the variance, v > 0 is the smoothness parameter, w >> 0 is the lengthscale and K, is the modified Bessel function of the second kind. GPs with Matérn correlations
2

---

Under review as a conference paper at ICLR 2023
form a rich family with finite smoothness; their sample paths are [v-1] times differentiable (Satner ct al. 2003). By virtue its flexibility, Matérn family is deemed a popular choice of correlation functions spatial statistics (Diggle et al. (2003) geostatistics (Curriero 2006; Pardo-lguzquiza &
Chica-Olmo 2008). image analysis (Zafari et al. 2020 Okhrin et al. 2020), and other applications.
A common choice of multi-dimensional correlation structure is the "separable" or "product" corre- lations given by
K K(x,x' x') II Ky (z),
(4)
for any X, E Rd where Kj is the one-dimensional Matérn correlation function for each j. Although the product of Matérn correlations doesn't 't have the same smoothness properties with multi-dimensional Matérn correlations, the assumption separability used extensively spatio- temporal statistics (Gneiting et al. 2006; Genton! 2007 Constantinou et al. (017)) because it allows for a simple construction of valid space-time parametric models and facilitates the computational procedures for large datasets in inference and parameter estimation.
2.2 SAMPLING
The goal is sample f(-) GP(((.)), ((+))) To achieve finite representation, we discretize the input space and consider the function values over a set of grid points Z = }z,)e_ and the objective is to generate samples f, (f(Z))..... (...),.... [(zp)) from multivariate normal distribution (Fu(Z), K (Z, = (p K.p)). The standard sampling method of multivariate normal dis- tribution is as follows: 1) generate vector of samples fp.0 whose entries are independent and identically distributed normal, 2) employ the Cholesky decomposition (Golub & Van Loan 2013) to factorize the covariance matrix Ky,p as C₁ Kp.p. 3) generate the output sample r, as
fp f- Cpfp.0 + Mp.
Sampling posterior GP can be done in similar manner. Suppose have observations y (V1. ... yn) on n distinct points X = )ix) where y; = fx,) + c, with noise & ~ NNO0 o?) and f (P((()), K):,))). The goal
generate posterior samples f. ly at m test points X* J.m {x]} Because the posterior samples f.ly ~ (( (..ii, K.,+jo) are also multivariate normal
distributed according to [ and (2) in Section 2.1 can do Cholesky decomposition C.ja CT ./n
Kn,,/n and generate GP posterior sample as
f,ly < Cujn Fn + Malms
(6)
where fn. N...I,,)
2.3 KERNEL PACKETS
In this section, we review the theory and methods of kemel packets by Chen et al. (002)) Suppose K (...) is a one-dimensional Matérn kernel defined in (3). and one-dimensional input points X is Then there exists a collection of linearly inde- pendent functions { K., such that each Φ, = Σグ= A,)((r; Then covariance matrix Kn = K(X, X) connected to a sparse matrix ((X) in the following way:
ordered and distinct. Let I. ((,,n
K(x, I.y) has compact support.
Ax (XX),
where both Ax and (xX) are banded matrices, the (I, t) entry of ((X) is de(x))). The matrix
Ax consists of the coefficients construct the KPs, and specifically, A is the (j, entry of
Ax. In view of the sparse representation and the compact supportedness of dj. the bandwidth of Ax is (k- - 1)/2, the bandwidth of (XX) is (A: 3)/2, k == 2v + 2, is the smoothness parameter in (3) We defer the detailed algorithm to establish the factorization (7) to Appendix A.I and the connections to the state-space GP to Appendix A.3 Here we only emphasize that the algorithm to find A and (xX takes only (v²))) operations and (vn) storage.
3

---

Under review as a conference paper at ICLR 2023
Based on 7 follows:
can substitute Ax and Φ(.) for K(;,)) in I (
and rewrite the equations as
MaIn = pe -l. -x'(x')) + 0²
Ax 1 un)
(8)
+ o²
,| = -
σ²
Because (XX and Ax are both banded matrices, the summations ( + Ax and (X) + of Ax x are also banded matrices. Therefore, the time complexity of matrix inversion via KPs is only O(v")).
SAMPLING WITH KERNEL PACKETS
In this section, we present the main algorithms for sampling from Gaussian process priors and pos- teriors based the KP technique. We introduce the algorithms of sampling from -ie-dimeniional GP and multi-dimensional GP respectively in sections and 3.2
3.1 SAMPLING FROM ONE-DIMENSIONAL GPs
Consider one-dimensional GPs with Matérn correlations for half integer smoothness V.
Sampling from the prior distribution For set of grid points Z = {z, K= we first compute the sparse matrices Az and Φf (ZZ) in (7). Now, instead of direct Cholesky decomposition of the covariance matrix K p.p = K (Z Z), we consider the Cholesky decomposition of the symmetric positive definite matrix Rz == (ZZ) = P. This shifting makes significant dif- ference: Kp,p is a dense matrix, and the Cholesky decomposition takes O(n²) time. In contrast, Rz, a multiplication of two banded matrices, is also banded matrix with bandwidth of no more than 2v. It is well known that the Cholesky factor of a banded matrix is also banded, and the computation can be done O(v"n) time. Suppose Q2 by a vector fp.o. which is standard multivariate normal distributed. It's not hard to show that (AZ) -2Q22p + Mo ~ (Hpp Kp.)). Hence, we obtain the samples which have same distribu- tion with fp. practical calculation, we may first compute the multiplication Qzfp. then compute reduce the time complexity since both Az and Qz are banded matrices. The entire algorithm only costs O()ân) time and (vn) storage. summary given Algorithm [
Rz, then multiply (AZ)-1Q) 'Qz
Algorithm/ Sampling from one- nne-dimensional GP priors.
Input: Ordered input dataset Z defined in section 2.1
1: Compute banded matrices Az and (ZZ) with respect to the input data Z
2: Obtain a banded matrix: Rz == A (Z))
3: Apply Cholesky decomposition Rz, get the lower triangular matrix Qz satisfying Qz Rz:
4: Generate samples f,,0 ~ I))
5: Compute fxp : (AZ) -1Q2fp + Mp.
Output: fxp
Sampling from one-dimensional GP regression posteriors. This can be done by combining the KP technique with the Matheron's 's rule. The Matheron's rule was popularized by Journel & Huijbregts (1976) to geostatistics field. Recently Wilson et al. (2020) rediscovered and exploited it to develop GP posterior sampling approach. The Matheron's 's rule stated as follows. Let and be jointly Gaussian random vectors. Then the random vector a, conditional on b B, is equal in distribution to
(a|b = 3) a + Cov(a, b)Cov(b, b) ( b),
(10)
where Cov(a, b) is the covariance of (a, b).
4

---

Under review as a conference paper at ICLR 2023
By (0). the exact GP posteriors can be sampled by two jointly Gaussian random variables, and obtain
f.ly - f. + K., Kn,
σ² - (y
(11)
where f, and f, are jointly Gaussian random variables from the prior distribution, noise variates (0, I.). Clearly, the joint distribution of (fn, f.) follows the multivariate normal distribution as follows:
)
(fn, f.) ~ .N[
Fnn [/.
K Knn
K
).,,.
(12)
K,,
K
We may apply KP to (11 and get the following corollary:
fy - f. + Φ xx(xx) [(xx+ of
r² Ax ff cc.
f,,
(13)
Given that Algorithm D requires distinct and ordered data points, it's reasonable to assume training set X doesn't t coincide with test set X Also, we can rearrange the combined set of the training X and test set X an ordered set and record the index of reordering. Next, we can utilize the Algorithm [ to draw a reordered vector in Ru+m from the GP prior and obtain jointly samples f. and f, by recovering the reordered vector the original sequence. Finally, we plug f, and f, into formula (13) to calculate the posterior samples. It's obvious the time complexity of this approach is also (v"") due to the sparsity of the matrices to Cholesky decomposed and inverted.
3.2 SAMPLING FROM MULTI-DIMENSIONAL GPs
For multi matti-dimensional GPs, suppose that the points X are given by full grid, defined X , where each X is a set of one-dimensional data points, A B denotes the Cartesian product of A and B. Based on the separable structure defined in (4) and a full grid XFG can transform the multi-dimensional problem to mme-dimensional problem since the correlation matrix K FFG can be represented by Kronecker products (Henderson et al. 1983. Saatçi
2012
Wilson & Nickisch, 2015) of matrices over each input dimension j:
K(XFG XFG)
K2(X (6)))
(14)
Prior To sample from GP priors over a full grid design ZG = x zú), it's easy to compute that tensor product /²= (ATT 2 Qz(n is the Cholesky factor of the correlation matrix with full
grid FG For each dimension j. we can generate matrices Azo o(Z²)), and Cholesky factor QZa with respect to one-dimensional dataset Z(), To get the GP prior samplings over full grids,
need apply the tensor product j=1 A Z Qzoy) to g.0 More specifically, the samples from GP priors can be computed via
FFG
A (AZU Qzw)
PFG + 1pp
N K(Z (20))
(15)
Accordingly, on full grids, we only need to perform for loop over dimensions, multiply banded matrices in each dimension, and calculate tensor product. The total method is also ( II ny) because we are also dealing with sparse matrices, here n, is the number of data points X (1)
Posterior With regard to GP posteriors, it's impossible employ the prior sampling scheme to draw jointly distributed samples (f!T fFG) since we cannot order sequence of multi-dimensional data. We consider directly Cholesky factorizing the posterior correlation matrix KF to obtain
*,*/n
CFG then use (6) to generate samples. The calculation of posterior mean 1F5 and posterior vari-
ance KFG * only costs O I ny) time by using equations 25 and 26 in Appendix A.2
Therefore the time complexity of the entire scheme requires O
nz))
5

---

Under review as a conference paper at ICLR 2023
Log of Time
2-Wasserstein
3.0
RFF Cholesky
2.5
1.0
KP
RFF
0.5
Cholesky
-4
0.0
10a 5- 10t 5-10 10a 5-10²
10*
caso
C3S0
Figure I: Time and accuracy of different algorithms for sampling from one-dimensional GP priors with Matérn 3/2. We denote approach (KP) by blue dots, random Fourier features (RFF) by orange stars, direct Cholesky decomposition by green crosses. Left: Logarithm time taken to generate draw from GP priors, z-axis the number of grid points p. Right: 2-Wasserstein distances between priors and true distributions over ten points for three different cases { [2]12200 {2:) 500
491 when 500.
EXPERIMENTS
In this section, we will demonstrate the computational efficiency and accuracy of the proposed sampling algorithm. We first generate samples from GP priors and posteriors with one-dimensional space section 4.T and two-dimensional full grid designs in section 4.2 Then we consider the same applications as in (Wilson et al.) 2020) and perform approach to real problems in section 4.3 For prior samplings, we use the random Fourier features (RFF) with 1024 features (Rahimi & Recht, 2007) and the Cholesky decomposition method benchmarks. For posterior samplings, decou- pled (Wilson et al.] 2020) algorithm with exact Matheron's 's update and the Cholesky decomposition method are used as benchmarks. We consider Matérn correlations in (3) with lengthscale w = V2v, smoothness 3/2 and 5/2. 2. In -wo-dimensional problems, choose "separable" Matérn corre- correlations and parameters
lations mentioned in Section 2.1 with the dimension. We set the variance as o² = for all experiments. We set seed value as 99 and perform 1000 replications for each experiment. All plots regarding 3/2 are given in the main article and these associated with 5/2 are deferred to Appendix A.5
/2v in each
4.1 ONE-DIMENSIONAL EXAMPLES
Prior Sampling We generate one-dimensional prior samples on uniformly distributed points Z {zHPp over interval [0, 10] with p 10, 50, 100, 500, 1000, 5000, 10000. Left plots in Figure and and Figure show the time taken in sampling schemes for different algorithms over the different number of points p.
observe that KP algorithm costs much less than other algorithms for both Matérn 3/2 and Matérn 5/2 correlations especially when = 5000, 10000. Also, the curves of Cholesky decomposition are incomplete due to the limit of storage. To test the accuracy, select three subsets of size ten: ]- {z,} 259 i=200- { when 500 and compute the 2-Wasserstein distances between empirical priors and true distributions over these three subsets (called Cases 1-3 thereafter). The 2-Wasserstein distance measures the similarity of two distributions. Let fi ~ N(p1, Σn and f2 N[12, Σ2), the 2- Wasserstein distance between the Gaussian distributions f1, f2 on I,2 norm is given by (Dowson & Landau 1982, Gelbrich 1990, Mallasto & Feragen] 2017)
500
W2(f1, (iff..fo f2) = (114 22ll + tr (E) + +2 2( ))))
(16)
For the empirical distributions, the parameters are estimated from the replicated samples. From right plots in Figure and Figure 7 we can observe that KP algorithm outperforms RFF and Cholesky method greatly in Case 1 and Case 3, which are boundary points of It may be because of the numerical precision that Cholesky method has the lower 2-Wasserstein distances in Case 3.
Posterior Sampling We investigate the performance of sampling over m = 1000 test points on the one-dimensional test function, the shifted and dilated Gramacy & Lee function (Gramacy &

---

Under review as a conference paper at ICLR 2023
Log of Time
Log of MSE
-2-
KP
Decoupled Cholesky
Decoupled Cholesky
-10
-12
-14
10a 50 10t 5- 10a 5-
10*
10t 5
10a
10
5-
10a
5 10" 10
Figure 2: Time and accuracy of different algorithms for sampling from one-dimensional GP poste- riors with Matérn 3/2. We denote decoupled method by red triangles. Left: Logarithm of time taken to generate a draw from GP posteriors over m = 1000 points, r-axis is the number of observations n. Right: Logarithm of MSE over m 1000 points.
Log of Time
2 2-Wasserstein
KP
RFF Cholesky
RFF
Cholesky
0.2
0.2
2
2
210
2
0350
cate
Ca50
Figure 3: Time and accuracy of different algorithms for sampling from GP priors full grids with Matérn 3/2. Left: Logarithm of time taken to generate a draw from GP priors. Right: 2- Wasserstein distances between priors and true distributions over nine points for three different cases 15. X {2 131 when {2 = ={ {-5 + 10 10.2-- 2 , 5 1-10---- 10-2-53.
x (2)
Lee 2012 f(x) sin + 5T) 1) (0.2z 0.5)4 over the interval [0, 10) with
10, 50, 100, 500, 1000, 5000, 10000 respectively given a, = 10-3 33 Figure 2 and Figure S show the time cost of posterior GP sampling and Mean Squared Error (MSE) between the posterior samplings and the true function values with the above settings. It's clear that KP consumes less time to achieve better accuracy compared with decoupled approach and direct Cholesky decomposition method. The performance curves of Cholesky factorization in Figure 2 and Figure g in Appendix A.5 are incomplete, because the Cholesky decomposition function fails work for such a huge sample size and gives a runtime error.
4.2 TWO-DIMENSIONAL EXAMPLES
Prior Sampling We generate prior samples over level-7) full grid design: 2,5G
(-=1(-5 +
10 2- -5 + 2 10 2- 2- 10-2-") with ")
6 and 2. Likewise,
Wassersst distances defined in (16) to evaluate the accuracy of prior samplings in three cases, when ]31 131 { {-5 + 10 2-5
10 2 2-5). Left plots in Figure 3 and Figure 9 in Appendix A.5 show that KP takes less time than RFF and the Cholesky method especially for 2 grid points. Right plots in Figures 3 and Figure 9 show that KP has higher accuracy than RFF for both Matérn 3/2 and 5/2 correlations but has lower accuracy than the Cholesky method for Matérn 5/2 correlations.
x x==1==)_11 (z)11=1) x 2z
127
Posterior Sampling
We choose the Griewank function Griewank 1981
d
f(x) = Σ
z
E cos
+1,
x 5)
4000
1=1

---

Under review as a conference paper at ICLR 2023
Log Time
Log MSE
-2.5
Decoupled Cholesky
-5.0
Decoupled
Cholesky
-7.5
-12.5
25
210
2
24
2
2
2
28
2
2
214
2
2
Figure 4: Time and accuracy of different algorithms for sampling from GP posteriors over full grids with Matérn 3/2. Left: Logarithm of time taken generate draw from GP posteriors over 1024 points. Right: Logarithm of MSE over 1024 points.
as our test function and level-r full grid design: X = (-511-----00
10 "] with 17 = 3, 4, ... and = 2 as our design of experiment. We then investigate the av- erage computational time and MSE over Figure 4 and Figure [0 in Appendix A.5 illustrate the performance of different sampling strate- gies,
1024 random test points for each sampling method.
observe that both direct Cholesky decomposition and decoupled algorithm only generate posterior samples from at most 212 observations due to the limit of storage, however, KP algorithm can sample from 2 observations because the space complexity of KP-based computation only requires (vn). Although the accuracy of the KP method is lower than the direct Cholesky decomposition method and decoupled algorithm, is still of high accuracy and the time cost is in much shorter period of time compared with other algorithms.
4.3 APPLICATIONS
Thompson Sampling Thompson Sampling (TS) (Thompson, 1933 is a classical strategy for decision-making by selecting actions E X' that minimize a black-box function f X - R. In round t, TS selects Tz+1 arg (fl minimizer, we may obtain Ut+11 by evaluating at 21+1, and then add (I2+1, Mt+1) to the training set.
y is the observation set. Upon finding the
In this experiment, we consider one-dimensional cut of the Ackley function Ackkey 2012)
f(x) 20 expp -02 v0.5 x²} -expf expf-0.5(cos(222)) + 1)} + exp(1) + 20.
(17)
The goal of this experiment to find the global minimizer of function in (7). We start with & 2v + 2 samples before the optimization, then at each round of TS. we draw posterior sample fly on 1000 uniformly distributed points over the interval [-5, 5] given the observation set. Next, we pick the smallest posterior sample at this round and add it to the training set, and repeat the above process. After some steps, we are able to get closer to the global minimum. In Figure we compare the logarithm of total regret of different sampling algorithms, both the proposed approach (KP) and the decoupled method can find the global minimum within 15 rounds, which outperform the direct Cholesky factorization sampling scheme.
Simulating Dynamical Systems Gaussian process posteriors are also commonly used in dynam- ical systems when we don't have sufficient data. Consider a one-dimensional ordinary differential equation
x' = f(s,z(s)) xs)),
(18)
can discretize the system's equation to a difference equation in the following formula by Euler method (Butcher 2016).
Vt T(+1 Fz (8t,
(19)
where is the fixed step size. We aim simulate the state trajectories of this dynamical system. First, we set an initial point Fa. then at iteration t, we draw a posterior GP sampling from the conditional distribution p(y, D, where D. denotes the set of the data {(x,, H= and the
current trajectory ((r), My) )=) In our implement, choose the model [8, 0.5z
8

---

Under review as a conference paper at ICLR 2023
6.5
Decoupled
Cholesky
Log et reggee
Random Search
Decoupled
5.0
5.00
Cholesky
4.5
Random Search
4.0
40 NUT evaruations
40
60
Num evoluations
Figure 5: Logarithm of regret of Thompson sampling methods when optimizing Ackley function with Matérn 3/2 (left) and Matérn 5/2 (right).
0.05z with step size 0.25 and initial point To -4.5. The training set was generated by evaluating (19) for n = 1000 points {a:)== uniformly from the interval [-5, 5]. Variations in each step were simulated by independent Matérn 3/2 GPs with w = v3. 6 state trajectories of each algorithm for T 480 steps, time cost each iteration, and logarithm of MSE between Ie obtained from GP-based simulations and obtained by directly performing the Euler method in (19) at each iteration The left plot in Figure 6 shows that the KP algorithm can accurately characterize state trajectories of this dynamical system. The middle and right plots in Figure indicate that the KP algorithm takes much less time and yields high accuracy each step. It achieves nearly the same performance as the decoupled method and outperforms the Cholesky decomposition method.
the
Walttime(seconds)
Log of MSE
3300
KP Decoupled
KP
3325
Decoupled
Cholesky
Cholesky
KP
-.00
Decoupled
-4.25
Chalesky
0.50
Ground Truth
0.25
-00
100
200
300
400
500
100
200
200
400
500
100
200
300
400
500
Figure 6: Simulations of an ordinary differential equation. Left: Trajectories generated via different algorithms. Middle: Time cost at each iteration. Right: Logarithm of MSE between simulations and ground truth state trajectories each iteration.
DISCUSSION
In this work, we propose scalable and exact algorithm for one- mne-iimensiona Gaussian process sam- pling with Matérm correlations for malf-integen smoothness v, which only requires O(z²n) time and (un) space. The proposed method can be extended to some multi-dimensional problems such as noiseless full grid designs by using tensor product techniques. If the design not grid- the proposed algorithm is not applicable, we may use the method in (Ding et al. 2020 to devise approximation algorithms for sampling.
While the proposed method is theoretically exact and scalable algorithm, we observe some nu- merical stability issues (see Appendix A.4) in our experiments. This explains why sometimes the proposed method is not accurate the Cholesky decomposition method. Improvements and extensions the proposed algorithm to overcome the stability issues and accommodate general multi-dimensioon observations will be considered in future work.

---

Under review as a conference paper at ICLR 2023
REFERENCES
David Ackley. A connectionist machine for genetic hillclimbing, volume 28, Springer Science & Business Media, 2012.
Christopher M Bishop et al. Neural networks for pattern recognition. Oxford university press, 1995.
John Charles Butcher. Numerical methods for ordinary differential equations. John Wiley & Sons, 2016.
Haoyuan Chen, Liang Ding, and Rui Tuo. Kernel packet: An exact and scalable algorithm for gaussian process regression with matérn correlations. Journal of Machine Learning Research, 23 (127):1-32, 2022.
Panayiotis Constantinou, Piotr Kokoszka, and Matthew Reimherr. Testing separability of space-time functional processes. Biometrika, 104(2):425-477, 2017.
Noel Cressie.
Statistics for spatial data. John Wiley & Sons, 2015.
Frank C Curriero. On the use of non-euclidean distance measures in geostatistics. Mathematical Geology, 38(8):907-926 2006.
Peter Diggle, Paulo Ribeiro, and Ole F Christensen..
An introduction to model-based geostatistics. In Spatial statistics and computational methods, 43-86. Springer, 2003.
Liang Ding, Rui Tuo, and Shahin Shahrampour. Generalization guarantees for sparse kernel approx- imation with entropic optimal features. In International Conference on Machine Learning, pp- 2545-2555. PMLR, 2020.
Liang Ding, Rui Tuo, and Shahin Shahrampour. sparse expansion for deep gaussian processes. arXiv preprint arXiv:21/2.0588 2021.
DC Dowson and BV660117 Landau. The fréchet distance between multivariate normal distributions. Journal of multivariate analysis, 12(3):450-455, 1982.
Yaakov Engel, Shie Mannor, and Ron Meir. Reinforcement learning with gaussian processes. In Proceedings of the 22nd international conference on Machine learning, pp. 201-208, 2005.
Peter Frazier. Bayesian optimization. In Recent advances in optimization and modeling of contem- porury problems, pp. 255-278. Informs, 2018a
Peter Frazier. A tutorial on bayesian optimization. arXiv preprint arXivv 2018b. 2018b.
Matthias Gelbrich. On formula for the 12 wasserstein metric between measures on euclidean and hilbert spaces. Mathematische Nachrichten, 147(1):185-200 1990.
Mare G Genton. Separable approximations of space-time covariance matrices. Emvirommetrics: The official journal of the International Environmetrics Sociery, 18(7):681-66) 2007.
Tilmann Gneiting, Mare G Genton, and Peter Guttorp. Geostatistical space-time models, stationar- ity, separability, and full symmetry. Monographs On Statistics and Applied Probability, 107:151, 2006.
Gene H Golub and Charles F Van Loan. Matrix computations. JHU press, 2013.
Robert B Gramacy and Daniel W Apley. Local gaussian process approximation for large computer experiments. Journal of Computational and Graphical Statistics, 24(2):561-578, 2015.
Robert B Gramacy and Herbert KH Lee. Cases for the nugget in modeling computer experiments. Statistics and Computing, 22(3):713-722 2012.
Robert Grande, Thomas Walsh, and Jonathan How. Sample efficient reinforcement learning with gaussian processes. In International Conference on Machine Learning, pp. 1332-1348 1340. PMLR, 2014.
10

---

Under review as a conference paper at ICLR 2023
Andreas O Griewank. Generalized descent for global optimization. Journal of optimization theory and applications, 34(1):11-39) 1981.
Alexander Grigorievskiy, Neil Lawrence, and Simo Särkkä. Parallelizable sparse inverse formula- tion gaussian processes (spingp). 2017 IEEE 27th International Workshop on Machine Learn- ing for Signal Processing (MLSP), pp. 1-6. IEEE, 2017.
Harold Henderson, Friedrich Pukelsheim, and Shayle R Searle. the history the kronecker product. Linear and Multilinear Algebru, 14(2):113-120, 1983.
James Hensman, Nicolo Fusi, and Neil D Lawrence. Gaussian processes for big data. arXiv preprint arXiv:1309.6835,, 2013.
James Hensman, Alexander Matthews, and Zoubin Ghahramani. Scalable variational gaussian pro- cess classification. In Artificial Intelligence and Statistics, pp. 351-360. PMLR, 2015.
Andre G Journel and Charles Huijbregts. Mining geostatistics. 1976.
Matthias Katzfuss and Joseph Guinness. A general framework for vecchia approximations of gaus-
sian processes. Statistical Science, 36(1):124 1:::::2244::444 2021.
Malte Kuss and Carl Rasmussen. Gaussian processes in reinforcement learning. Advances in neural information processing systems, 16, 2003.
Malte Kuss, Carl Edward Rasmussen, and Ralf Herbrich. Assessing approximate inference for
binary gaussian process classification. Journal of machine learning research, 6(10), 2005.
David JC MacKay, David JC Mac Kay, al. Information theory, inference and learning algorithms. Cambridge university press, 2003.
Wesley J Maddox, Maximilian Balandat, Andrew G Wilson, and Eytan Bakshy. Bayesian optimiza- tion with high- high-dimensional outputs. Advances in Neural Information Processing Systems, 34: 19274-19287, 2021.
Anton Mallasto and Aasa Feragen. Learning from uncertain curves: The 2-wasserstein metric for gaussian processes. Advances in Neural Information Processing Systems, 30, 2017.
Youssef M Marzouk and Habib N Najm. Dimensionality reduction and polynomial chaos accel- eration of bayesian inference in inverse problems. Journal of Computational Physics, 228(6): 1862-1902, 2009.
Roderick Murray Smith and Barak A Pearlmutter. Transformations of gaussian process priors. In International Workshop on Deterministic and Statistical Methods in Machine Learning, pp. 110- 123. Springer, 2004.
Radford M. Neal. Bayesian learning for neural networks, volume 118. Springer Science & Business Media, 2012.
Vu Nguyen, Mare Peter Deisenroth, and Michael A Osborne. Gaussian process sampling and opti- mization with approximate upper and lower bounds. arXiv preprint arXiv:2/ 10. /2087, 2021.
Hannes Nickisch and Carl Edward Rasmussen. Approximations for binary gaussian process classi- fication. Journal of Machine Learning Research, (0cc): 2008.
Hannes Nickisch, Arno Solin, and Alexander Grigorevskiy. State space gaussian processes with gaussian likelihood, In International Conference on Machine Learning, pp. 3789-3798. PMLR, 2018.
Anthony O' O'Hagan. Curve fitting and optimal design for prediction. Journal of the Royal Statistical Society: Series B (Methodological), 40(1):1-24, 1978.
Yarema Okhrin, Wolfgang Schmid, and Ivan Semeniuk. New approaches for monitoring image data. IEEE Transactions Image Processing, 30:921-933, 2020.

---

Under review as a conference paper at ICLR 2023
Eulogio Pardo-Iguzquiza and Mario Chica-Olmo. Geostatisties with the matern semivariogram model: library of computer programs for inference, kriging and simulation. Computers & Geosciences, 34(9):1073-1079, 2008.
Joaquin Quinonero-Candela and Carl Edward Rasmussen. A unifying view of sparse approximate gaussian process regression. The Journal of Machine Learning Research, 6:1939-19955 2005.
Ali Rahimi and Benjamin Recht. Random features for large-scale kernel machines. Advances in neural information processing systems, 20, 2007.
Carl Edward Rasmussen. Gaussian processes in machine learning. In Summer school on machine learning, pp. 63-71. Springer, 2003.
Yunus Saatçi. Scalable inference for structured Gaussian process models. PhD thesis, Citeseer, 2012.
Thomas J Santner, Brian J Williams, William Notz, and Brain Williams, The design and analysis of computer experiments, volume 1. Springer, 2003.
Jasper Snoek, Hugo Larochelle, and Ryan P Adams. Practical bayesian optimization of machine learning algorithms. Advances in neural information processing systems, 25, 2012.
Niranjan Srinivas, Andreas Krause, Sham M Kakade, and Matthias Seeger. cess optimization in the bandit setting: No regret and experimental design. arXiv:0912.3995, 2009.
Gaussian pro- arXiv preprint
Michael L Stein. Interpolation of spatial data: some theory for kriging. Springer Science & Business Media, 1999.
Aretha L Teckentrup. Convergence of gaussian process regression with estimated and applications bayesian inverse problems. SIAM/ASA Journal on Uncertainty Quantification, 8(4):1310-1337, 2020.
William R Thompson. On the likelihood that one unknown probability exceeds another in view of the evidence of two samples. Biometrika, 25(3 (1:3-4):285-299 294, 1933.
Michalis Titsias. Variational learning of inducing variables in sparse gaussian processes. In Artificial intelligence and statistics, pp. 567-574. PMLR, 2009.
Christopher KI Williams and Carl Edward Rasmussen. Gaussian processes for machine learning, volume 2. MIT press Cambridge, MA, 2006.
Andrew Wilson and Hannes Nickisch. Kernel interpolation for scalable structured gaussian pro- cesses (kiss-gp). In International conference on machine learning, pp. 1775-1784 PMLR, 2015.
James Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Mare Deisenroth. Efficiently sampling functions from gaussian process posteriors. In International Conference Machine Learning. pp. 10292-10302. PMLR, 2020.
James Wilson, Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Mare Peter Deisenroth. Pathwise conditioning of gaussian processes. J. Mach. Learn. Res., 22:105-1 2021.
Sahar Zafari, Mariia Murashkina, Tuomas Eerola, Jouni Sampo, Heikki Kälviäinen, and Heikki Haario. Resolving overlapping convex objects in silhouette images by concavity analysis and gaussian process. Journal of Visual Communication and Image Representation, 73:102962, 2020.
12

---

Under review as a conference paper at ICLR 2023
A
APPENDIX
A.1
CONSTRUCTION OF BANDED MATRICES IN 7
Intermediate KPs Let a (01, .... ax) _ be vector with dn < ... dg, then intermediate KPs are defined as
(x(x) = Σ AzK(z,
(20)
and the coefficients Aj 's can be obtained by solving the (k - 1) x linear system
Σ
Asaz exp(őca;) O,
(21)
j=1
with O...., 32/2 and 8 +1.
One-sided KPs As before, let (aI, a²)² be vector with ai < ... < ag. one-sided KP given by
(x( az)(
(22)
with (k+ + 1)/2 < k: 1. For right-sided KPs, we can get coefficients Ay's by solving
Σ A,a,
= 0,
Σ ==1
Aja explcay =0
(23)
where 0, ..., (k- - 32/2 and second term equation 23 comprises auxiliary equations for
the case (k + 3)/2 with r 0...., (k + 3)/2. Similar to equation 21 equation 23 is an (s 1) x 8 linear system. Left-sided KPs are constructed similarly by solving the following equations:
A,a exp(ca,} = 0, Σ 1=1
Σ 1=1
aja exp{-ca,} = O,
(24)
where / = .,..., (k + 3)/2 with O,..., - (k + 32/2.
. (k: 32/2 and the second term comprises auxiliary equations for the case
>
KP Basis Let X
{ be the miedimmesioonal input data satisfying ž1 < ... En and K
a Matérn correlation function with half-integer smoothness v. Suppose n 2 k := 2y + 2. We can construct the following n functions, subset of linear space K. := span K (,
1. d1. .2..... Ok-- 1)/2. defined as left-sided KPs (( 0 (a d (1.....4-1)>>
11122//1111111111111 13/2111)
2. (1++1//2 01R+1)/211 (1-/221) defined as KPs ((e)...e)) (ep.... ( ".)
3. (n- Pn-1, On. ((5--(0-1)/2-1
defined as right
KPs
Here functions (0, im linearly independent in K., together with the fact that the dimension of K. is n, implies that ④ 1 forms basis for K,, referred to the KP basis.
13

---

Under review as a conference paper at ICLR 2023
In equation 7 the (1. entry of ((X) is (z(2)) In view of the compact supportedness of Øj Φ (X) is banded matrix with bandwidth (k: -3) -3)/2:
2 (It 1
((X)
k-a (zz))
@ j+ k- (Iz)
0 -يا (sr (1)+2
The matrix of Ax consists of the coefficients to construct the KPs. Ax is a banded matrix with bandwidth (h: - 1)/2:
A J- 2
Ax
Ajs
A
A
A.2 SAMPLING FROM MULTI-DIMENSIONAL GP POSTERIORS
Combine ① and 22) with the properties of the Kronecker product, the posterior mean and variance
for a training set [XPG, PG with full grid designs FG xd x can be computed via the following equations:
FG MMmm 5. K (Axx)((())))) fG - HG))
(25)
KFG o² Kr K, K,n X (Axx)(((()))))))
(26)
4,0/n
A.3 KERNEL PACKETS AND STATE-SPACE
Despite kernel packets (Chen et al. 2022 and state-space formulation (Grigorievskiy et al. 2017) can train one-dimensional Gaussian process regression with Matérn correlations exactly linear time, there is no obvious clue to establish any mathematical connections between these methods. It is worth noting that the state-space approach a sequential algorithm, which needs significantly more efforts to implement parallel computing; in contrast, the kernel packets approach can be paralleled trivial way. Besides, like most other supervised learning algorithms, the kernel packets method able to separate the learning task into two independent steps: training and prediction, and the computational time for each prediction point is only O(1) when the inputs are uniformly spaced. However, the space formulation cannot make this separation (Grigorievskly al. (017). and users have to either provide the prediction points during the training step, or cost ((n) time for each prediction point.
14

---

Under review as a conference paper at ICLR 2023
Log of Time
2-Wasserstein
5.0
RFF
Cholesky
2,5
0.0
-2.
KP
-5.0
10-
RFF
Cholesky
0.5
-7.5
10a
5-
10t
5- 10a 5- 10
1
caso
C200
Figure 7: Time and accuracy of different algorithms for sampling from one-dimensional GP priors with Matérn 5/2. Left: Logarithm of time taken to generate draw from GP priors, 22-axis is the number grid points p. Right: 2- 2-Wasserstein distances between priors and true distributions over points for three different cases {2+ )30 {2s (=2000 {24 when p 500.
Log of Time
Log of MSE
KP
KP
Decoupled
Decoupled
Cholesky
Cholesky
-10
-12
-14
5-10
10a 5-10a
10
5-10
10a
Figure 8: Time and accuracy of different algorithms for sampling from one-dimensional GP poste- riors with Matérn 5/2. Left: Logarithm time taken to generate draw from GP posteriors over 1000 points, z-axis is the number of observations Right: Logarithm of MSE over 1000 points.
A.4 NUMERICAL INSTABILITY
In real computation, the matrix Rz in Algorithm ! may not be numerically symmetric positive- definite when the distances between the input points are very small. Since each column of the matrix Az is obtained by solving null space, each column of Az can be up to a scalar, more work needs to be done to enhance the numerical stability the algorithm.
A.5 FIGURES
Figure and Figure show the performance of mie-dmmeniional examples with Matérn 5/2 in section Figure 9 and Figure I0 show the performance of two-dimensional examples with Matérn 5/2 in section
15

---

Under review as a conference paper at ICLR 2023
Log of Time
2 2-Wasserstein
KP
0.30
Cholesky
0.25
0,20
0.10
Cholesky
2
2
0.05
Figure 9: Time and accuracy of different algorithms for sampling from GP priors over full grids with Matérn 5/2. Left: Logarithm of time taken to generate a draw from GP priors. Right: 2- Wasserstein distances between priors and true distributions over nine points for three different cases
x x== f {2 when { = (-5 2 ,5 - 5-10- =].
Log Time
Log MSE
Decoupled
-5
Decoupled Cholesky
Cholesky
-10
-15
-20
28
2350
212
244
236
238
2a
2a
2.0
232
234
2
2
Figure 10: Time and accuracy of different algorithms for sampling from GP posteriors over full grids with Matérn 5/2. Left: Logarithm of time taken to generate draw from GP posteriors over 1024 points. Right: Logarithm of MSE over 1024 points.
16


Reviews:
[{"review_id": "sMhsDD75_g", "reviewer": "Reviewer_fxzb", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "The kernel packet approach allows fast solution for kernel systems with a 1D Matern kernel by reducing linear algebra with the dense kernel matrix to computations with a pair of dense matrices.  On tensor product grids, the tensor product of 1D Matern kernels can be treated the same way.  This approach can be used for fast inference (previous work) or fast sampling (the current paper).", "summary_of_the_review": "The paper takes the previously-introduced kernel packet decomposition for 1D Matern kernels and applies it not to the problem of training, but rather to the problem of generating sample draws.  The paper is well-written and easy to understand.  Much of the novelty is the previously-introduced kernel packet factorization, but that factorization has not been used in the context of sampling before.  The limitation to tensor products of 1D Matern kernels on tensor product grids is significant, though there are certainly settings where this combination is very natural.", "strength_and_weaknesses": "This seems like a reasonable algorithm for fast sampling at large numbers of points, though the restriction to tensor products of 1D Matern kernels (and tensor product grids in the multi-dimensional case) is significant. \n\nThe comparisons here are to methods like Random Fourier Features and dense methods.  1D Matern matrices are also very amenable to fast solves with rank-structured matrix techniques (using HSS solvers, HODLR solvers, etc), and these might make a more interesting point of comparison.  Such solvers are available online from Xia, Martinsson, Chandresekaran, and others.\n\nThe authors mention a numerical stability issue.  To the extent that they understand the source of this issue, it would be good to clarify what is happening.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "Not applicable"}, "extra_text": {}}, {"review_id": "g1YkAqkscE", "reviewer": "Reviewer_pohu", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "The manuscript proposes a scalable algorithm to draw sample functions from a univariate Gaussian process prior and posterior under the constraint that a Matérn covariance and a Gaussian likelihood is used.\nAlgorithmically, the method uses the decomposition of the covariance matrix into a product of two banded matrices whose bandwidth depends on the smoothness of the underlying GP. A set of comparative experiments illustrates some properties of the method.", "summary_of_the_review": "Even though, the \"kernel packet\" view has not been used for sampling before, the manuscript needs to be improved before getting published. Exact computations and state space computations need to be included in the experiments. The relation between the state space representation needs to be better analysed. Numerical issues are not yet well understood.", "strength_and_weaknesses": " + The manuscript is rather well written and simple to understand.\n + The manuscript studies properties of a widely used and well-understood model and provides a scalable sampling algoritm.\n - The method is only applicable to a very small subset of GP models (1D, Matérn covariance, Gaussian regression, rather small smoothness parameter).\n - There is no connection made to the state space representation of GPs. In fact, a scalable algorithm with exactly the same space/time complexity can be obtained via the state space view (see SpInGP https://arxiv.org/abs/1610.08035 and SSGP https://arxiv.org/abs/1802.04846 for two papers where explicit sparse matrices are used). The SpInGP needs to be included as a baseline in the experiments.\n - The manuscript does not make a step forward to better understand the theory behind the \"kernel packets\" and does not provide insights to allow a statement whether the \"kernel packets\" are a rediscovery of the state space representation from a different angle or whether the \"kernel packets\" are something different or more generally applicable beyond the Matérn covariance. Even, numerical comparisons could have given insights.\n - It is well known that SpInGP comes with numerical instabilities. Also the manuscript mentions numerical issues but does not provide a detailed and insightful theoretical or empirical analysis.\n - The numerical experiments do not cover large scale datasets.\n - A code base to replicate the experiments is missing.", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "OP-zWUs9j5", "reviewer": "Reviewer_4pJF", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "The authors consider the problem of sampling from Gaussian processes defined using the Matérn kernel (both the prior and posterior) in one-dimension. Building on the \"kernel packet\" approach of Chen et al, 2022, the authors show that samples can be drawn in time that is nearly linear in the number of points sampled. Further, these samples are \"exact\" in the sense that any error is due to machine precision. Samples from the posterior can be drawn by combining this approach with Matheron's rule, similar to as in Wilson et al, 2020. The approach can be extended to several dimensions if products of one-dimensional Matérn kernels are considered and points are sampled on a grid. The authors provide several small scale experiments to validate the proposed approach.", "summary_of_the_review": "While narrow in scope, the method seems close to being practical for an interesting task, which is a strong indicator in favor of acceptance. However, I would like to see this practicality demonstrated more convincingly in experiments. The paper falls short in this regard, and doesn't appear to provide sufficient detail to reproduce them (particularly for baseline methods). Crucially, a method must be numerically stable as well as fast to be practical, and I think the paper would be strengthened by some discussion of why there is often such a large discrepancy between the method presented and the Cholesky based method, which is also \"exact\" ignoring numerical issues. \n\nIt is also theoretically interesting that samples can be drawn from these GP priors/posteriors in (nearly) linear time, although I think most of the (theoretical) novelty in this regard comes from Chen et al 2022.\n\nI am currently leaning toward rejection of the paper on balance, but I think it is borderline\n", "strength_and_weaknesses": "## Strengths\n- The proposed idea seems practical.\n- Even in one-dimension (where things are inevitably simpler) achieving (near) linear time exact sampling is quite desirable. While the method is quite specific, Matérn kernels are at least widely used as the authors point out.\n\n## Weaknesses\n- The contribution beyond work contained in Chen, 2022 is reasonably small. A nice linear-algebraic observation is needed to obtain samples (section 3.1), but I am unsure if this on its own is convincing for a paper without showing the method more convincingly in an interesting application (more discussion of experiments below).\n- The scope of the method is narrow. The limitation to Matérn kernels on its own is already quite narrow. This is made more restrictive by the assumption of one-dimensional data. Note also that the multi-dimensional Matérn kernel that is generally used is not a product of one-dimensional  Matérn kernels (this should be made explicit in the text). The product of Matérn kernel considered does not have the same smoothness properties generally associated to Matérn kernels, which seems to be a disadvantage. A bit more explanation in the text as to why the method cannot be generalized to a wider class of kernels would be useful.\n- The experiments are not described as thoroughly as they should be, and could certainly be strengthened, more specific comments are below.\n", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "KM36Gb0wRZr", "reviewer": "Reviewer_xLQn", "summary": "", "questions": "", "limitations": "", "rating": 3, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "The paper proposes an exact and scalable sampling algorithm for Gaussian processes when using Matern correlation functions. The authors make use of the recent kernel packets formalism that enables a sparse representation of the covariance matrix and reduced computational complexity. Numerical results that show the advantages of the proposed algorithm complete the paper. ", "summary_of_the_review": "The paper proposes a sampling algorithm derived as mix of a recent theoretical discovery, Kernel Packets, and a recently re-discovered rule, Matheron's rule. The novelty factor of the paper is low and there are some issues with the numerical experiments. Leaning towards reject.", "strength_and_weaknesses": "The proposed algorithm is interesting as it enables efficient sampling in the case of Matern correlation functions. The paper can be overall seen as a mix of the results in [Chen et al.] and [Wilson et al.]: kernel packets from [Chen et al.], Matheron's rule from [Wilson et al.]. There's nothing inherently wrong in doing so, however, the novelty factor is low.\n\nA minor weakness is that not enough credit is given to the two source papers. Take for example the numerical applications, why was the Griewank function chosen section 4.2 ? Or why was the Thompson sampling considered in section 4.3 ? There is nothing wrong in considering the same examples as other papers, more so as you compare to them. There is nothing wrong with a sentence the likes of \"we consider the same application as in [insert paper citation here]\".\n\nA moderate weakness is the lack of interpretation of the numerical results. You content most of the times in just pointing out that the results are better. When it comes for computational cost, there is no issue in interpreting the results. However, when it comes to error analysis, it is not that obvious or not all the implications are obvious. Consider for example the results in figure 1. For case 1 and 3 you have a lower 2-Wasserstein distance than Cholesky decomposition. For me it is not clear why an exact method, your approach, is better than another exact method, Cholesky decomposition approach. Such an explanation would have been more than welcome. You also mention in the conclusion part about some numerical stability issues. It would have been beneficial to mention those issues when the results were presented. It would have been even better to give at least a hint at the cause.\n\n[Chen et al.] - \"Kernel Packet: An Exact and Scalable Algorithm for Gaussian Process Regression with Matérn Correlations\", Journal of Machine Learning Research 23 (2022) 1-32\n\n[Wilson et al.] - \"Efficiently sampling functions from gaussian process posteriors\", Proceedings of the 37 th International Conference on Machine Learning, Online, PMLR 119, 2020", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Decision: Reject