## Gaussian Processes Sampling with Sparse Grids under Additive Schwarz Preconditioner

Haoyuan Chen 1 ∗ Rui Tuo 1

1 Department of Industrial and Systems Engineering, Texas A&amp;M University College Station, TX 77843 chenhaoyuan2018@tamu.edu ruituo@tamu.edu

## Abstract

Gaussian processes (GPs) are widely used in non-parametric Bayesian modeling, and play an important role in various statistical and machine learning applications. In a variety tasks of uncertainty quantification, generating random sample paths of GPs is of interest. As GP sampling requires generating high-dimensional Gaussian random vectors, it is computationally challenging if a direct method, such as the Cholesky decomposition, is used. In this paper, we propose a scalable algorithm for sampling random realizations of the prior and posterior of GP models. The proposed algorithm leverages inducing points approximation with sparse grids, as well as additive Schwarz preconditioners, which reduce computational complexity, and ensure fast convergence. We demonstrate the efficacy and accuracy of the proposed method through a series of experiments and comparisons with other recent works.

KeywordsGaussian processes, sparse grids, additive Schwarz preconditioner

## 1 Introduction

Gaussian processes (GPs) hold significant importance in various fields, particularly in statistical and machine learning applications (Rasmussen, 2003; Santner et al., 2003; Cressie, 2015), due to their unique properties and capabilities. As probabilistic models, GPs are highly flexible and able to adapt to a wide range of data patterns, including those with complex, nonlinear relationships. GPs have been successfully applied in numerous domains including regression (O'Hagan, 1978; Bishop et al., 1995; MacKay et al., 2003), classification (Kuss et al., 2005; Nickisch and Rasmussen, 2008; Hensman et al., 2015), time series analysis (Girard et al., 2002; Roberts et al., 2013), spatial data analysis (Banerjee et al., 2008; Nychka et al., 2015), Bayesian networks (Neal, 2012), optimization (Srinivas et al., 2009), and more. This wide applicability showcases their versatility in handling different types of problems. As the size of the dataset increases, GP inference suffers from computational issues of the operations on the covariance matrix such as inversion and determinant calculation. These computational challenges have motivated the development of various approximation methods to make GPs more scalable. Recent advances in scalable GP regression include Nyström approximation (Quinonero-Candela and Rasmussen, 2005; Titsias, 2009; Hensman et al., 2013), random Fourier features (Rahimi and Recht, 2007; Hensman et al., 2018), local approximation (Gramacy and Apley, 2015; Cole et al., 2021), structured kernel interpolation (Wilson and Nickisch, 2015; Stanton et al., 2021), state-space formulation (Hartikainen and Särkkä, 2010; Grigorievskiy et al., 2017; Nickisch et al., 2018), Krylov methods (Gardner et al., 2018; Wang et al., 2019), Vecchia approximation (Katzfuss and Guinness, 2021; Katzfuss et al., 2022), sparse representation (Ding et al., 2021; Chen et al., 2022; Chen and Tuo, 2022), etc.

The focus of this article is on generating samples from GPs, which includes sampling from both prior and posterior processes. GP sampling is a valuable tool in statistical modeling and machine learning, and is used in various applications, such as in Bayesian optimization (Snoek et al., 2012; Frazier, 2018b,a), in computer

∗ Submitted to the editors on January 31, 2024.

model calibration (Kennedy and O'Hagan, 2001; Tuo and Wu, 2015; Plumlee, 2017; Xie and Xu, 2021), in reinforcement learning (Kuss and Rasmussen, 2003; Engel et al., 2005; Grande et al., 2014), and in uncertainty quantification (Murray-Smith and Pearlmutter, 2004; Marzouk and Najm, 2009; Teckentrup, 2020). Practically, sampling from GPs involves generating random vectors from a multivariate normal distribution after discretizing the input space. However, sampling from GPs poses computational challenges particularly when dealing with large datasets due to the need to factorize the large covariance matrices, which is an O ( n 3 ) operation for n inputs. There has been much research on scalable GP regression, yet sampling methodologies remain relatively unexplored.

Scalable sampling techniques for GPs are relatively rare in the existing literature. A notable recent development in the field is the work of Wilson et al. (2020), who introduced an efficient sampling method known as decoupled sampling. This approach utilizes Matheron's rule and combines the inducing points approximation with random Fourier features. Furthermore, Wilson et al. (2021) expanded this concept to pathwise conditioning based on Matheron's update. This method, which requires only sampling from GP priors, has proven to be a powerful tool for both understanding and implementing GPs. Inspired by this work, Maddox et al. (2021) applied Matheron's rule to multi-task GPs for Bayesian optimization. Additionally, Nguyen et al. (2021) applied such techniques to tighten the variance and enhance GP posterior sampling. Wenger et al. (2022b) introduced IterGP, a novel method that ensures consistent estimation of combined uncertainties and facilitates efficient GP sampling by computing a low-rank approximation of the covariance matrix. Lin et al. (2023) developed a method for generating GP posterior samples using stochastic gradient descent, targeting the complementary problem of approximating these samples in O ( n ) time. It is important to note, however, that none of these methods have explored sampling within a structured design, an area that possesses many inherent properties.

In this study, we propose a novel algorithms for efficient sampling from GP priors and GP regression posteriors by utilizing inducing point approximations on sparse grids. We show that the proposed method can reduce the time complexity to O ( ( η + n n ) sg ) for input data of size n and inducing points on sparse grids of level η and size n sg . This results in a linear-time sampling process, particularly effective when the size of the sparse grids is maintained at a reasonably low level. Specifically, our work makes the following contributions:

- · We propose a novel and scalable sampling approach for Gaussian Processes using product kernels on multi-dimensional points. This approach leverages the inducing point approximation on sparse grids, and notably, its computational time scales linearly with the product of the sizes of the inducing points and the sparse grids.
- · Our method provides a 2 -Wasserstein distance error bound with a convergence rate of O ( n -ν sg (log n sg ) ( ν +2)( d -1)+ +1 d ) for d -dimensional sparse grids of size n sg and product Matérn kernel of smoothness ν .
- · We develop a two-level additive Schwarz preconditioner to accelerate the convergence of conjugate gradient during posterior sampling from GPs.
- · We illustrate the effectiveness of our algorithm in real-world scenarios by applying it to areas such as Bayesian optimization and dynamical system problems.

The structure of the paper is as follows: Background information is presented in Section 2, our methodology is detailed in Section 3, experimental results are discussed in Section 4, and the paper concludes with discussions in Section 5. Detailed theorems along with their proofs are provided in Appendix A, and supplementary plots can be found in Appendix B.

## 2 Background

This section provides an overview of the foundational concepts relevant to our proposed method. Section 2.1 introduces GPs, GP regression, and the fundamental methodology for GP sampling. The underlying assumptions of our work are outlined in Section 2.2. Additionally, detailed explanations of the inducing points approximation and sparse grids, which are key components of our sampling method, are presented in Section 2.3 and Section 2.4, respectively.

## 2.1 GP review

GPs AGPis a collection of random variables uniquely specified by its mean function µ ( ) · and kernel function K , ( · · ) , such that any finite collection of those random variables have a joint Gaussian distribution (Williams and Rasmussen, 2006). We denote the corresponding GP as GP ( µ ( ) · , K ( · , · )) .

GP regression Let f : R d → R be an unknown function. Suppose we are given a set of training points X = { x i } n i =1 and the observations y = ( y , . . . , y 1 n ) ⊤ where y i = f ( x i ) + ϵ i with the i.i.d. noise ϵ i ∼ N (0 , σ 2 ϵ ) .

The GP regression imposes a GP prior over the latent function as f ( ) · ∼ GP ( µ ( ) · , K ( · , · )) . Then the posterior distribution at m test points X ∗ = { x ∗ i } m i =1 is

$$p ( f ( \mathbf X ^ { * } ) | \mathbf y ) = \mathcal { N } ( \mu _ { * | \mathbf y }, \mathbf K _ { *, * | \mathbf y } ),$$

with

$$\mu _ { * | y } & = \mu _ { * } + \mathbf K _ { *, X } \left [ \mathbf K _ { X, X } + \sigma _ { \epsilon } ^ { 2 } \mathbf I _ { n } \right ] ^ { - 1 } \left ( y - \mu _ { X } \right ), & ( 2 a )$$

$$K _ { *, * | y } & = K _ { *, * } - K _ { *, X } \left [ K _ { X, X } + \sigma _ { \epsilon } ^ { 2 } \mathfrak { I } _ { n } \right ] ^ { - 1 } K _ { X, * }, & ( 2 b )$$

where σ 2 ϵ &gt; 0 is the noise variance, I n is a n × n identity matrix, K ∗ , X = K ( X ∗ , X ) = [ K ( x ∗ i , x j ) ] i = m,j = n i =1 ,j =1 = K ( X X , ∗ ) ⊤ = K ⊤ X , ∗ , K X X , = [ K ( x i , x j ) ] n i,j =1 , K ∗ ∗ , = [ K ( x ∗ i , x ∗ j ) ] m i,j =1 , µ X = ( µ ( x 1 ) , · · · , µ ( x n ) ) ⊤ and µ ∗ = ( µ ( x ∗ 1 ) , · · · , µ ( x ∗ m ) ) ⊤ .

Sampling from GPs The goal is to sample a function f ( ) · ∼ GP ( µ ( ) · , K ( · , · )) . To represent this in a finite form, we discretize the input space and focus on the function values at a set of grid points Z = { z i } n s i =1 . Our goal is to generate samples f Z = ( f ( z 1 ) , · · · , f ( z n s ) ) ⊤ from the multivariate normal distribution N ( µ ( Z ) , K ( Z Z , )) = N ( µ Z , K Z Z , ) . The standard method for sampling from a multivariate normal distribution involves the following steps: (1) generate a vector of samples f I whose entries are independent and identically distributed normal, i.e. f I ns ∼ N ( 0 I , n s ) ; (2) Use the Cholesky decomposition (Golub and Van Loan, 2013) to factorize the covariance matrix K Z Z , as L L Z ⊤ Z = K Z Z , ; (3) generate the output sample f Z as

$$f _ { Z } \leftarrow L _ { Z } f _ { I _ { n _ { s } } } + \mu _ { Z }.$$

Sampling a posterior GP follows a similar procedure. Suppose we have observations y = ( y , 1 · · · , y n ) ⊤ on n distinct points X = { x i } n i =1 , where y i = f ( x i ) + ϵ i with noise ϵ i ∼ N (0 , σ 2 ϵ ) and f ( ) · ∼ GP ( µ ( ) · , K ( · , · )) . The goal is to generate posterior samples f ( X ∗ ) | y at m test points X ∗ = { x ∗ i } m i =1 . Since the posterior samples f ( X ∗ ) | y ∼ N ( µ ∗| y , K ∗ ∗| , y ) also follow a multivariate normal distribution, as detailed in eq. (2a) and eq. (2b) in Section 2.1, we can apply Cholesky decomposition L ∗| X L ⊤ ∗| X = K ∗ ∗| , X to generate GP posterior sample f ∗ = f ( X ∗ ) = ( f ( x ∗ 1 ) , · · · , f ( x ∗ m ) ) ⊤ as

$$f _ { * } \, | \, y \leftarrow L _ { * | X } f _ { I _ { n } } + \mu _ { * | X },$$

where f I n ∼ N ( 0 I , n ) .

## 2.2 Assumptions

A widely used kernel structure in multi-dimensional settings is 'separable' or 'product' kernels, which can be expressed as follows:

$$K ( { \mathbf x }, { \mathbf x } ^ { \prime } ) = \sigma ^ { 2 } \prod _ { j = 1 } ^ { d } K _ { 0 } ^ { ( j ) } ( { \mathbf x } ^ { ( j ) }, { \mathbf x } ^ { \prime ( j ) } ),$$

where x ( j ) and x ′ ( j ) are the j -th components of x and x ′ respectively, K ( j ) 0 ( · , · ) denotes the one-dimensional correlation function for each dimension j with the variance set to one, d denotes the dimension of the input points. While the product of one-dimensional kernels does not share identical smoothness properties with multidimensional kernels, the assumption of separability is widely adopted in the field of spatio-temporal statistics (Gneiting et al., 2006; Genton, 2007; Constantinou et al., 2017). Its widespread adoption is mainly due to its efficacy in enabling the development of valid space-time parametric models and in simplifying computational processes, particularly when handling large datasets for purposes such as inference and parameter estimation.

In this article, we consider the inputs to be d -dimensional points X = { x i } n i =1 where each x i ∈ R d and d ∈ N + , and we focus on product Matérn kernels that employ identical base kernels across every dimension. Matérn kernel (Genton, 2001) is a popular choice of covariance functions in spatial statistics (Diggle et al., 2003), geostatistics (Curriero, 2006; Pardo-Iguzquiza and Chica-Olmo, 2008), image analysis (Zafari et al., 2020; Okhrin et al., 2020), and other applications. The Matérn kernel is defined as the following form (Stein, 1999):

$$K ( x, x ^ { \prime } ) = \sigma ^ { 2 } \frac { 2 ^ { 1 - \nu } } { \Gamma ( \nu ) } \left ( \sqrt { 2 \nu } \frac { | x - x ^ { \prime } | } { \omega } \right ) ^ { \nu } K _ { \nu } \left ( \sqrt { 2 \nu } \frac { | x - x ^ { \prime } | } { \omega } \right ),$$

for any x, x ′ ∈ R , where σ 2 &gt; 0 is the variance, ν &gt; 0 is the smoothness parameter, ω &gt; 0 is the lengthscale and K ν is the modified Bessel function of the second kind.

## 2.3 Inducing points approximation

Williams and Seeger (2000) first applied Nyström approximation to kernel-based methods such as GPs. Following this, Quinonero-Candela and Rasmussen (2005) presented a novel unifying perspective, and reformulated the prior by introducing an additional set of latent variables f u , known as inducing variables . These latent variables correspond to a specific set of input locations u = ( u , 1 · · · , u q ) ⊤ , referred to as inducing inputs . In the subset of regressors (SoR) (Silverman, 1985) algorithm, the kernels can be approximated using the following expression:

$$K _ { S O R } ( x _ { i }, x _ { j } ) = K ( x _ { i }, u ) K ( u, u ) ^ { - 1 } K ( u, x _ { j } ).$$

Employing the SoR approximation detailed in eq. (7), the form of the SoR predictive distribution q SoR ( f ∗ | y ) can be derived as follows:

$$q _ { \text{SoR} } ( \mathfrak { f } _ { * } \, | \, y ) = \mathcal { N } ( \mu _ { * } + \sigma _ { \epsilon } ^ { - 2 } \mathbf K _ { *, u } \Sigma _ { u } ^ { - 1 } \mathbf K _ { u, x } ( y - \mu _ { x } ), \, \mathbf K _ { *, u } \Sigma _ { u } ^ { - 1 } \mathbf K _ { u, * } ), \, \quad ( 8 )$$

where we define Σ u = [ K u u , + σ -2 ϵ K u X , K X u , ] , K ∗ , u = K ( X ∗ , u ) = [ K ( u X , ∗ ) ] ⊤ = K ⊤ ∗ u , , K u u , = [ K ( u i , u s ) ] q i,s =1 , µ X = ( µ ( x 1 ) , · · · , µ ( x n ) ) ⊤ , and µ ∗ = ( µ ( x ∗ 1 ) , · · · , µ ( x ∗ m ) ) ⊤ . Consequently, the computation of GP regression using inducing points can be achieved with a time complexity of O ( nq 2 ) for q inducing points.

To quantify the approximation power in the inducing points approach, we shall use the Wasserstein metric . Consider a separable complete metric space denoted as ( M,d ) . The p -Wasserstein distance between two probability measures µ 1 and µ 2 , on the space M , is defined as follows:

$$W _ { p } ( \mu, \nu ) = \left ( \inf _ { \gamma \in \Gamma ( \mu _ { 1 }, \mu _ { 2 } ) } \mathbb { E } _ { ( x, y ) \sim \gamma } d ^ { p } ( x, y ) \right ) ^ { 1 / p },$$

for p ∈ [1 , + ∞ ) , where Γ( µ , µ 1 2 ) is the set of joint measures on M M × with marginals µ 1 and µ 2 . Theorem A.1 in Appendix A shows that the inducing points approximation of a GP converges rapidly to the target GP, with a convergence rate of O ( n -ν ) for the Matérn kernels with smoothness ν , provided that the inducing points are uniformly distributed. It is worth noting that, inequality presented in eq. (30) is also applicable to multivariate GPs.

Therefore, the convergence rate for multivariate inducing points approximation can be established, assuming the corresponding convergence rate of GP regression is known. A general framework to investigate the latter problem is given by Wang et al. (2020). Inspired by Theorem A.1, we assume that the vector y = ( y , 1 · · · , y n ) ⊤ is observed on n distinct d -dimensional points X = { x i } n i =1 . The inducing inputs U , as defined in Section 2.3, are given by a full grid, denoted as U = × d j =1 u ( j ) , where each u ( j ) is a set of q j one-dimensional data points,

A × B denotes the Cartesian product of sets A and B , accordingly the grid U has q = ∏ d j =1 q j points. Given the separable structure as defined in eq. (5) and utilizing the full grid U , we can transform the multi-dimensional problem to a one-dimensional problem since the covariance matrix K U U , : = K ( U U , ) can be represented by Kronecker products (Henderson et al., 1983; Saatçi, 2012; Wilson and Nickisch, 2015) of matrices over each input dimension j :

$$K _ { U, U } = K ( U, U ) = \bigotimes _ { j = 1 } ^ { d } K _ { j } ( u ^ { ( j ) }, u ^ { ( j ) } ).$$

However, the size of the inducing points q grows exponentially as d gets large, therefore we consider using sparse grids (Garcke, 2013; Rieger and Wendland, 2017) to solve the curse of dimensionality.

## 2.4 Sparse grids

A sparse grid is built on a nested sequence of one-dimensional points ∅ = U j, 0 ⊆ . . . ⊆ U j,t ⊆ U j,t +1 ⊆ . . . for each dimension j , j = 1 , . . . , d . A sparse grid with level η and dimension d is a design as follows:

$$\mathcal { U } ( \eta, d ) = \bigcup _ { \vec { t } \in G ( \eta ) } \mathcal { U } _ { 1, t _ { 1 } } \times \mathcal { U } _ { 2, t _ { 2 } } \times \cdots \times \mathcal { U } _ { d, t _ { d } },$$

where # » t = ( t 1 , . . . , t d ) , G ( η ) = { # » t ∈ N d | ∑ d j =1 t j = η } . Figure 1 and Figure 9 show sparse grids corresponding to the hyperbolic cross points (bisection) for two dimensional points at levels η = 3 4 5 6 , , , and for three dimensional points at η = 4 5 6 7 , , , respectively. Hyperbolic cross points sets, defined over the interval [0 , 1] , are detailed in eq. (11):

$$\psi \circ \psi ^ { \prime } \quad u _ { j, t _ { j } } = \left \{ \frac { i } { 2 ^ { t _ { j } } } \, \colon i = 1, \dots, 2 ^ { t _ { j } } - 1 \right \}.$$

An operator on the sparse grid can be efficiently computed using Smolyak's algorithm (Smolyak, 1963; Ullrich, 2008), which states that the predictor operator P ( η, d ) on function f with respect to the sparse grid U ( η, d ) can

Figure 1: Sparse grids U ( η, d ) over [0 , 1] d of level η = 3 4 5 6 , , , and dimension d = 2 .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-11854330-4bac-4971-b3ea-48ae26a5cb95.jpg)

be evaluated by the following formula:

$$\mathcal { P } ( \eta, d ) ( f ) = \sum _ { \vec { t } \in \mathbb { P } ( \eta, d ) } ( - 1 ) ^ { \eta - | \vec { t } | } \binom { d - 1 } { \eta - | \vec { t } | } \left ( \bigotimes _ { j = 1 } ^ { d } u _ { j, t _ { j } } \right ) ( f ),$$

where | # » t | = ∑ d j =1 t j , P ( η, d ) = { # » t ∈ N d | max( d, η -d +1) ≤ | # » t | ≤ η } , ( ⊗ d j =1 U j,t j ) here represents the operator over the Kronecker product points ⊗ d j =1 U j,t j . Plumlee (2014) and Yadav et al. (2022) also provided the linear time computing method for GP regression based on Smolyak's algorithm. Theorem A.3 in Appendix A demonstrates that for a d -dimensional GP, the inducing points approximation converges to the target GP at the rate of O ( n -ν (log n ) ( ν +2)( d -1)+ +1 d ) for the product Matérn kernels with smoothness ν , given that the inducing points are sparse grids U ( η, d ) of level η and dimension d , as defined in eq. (10).

## 3 Sampling with sparse grids

In this section, we elaborate the algorithm for sampling from multi-dimensional GPs for product kernels with sparse grids.

## 3.1 Prior

Sampling algorithm Our goal is to sample from GP priors f ( ) · ∼ GP ( µ ( ) · , K ( · , · )) over d -dimensional grid points Z = { z i } n s i =1 , z i ∈ R d . According to the conditional distributions of the SoR model in equations eq. (7), the distribution of f Z = f ( Z ) = ( f ( z 1 ) , · · · , f ( z n s )) ⊤ can be approximated by K Z , U K -1 U U , f U with f U ∼ N ( 0 K , U U , ) , where f U = f ( U ) , and U = U ( η, d ) is a sparse grid design with level η and dimension d defined in eq. (10). Then it reduces to generating samplings from K -1 U U , f U ∼ N ( 0 K , -1 U U , ) . Apply Smolyak's algorithm to samplings from GPs (see Theorem A.4 in Appendix A), we can infer that

$$\mathbf K _ { \mathcal { U }, \mathcal { U } } ^ { - 1 } \mathfrak { f } _ { \mathcal { U } } \stackrel { d } { = } \sum _ { \vec { t } \in \mathcal { P } ( \eta, d ) } ( - 1 ) ^ { \eta - | \vec { t } | } \binom { d - 1 } { \eta - | \vec { t } | } \bigotimes _ { j = 1 } ^ { d } \text{chol} \left ( K ( \mathcal { U } _ { j, t _ { j } }, \mathcal { U } _ { j, t _ { j } } ) ^ { - 1 } \right ) \mathfrak { f } _ { n _ { \vec { t } } } \delta _ { x } ( \mathcal { U } _ { \vec { t } } ), \quad ( 1 3 )$$

$$\delta _ { x } ( \mathcal { U } _ { \vec { t } } ) = \begin{cases} 0 & \text{ if } x \in \mathcal { U } \, \mathcal { U } _ { \vec { t } } \\ 1 & \text{ if } x \in \mathcal { U } _ { \vec { t } } \end{cases},$$

where chol ( K ( U j,t j , U j,t j ) -1 ) is the Cholesky decomposition of K ( U j,t j , U j,t j ) -1 , U = U ( η, d ) is a sparse grid design with level η and dimension d defined in eq. (10), f I n t # » ∼ N ( 0 I , n # » t ) is standard multivariate normal distributed, I n # » t is a n # » t × n # » t identity matrix, n # » t = ∏ d j =1 t j is the size of the full grid design U # » t = U 1 ,t 1 ×U 2 ,t 2 ×··· × U d,t d associated with # » t = ( t 1 , . . . , t d ) , δ x ( U # » t ) is a Dirac measure on the set U . Therefore, the GP prior f Z over grid points Z can be approximately sampled via

$$f _ { z } \leftarrow K _ { z, u } \left ( \sum _ { \vec { t } \in P ( \eta, d ) } ( - 1 ) ^ { \eta - | \vec { t } | } \binom { d - 1 } { \eta - | \vec { t } | } \bigotimes _ { j = 1 } ^ { d } \text{chol} \left ( K ( \mathcal { U } _ { j, t _ { j } }, \mathcal { U } _ { j, t _ { j } } ) ^ { - 1 } \right ) f _ { n _ { \vec { t } } } \delta _ { z } ( \mathcal { U } _ { \vec { t } } ) \right ) \\ + \mu _ { z }, \quad ( 1 5 )$$

where f I n t # » ∼ N ( 0 I , n # » t ) is standard multivariate normal distributed, δ x ( U # » t ) is a Dirac measure on set U as defined in eq. (14), chol ( K ( U j,t j , U j,t j ) -1 ) is the Cholesky decomposition of K ( U j,t j , U j,t j ) -1 .

Complexity Since the Kronecker matrix-vector product can be written as the form of vectorization operator (see Proposition 3.7 in (Kolda, 2006)), the Kronecker matrix-vector product in eq. (15) only costs O ( ∑ d j =1 t j ∏ d j =1 t j ) = O ( ηn # » t ) time, hence the time complexity of the prior sampling via eq. (15) is O ( ( η + n s ) ∑ # » t ∈ P ( η,d ) n # » t ) = O ( ( η + n s ) n sg ) , where ν is the smoothness parameter of the Matérn kernel, d is the dimension of the sampling points Z = { z i } n s i =1 , n s is the size of the sampling points, n sg is the size of inducing points on a sparse grid U ( η, d ) , n # » t = ∏ d j =1 t j is the size of the full grid design U # » t = × d j =1 U j,t j , # » t = ( t 1 , . . . , t d ) , P ( η, d ) = { # » t ∈ N d | max( d, η -d +1) ≤ ∑ d j =1 t j ≤ η } .

Note that we can obtain the upper bound of the number of points in the Smolyak type sparse grids. According to Lemma 3 in (Rieger and Wendland, 2017), given that |U j,t j | ≤ 2 t j -1 , when sparse grids are constructed on hyperbolic cross points as defined in eq. (11), the cardinality of a sparse grid U ( η, d ) defined in eq. (10), is thereby bounded by

$$\text{vann} \quad n _ { s g } = | \mathcal { U } ( \eta, d ) | & \leq 1 ^ { d } ( \frac { 2 - 1 } { 2 } ) ^ { d } 2 ^ { \eta } \binom { \eta - 1 } { d - 1 } \min \{ \frac { 2 } { 2 - 1 }, \frac { \eta } { d } \} \\ & \leq 2 ^ { \eta - d } \frac { \eta ^ { d - 1 } } { ( d - 1 )! } \min \{ 2, \frac { \eta } { d } \} \\ & \leq 2 ^ { \eta - d } \eta ^ { d } / ( d! ).$$

## 3.2 Posterior

Matheron's rule Sampling from the posterior GPs can be effectively conducted using Matheron's rule . The Matheron's rule was popularized in the geostatistics field by Journel and Huijbregts (1976), and recently rediscovered by Wilson et al. (2020), who leveraged it to develop a novel approach for GP posterior samplings. Matheron's rule can be described as follows: Consider a and b as jointly Gaussian random vectors. Then the random vector a , conditional on b = β , is equal in distribution to

$$( \mathbf a \, | \, b = \beta ) \overset { d } { = } \mathbf a + C o v ( \mathbf a, b ) C o v ( b, b ) ^ { - 1 } ( \beta - b ),$$

where Cov( a b , ) denotes the covariance of ( a b , ) .

According to eq. (17), exact GP posterior samplers can be obtained using two jointly Gaussian random variables. Following this technique, we derive

$$f _ { * } \ | \, y \stackrel { d } { = } f _ { * } + K _ { *, x } \left [ K _ { X, X } + \sigma _ { \epsilon } ^ { 2 } \mathbf I _ { n } \right ] ^ { - 1 } ( y - f _ { X } - \epsilon ),$$

where f X = f ( X ) and f ∗ = f ( X ∗ ) are jointly Gaussian random variables from the prior distribution, noise variates ϵ ∼ N ( 0 , σ 2 ϵ I n ) . Clearly, the joint distribution of ( f X , f ∗ ) follows the multivariate normal distribution as follows:

$$( f _ { X }, f _ { * } ) \sim \mathcal { N } \left ( \begin{bmatrix} \mu _ { X } \\ \mu _ { * } \end{bmatrix}, \begin{bmatrix} K _ { X, X } & K _ { X, * } \\ K _ { *, X } & K _ { *, * } \end{bmatrix} \right ).$$

Sampling algorithm We can apply SoR approximation eq. (7) to Matheron's rule eq. (17) and obtain

$$f _ { * } \, | \, y \stackrel { d } { \approx } f _ { * } + \sigma _ { \epsilon } ^ { - 2 } K _ { *, 4 } \Sigma _ { u } ^ { - 1 } K _ { u, x } ( y - f _ { x } - \epsilon ),$$

where Σ U = [ K U U , + σ -2 ϵ K U , X K X , U ] , f X = f ( X ) and f ∗ = f ( X ∗ ) are jointly Gaussian random variables from the prior distribution, noise variates ϵ ∼ N ( 0 , σ 2 ϵ I n ) . The prior resulting from the SoR approximation can be easily derived from eq. (7), yielding the following expression:

$$q _ { S o R } ( f _ { X }, f _ { * } ) & = \mathcal { N } ( \mu _ { X + * }, \, K _ { X + *, u } K _ { u, u } ^ { - 1 } K _ { u, x + * } ), & ( 2 1 )$$

where µ X + ∗ = [ µ X , µ ∗ ] ⊤ , K X + ∗ U , = [ K ( U , X ) , K ( U , X ∗ )] ⊤ . The same as above in eq. (15), ( f X , f ∗ ) can be sampled from

$$( \mathfrak { f } _ { X }, \mathfrak { f } _ { * } ) \leftarrow K _ { X + *, \mathcal { U } } \left ( \sum _ { \mathfrak { r } \in \mathcal { P } ( \eta, d ) } ( - 1 ) ^ { \eta - | \mathfrak { r } | } \binom { d - 1 } { \eta - | \vec { t } | } \bigotimes _ { j = 1 } ^ { d } \text{chol} \left ( K ( \mathcal { U } _ { j, t _ { j } }, \mathcal { U } _ { j, t _ { j } } ) ^ { - 1 } \right ) \mathfrak { f } _ { n _ { \vec { r } } } \delta _ { x } ( \mathcal { U } _ { \vec { r } } ) \right )$$

+

µ

X

+

∗

,

where f I n t # » , δ x ( U # » t ) , and chol ( K ( U j,t j , U j,t j ) -1 ) are the same as that defined in eq. (15).

(22)

Algorithm 1 The Preconditioned Conjugate Gradient (PCG) Algorithm.

- 1: Input: matrix Σ , preconditioning matrix P , vector v , convergence threshold ϵ , initial vector x 0 , maximum number of iterations T
- 2: Output: x ∗ ≈ Σ -1 v
- 3: r 0 := v -Σ x 0 ; z 0 := P -1 r 0
- 4: p 0 := z 0
- 5: for i = 0 to T do T
- 7: x i +1 := x i + α i p i

$$6 \colon \quad \alpha _ { i } \coloneqq \frac { r _ { i } ^ { T } z _ { i } } { P _ { i } ^ { T } \Sigma _ { P _ { i } } } \\ \pi, \quad \nu \cdots \quad \beta _ { i } \coloneqq \nu \cdots + \lambda \end{array}.$$

- 8: r i +1 := r i -α i Σp i
- 9: if ∥ r i +1 ∥ &lt; ϵ then
- 10: return x ∗ := x i +1
- 11: end if

$$1 2 \colon \ z _ { i + 1 } \coloneqq \mathbf P ^ { - 1 } \mathbf r _ { i + 1 }$$

$$\begin{smallmatrix} \dashrightarrow & & - \psi _ { i } \psi _ { i } \\ 1 3 \colon & & \beta _ { i } \coloneqq \frac { \mathbf r _ { i + 1 } ^ { T } \mathbf z _ { i + 1 } } { \mathbf r _ { i } ^ { T } \mathbf z _ { i } } \\ \vdots & & \vdots & \vdots & \vdots & \vdots & \end{smallmatrix}$$

$$1 4 \colon \ \mathfrak { p } _ { i + 1 } \coloneqq \overset { \cdot } { z } _ { i + 1 } + \beta _ { i } \mathfrak { p } _ { i }$$

- 15: end for

Conjugate gradient The only remaining computationally demanding step in eq. (20) is the computation of Σ -1 U v with v = K U , X ( y -f X -ϵ ) , requiring a time complexity of O ( n 3 sg ) where n sg = |U ( η, d ) | is the size of the sparse grid U ( η, d ) . To reduce this computational intensity, we consider using the conjugate gradient (CG) method (Hestenes et al., 1952; Golub and Van Loan, 2013), an iterative algorithm for solving linear systems efficiently via matrix-vector multiplications. Preconditioning (Trefethen and Bau III, 1997; Demmel, 1997; Saad, 2003; Van der Vorst, 2003) is a well-known tool for accelerating the convergence of the CG method, which introduces a symmetric positive-definite matrix P U ≈ Σ U such that P -1 U Σ U has a smaller condition number than Σ U . The entire scheme of the preconditioned conjugate gradient (PCG) takes the form in Algorithm 1.

Preconditioner choice The choice of the precondition P U has a trade-off between the smallest time complexity of applying P -1 U operation and the optimal condition number of P -1 U Σ U . Cutajar et al. (2016) first investigated several groups of preconditioners for radial basis function kernel, then Gardner et al. (2018) proposed a specialized precondition for general kernels based on pivoted Cholesky decomposition. Wenger et al. (2022a) analyzed and summarized convergence rates for arbitrary multi-dimensional kernels and multiple preconditioners. Due to the hierarchical structure of the sparse grid U , instead of using common classes of preconditioners, we consider two-level additive Schwarz (TAS) preconidtioner (Toselli and Widlund, 2004) for the system matrix Σ U , which is formed as an additive Schwarz (AS) preconditioner (Dryja and Widlund, 1994; Cai and Sarkis, 1999) coupled with an additive coarse space correction. Since the sparse grid U ( η, d ) defined in eq. (10) is covered by a set of overlapping domains {U # » t } # » t ∈ G ( η ) where G ( η ) = { # » t ∈ N d | ∑ d j =1 t j = η } , the TAS precondition for the system matrix Σ U is then defined as

$$P _ { T A S } ^ { - 1 } \coloneqq \sum _ { \vec { t } \, \in G ( \eta ) } S _ { \vec { t } } ^ { \top } P _ { \vec { t } } ^ { - 1 } S _ { \vec { t } } + S _ { c } ^ { \top } P _ { c } ^ { - 1 } S _ { c }.$$

On the right hand side of eq. (23), the first term is the one-level AS preconditioner, and the second term is the additive coarse space correction. S # » t ∈ R n # » t × n sg is the selection matrix that projects U # » t to U , P # » t = S # » t Σ S U ⊤ # » t ∈ R n # » t × n # » t is a symmetric positive definite sub-matrix defined on U # » t . S c ∈ R n c × n sg is the selection matrix whose columns spanning the coarse space, and P c = S Σ S c U ⊤ c ∈ R n c × n c is the coarse space matrix. n c is the dimension of the coarse space, n sg := |U ( η, d ) | is the size of the sparse grid U ( η, d ) , n # » t := |U # » t | = ∏ d j =1 t j is the size of the full grid design U # » t = × d j =1 U j,t j , # » t = ( t 1 , . . . , t d ) . Figure 2 illustrates an example of the selection matrices { S # » t } # » t ∈ G ( η ) for the sparse grid U ( η = 3 , d = 2) .

Some recent work (Toselli and Widlund, 2004; Spillane et al., 2014; Dolean et al., 2015; Al Daas et al., 2023) provide robust convergence estimates of TAS preconditioner, which states that the condition number of P -1 TAS Σ U is bounded by the number of distinct colors required to color the graph of Σ U so that span { S # » t } # » t ∈ G ( η ) of the same color are mutually Σ U -orthogonal, the maximum number of subdomains that share an unknown, and a user defined tolerance τ . Specifically, given a threshold τ &gt; 0 , a GenEO coarse space (see Section 3.2 in (Al Daas et al., 2023)) can be constructed, and the following inequality holds for the preconditioned matrix P -1 TAS Σ U :

$$\kappa ( \mathbf P _ { \text{TAS} } ^ { - 1 } \Sigma _ { \mathcal { U } } ) \leq ( k _ { c } + 1 ) \left ( 2 + ( 2 k _ { c } + 1 ) \frac { k _ { m } } { \tau } \right ),$$

where k c is the minimum number of distinct colors so that each two neighboring subdomains among {U # » t } # » t ∈ G ( η ) have different colors, k m is the coloring constant which is the maximum number of overlapping subdomains

Figure 2: Selection matrices with respect to the sparse grid U ( η = 3 , d = 2) .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-572865b5-2753-48a0-b53c-405541d4293d.jpg)

Figure 3: Residuals of PCG for Σ -1 U v for tolerance τ = 10 -3 with different preconditioners where Σ U = [ K U U , + σ -2 ϵ K U , X K X , U ] , U = U ( η, d ) is the sparse grid defined in eq. (10), vector v is randomly generated. x -axis is the iteration step, y -axis is the logarithm of the residuals norm. We denote PCG with TAS preconditioners by red dots, PCG with AS preconditioners by orange crosses, PCG with Jacobi preconditioners by green triangles, CG without preconditioning by blue diamonds. Left: The sparse grid U ( η = 12 , d = 2) . Right: The sparse grid U ( η = 10 , d = 4) .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-92150577-eee2-482a-b200-fa7e5b0826cc.jpg)

sharing a row of Σ U . Since each point in U ( η, d ) belongs to at most | G ( η ) | of the subdomains, k c = | G ( η ) | . Coloring constant k m can be viewed as the number of colors needed to color each subdomain in such a way that any two subdomains with the same color are orthogonal, therefore k m = | G ( η ) | .

However, the conventional approach for constructing a coarse space to bound the smallest eigenvalue of the P -1 TAS Σ U relies on the algebraic local symmetric positive semidefinite (SPSD) splitting (Al Daas and Grigori, 2019), which is computationally complicated. To select a coarse space that is practical and less complicated, we propose constructing it based on the sparse grid U ( η , d c ) with a lower level η c &lt; η . In practice, we set η c = max {⌈ η/ 2 ⌉ , d } . Although the TAS preconditioners with our proposed coarse space do not theoretically guarantee robust convergence, empirical results presented in Figure 3 indicate that TAS tends to converge faster than one-level AS preconditioner and CG method without preconditioning. Figure 10 shows the cases where Σ U is ill-conditioned, resulting in the instability (blow-up) of the CG method without any preconditioning.

Complexity Note that S # » t is a sparse matrix with n # » t nonzeros and all of its nonzeros are one, hence computing P # » t requires O ( nn 2 # » t ) operations, and the computation of the first term of P -1 TAS costs O ( ∑ # » t ∈ G ( η ) nn 2 # » t ) . The computational complexity of the preconditioner P -1 PAS is then O ( ∑ # » t ∈ G ( η ) nn 2 # » t + n 2 c ) = O ( n | G ( η ) ( | η/d ) 2 d ) = O ( nn sg ( η/d ) 2 d ) where n c = |U ( η , d c ) | is the size of the sparse grid according to the former setting.

The complexity of the matrix-vector multiplication Σ x 0 in Algorithm 1 is O ( nn sg ) if Σ = Σ U . In theory, the maximum number of iterations of the PCG with error factor ϵ is bounded by ⌈ 1 2 √ κ ln ( 2 ϵ ) ⌉ (Shewchuk et al., 1994), where κ is the condition number of the preconditioned matrix P -1 TAS Σ U , so the time complexity of Algorithm 1 for Σ = Σ U is O (( √ κ +( η/d ) 2 d ) n sg n ) where κ = κ ( P -1 TAS Σ U ) . Therefore the time complexity of the entire posterior sampling scheme is O ( ( η +( √ κ +( η/d ) 2 d + m n n ) ) sg ) , where ν is the smoothness parameter of the Matérn kernel, n sg is the size of inducing points on a sparse grid U ( η, d ) , n is size of observations X = { x i } n i =1 , m is the size of test points X ∗ = { x ∗ i } m i =1 , d is the dimension of the observations.

## 4 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 for dimension d = 2 and d = 4 in Section 4.1. Then we apply our method to the same real-world problems as those addressed in (Wilson et al., 2020), described in Section 4.2. For prior samplings, we use the random Fourier features (RFF) with 2 6 = 64 features (Rahimi and Recht, 2007) and the Cholesky decomposition method as benchmarks. For posterior samplings, we compare our approach with the decoupled algorithm (Wilson et al., 2020) using RFF priors and the exact Matheron's update, as well as the Cholesky decomposition method. We use Matérn kernels as defined in eq. (6) with the variance σ 2 = 1 , lengthscale ω = √ 2 ν and smoothness ν = 3 2 / . Our experiments employ 'separable' Matérn kernels as specified in eq. (5), with the same base kernels and the same lengthscale parameters ω = √ 2 ν in each dimension. For the proposed method, the inducing points are selected as a sparse grid U ( η, d ) as defined in eq. (10) for ( η = 5 , d = 2) , and ( η = 6 , d = 4) respectively. We set the same nested sequence ∅ = U j, 0 ⊆ . . . ⊆ U j,t ⊆ U j,t +1 ⊆ . . . for each dimension j , and set U j,t j as hyperbolic cross points defined in eq. (11) for all j . We set the noise variates ϵ ∼ N ( 0 , σ 2 ϵ I ) with σ 2 ϵ = 10 -4 for all experiments. The seed value is set to 99 , and each experiment is replicated 1000 times.

## 4.1 Simulation

Prior sampling We generate prior samples Z = { z i } n s i =1 over the domain [0 , 1] d using Matérn 3/2 kernels, with n s = 2 6 , . . . , 2 13 respectively. Left plots in Figure 4 and Figure 11 show the time required by different algorithms to sample from the different number of points n s . We can observe that the proposed algorithm (InSG) is comparable in efficiency to the RFF method. To evaluate accuracy, we compute the 2-Wasserstein distances between empirical and true distributions. The 2-Wasserstein distance measures the similarity of two distributions. Let f 1 ∼ N ( µ , 1 Σ ) 1 and f 2 ∼ N ( µ , 2 Σ ) 2 , the 2-Wasserstein distance between the Gaussian distributions f 1 , f 2 on L 2 norm is given by (Dowson and Landau, 1982; Gelbrich, 1990; Mallasto and Feragen, 2017)

$$W _ { 2 } ( f _ { 1 }, f _ { 2 } ) \coloneqq & \left ( \left | \| \mu _ { 1 } - \mu _ { 2 } \| ^ { 2 } + \text{tr} \left ( \Sigma _ { 1 } + \Sigma _ { 2 } - 2 ( \Sigma _ { 1 } ^ { \frac { 1 } { 2 } } \Sigma _ { 2 } \Sigma _ { 1 } ^ { \frac { 1 } { 2 } } ) ^ { \frac { 1 } { 2 } } \right ) \right ) ^ { \frac { 1 } { 2 } }.$$

For the empirical distributions, the parameters are estimated from the replicated samples. The right plots in Figure 4 and Figure 11 demonstrate that the proposed algorithm performs with nearly the same precision as Cholesky decomposition and RFF.

Posterior sampling We use the Griewank function (Griewank, 1981) as the test function, defined as:

$$f ( { \mathbf x } ) = \sum _ { j = 1 } ^ { d } \frac { x _ { j } ^ { 2 } } { 4 0 0 0 } + \prod _ { j = 1 } ^ { d } \cos ( \frac { x _ { j } } { \sqrt { j } } ) + 1, \ \mathbf x \in [ - 5, 5 ] ^ { d }.$$

The sparse grid U ( η, d ) is set over [ -5 5] , d with configurations ( η = 5 , d = 2) and ( η = 6 , d = 4) for the design of our experiment. We then evaluate the average computational time and 2-Wasserstein distance over m = 1000 random test points for each sampling method. Figure 5 and Figure 12 illustrate the performance of different sampling strategies. It's observed that our proposed algorithm is significantly more time-efficient compared to the Cholesky and decoupled algorithms, particularly when n is larger than 2 10 . Moreover, our proposed method, the inducing points approximation on the sparse grid (InSG) demonstrates comparable accuracy to the decoupled method.

## 4.2 Application

Thompson sampling Thompson Sampling (TS) (Thompson, 1933) is a classical strategy for decisionmaking by selecting actions x ∈ X that minimize a black-box function f : X → R . At each iteration t , TS determines x t +1 ∈ arg min x ∈X ( f | y )( x ) , where y represents the observation set at current iteration. Upon finding the minimizer, the corresponding y t +1 is obtained by evaluating f at x t +1 , and then the pair ( x t +1 , y t +1 )

Figure 4: Time and accuracy of different algorithms for sampling from GP priors with Matérn 3/2 in dimension d = 2 . x -axis is the number of the sampling points n s . The sparse grid level of the proposed method is set as η = 5 . We denote the proposed approach, the inducing points on the sparse grid (InSG), by red dots, random Fourier features (RFF) by green triangles, Cholesky decomposition (Chol) by blue diamonds. Left: Logarithm of time taken to generate a draw from GP priors. Right: Logarithm of 2-Wasserstein distances between priors and true distributions.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-656836e2-0499-49fa-9fd2-c2062ce7459b.jpg)

Figure 5: Time and accuracy of different algorithms for sampling from GP posteriors with Matérn 3/2 and dimension d = 2 . x -axis is the number of observations n . The sparse grid level of the proposed method is set as η = 5 . The decoupled method is configured to use RFF priors with 2 8 = 64 features and exact Matheron's update. Left: Logarithm of time taken to generate a draw from GP posteriors over m = 1000 points. Right: Logarithm of 2-Wasserstein distances between posterior samplings and ground truth GP posteriors over m = 1000 points.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-21c81f78-e227-4a23-9d62-e5ea2385946a.jpg)

is added to the training set. In this experiment, we consider the Ackley function (Ackley, 2012)

$$f ( x ) = - a \exp \left ( - b \sqrt { \frac { 1 } { d } \sum _ { j = 1 } ^ { d } x _ { j } ^ { 2 } } \right ) + a + \exp ( 1 ) - \exp \left ( \frac { 1 } { d } \sum _ { j = 1 } ^ { d } \cos ( c x _ { j } ) \right ),$$

with a = 20 , b = 0 2 . , c = 2 π . The goal of this experiment is to find the global minimizer of the target function f . We start with 3 initial samples, then at each iteration of TS, we draw a posterior sample f | y on 1024 uniformly distributed points over the interval [ -5 5] , d conditioned on the current observation set. Next, we pick the smallest posterior sample for this iteration, add it to the training set, and repeat the above process. This iterative approach helps us progressively approximate the global minimum. To neutralize the impact of the lengthscales and seeds, we average the regret over 8 different lengthscales [0 2 . , 0 4 . , 0 6 . , 0 8 . , 1 0 . , 1 2 . , 1 4 . , 1 6] . and 5 different seeds [9 , 99 999 9999 99999] , , , . In Figure 6, we compare the logarithm of regret at each iteration for different sampling algorithms, the proposed approach (InSG) demonstrates comparable performance to the decoupled method.

Simulating dynamical systems GP posteriors are also useful in dynamic systems where data is limited. An illustrative example is the FitzHugh-Nagumo model (FitzHugh, 1961; Nagumo et al., 1962), which intricately

Figure 6: Regret of Thompson sampling methods of different algorithms for optimizing Ackley function with d = 2 ( left ), d = 4 ( middle ) and d = 6 ( right ) for Matérn 3/2. The sparse grid design is set as U ( η = 5 , d = 2) , U ( η = 6 , d = 4) , and U ( η = 8 , d = 6) respectively. The decoupled method is configured to use RFF priors with 2 8 = 64 features and exact Matheron's update. x -axis is the iteration step, y -axis is the regret at each step.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-af666e57-ee42-4709-b06a-c4378fe201a0.jpg)

models the activation and deactivation dynamics of a spiking neuron. The system's equations are given by:

$$\dot { \mathbf x } = f ( \mathbf x, a ) = \begin{bmatrix} \dot { v } \\ \dot { w } \end{bmatrix} = \begin{bmatrix} v - \frac { v ^ { 3 } } { 3 } - w + a \\ \frac { 1 } { \gamma } ( v - \beta w + \alpha ) \end{bmatrix}.$$

Using the Euler-Maruyama method (Maruyama, 1955), the system's equations eq. (27) can be discretized into a stochastic difference equation as follows:

$$\cdot \\ y _ { t } = x _ { t + 1 } - x _ { t } = \tau f ( x _ { t }, a _ { t } ) + \sqrt { \tau } \epsilon _ { t }, \ \epsilon _ { t } \sim \mathcal { N } ( 0, \sigma _ { \epsilon } ^ { 2 } \mathbf I ),$$

where τ is the fixed step size. Write a FitzHugh-Naguomo neuron eq. (27) in the form eq. (28), we can obtain:

$$x _ { t + 1 } - x _ { t } = \tau \begin{bmatrix} v _ { t } - \frac { v _ { t } ^ { 3 } } { 3 } - w _ { t } + a _ { t } \\ \frac { 1 } { \gamma } ( v _ { t } - \beta w _ { t } + \alpha ) \end{bmatrix} + \sqrt { \tau } \epsilon _ { t }.$$

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-91e07b0c-0618-4feb-8364-92d6e3eecb84.jpg)

Figure 7: Trajectories of different algorithms for simulation of a stochastic FitzHugh-Nagumo neuron. Each iteration is generated over 1000 samples. The sparse grid level of the proposed method (InSG) is set as η = 5 . The decoupled method is configured to use RFF priors with 2 8 = 64 features and exact Matheron's update. Top: Voltage v t trajectories generated by different algorithms. Bottom: Recovery variable w t trajectories generated by different algorithms.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-6c34b39f-8f32-4e94-ac6c-20b5a8025f1b.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-11561a2e-83d0-41cb-aa1c-e24aa04a0bda.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-60316d9b-7093-438b-b1fb-449fcd128c73.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-b842aea4-9c31-4f36-979e-064c4444af75.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-a00bb7fd-00e5-49c8-97a7-1687541a4b5c.jpg)

Our objective is to simulate the state trajectories of this dynamical system. Starting with an initial point x 0 = ( v , w 0 0 ) T , at iteration t , we draw a posterior GP sampling from the conditional distribution p ( y t | D t -1 ) , where D t -1 denotes the set of the training data { ( x i , y i ) } n i =1 along with the current trajectory { ( x j , y j ) } t -1 j =1 . In our implementation, we choose the model eq. (29) with step size τ = 0 25 . , parameters α = 0 75 . , β = 0 75 . ,

Figure 8: Time and accuracy of different algorithms for simulation of a stochastic FitzHugh-Nagumo neuron. Left: Time cost at each iteration. Right: Logarithm of 2-Wasserstein distance between simulations and ground truth state trajectories at each iteration.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-4e76c819-1d67-4a07-9680-562f2750c8b5.jpg)

γ = 20 , noise variance σ 2 ϵ = 10 -4 , and initial point x 0 = ( -2 5 . , -1) ⊤ . The training set { ( x i , y i ) } n i =1 was generated by evaluating eq. (29) at n = 256 points { x i } n i =1 ( x i ∈ X ) , which are uniformly distributed in the interval X = [ -2 5 . , 2 5] . × -[ 1 2] , , with corresponding currents { a i } n i =1 ( a i ∈ A ) chosen uniformly from A = [0 1] , . Variations y t in each step were simulated using independent Matérn 3/2 GPs with ω = √ 3 and inducing points approximation on the sparse grid U ( η = 5 , d = 2) . We generate 1000 samples at each iteration. Figure 7 presents the state trajectories of voltage v t and recovery variable w t versus the iteration step t for each algorithm. The left plots in Figure 7 show that our proposed algorithm can accurately characterize the state trajectories of this dynamical system. Figure 8 illustrates the computational time required for each iteration and the 2-Wasserstein distance between the states x t derived from GP-based simulations and those obtained using the Euler-Maruyama method as described in eq. (28) at each iteration t . These results demonstrate that our proposed algorithm not only reduces the computational time but also preserves accuracy that is comparable to the decoupled method.

## 5 Conclusion

In this work, we propose a scalable algorithm for GP sampling using inducing points on a sparse grid, which only requires a time complexity of O ( ( η + n s ) n sg ) for n s sampling points, a sparse grid of size n sg and level η . Additionally, we show that the convergence rate for GP sampling with a sparse grid is O ( n -ν sg (log n sg ) ( ν +2)( d -1)+ +1 d ) for a d -dimensional sparse grid of size n sg and a product Matérn kernel of smoothness ν . For the posterior samplings, we develop a two-level additive Schwarz preconditioner for the matrix over the sparse grid, which empirically shows rapid convergence. In the numerical study, we demonstrate that the proposed algorithm is not only efficient but also maintains accuracy comparable to other approximation methods.

## Acknowledgments and Disclosure of Funding

This work is supported by NSF grants DMS-2312173 and CNS-2328395.

## A Theorems

## A.1 Rate of convergence

Theorem A.1. Let M = L q ([0 , 1]) for any q ∈ [1 , + ∞ ) . Let Z be a GP on [0 , 1] with a Matérn kernel defined in eq. (6) with smoothness ν , and ˆ Z be an inducing points approximation of Z under the inducing points { 0 1 , /n, 2 /n, . . . , 1 } for some positive integer n . Then for any p ∈ [1 , ∞ ) , as n → + ∞ , the order of magnitude of the approximation error is

$$W _ { p } ( Z, \hat { Z } ) = \mathcal { O } ( n ^ { - \nu } ).$$

Proof. Note that Z and ˆ Z already live in the same probability space, that is ( Z, Z ˆ ) ∈ Γ( Z, Z ˆ ) . Therefore,

$$W _ { p } ( Z, \hat { Z } ) \leq \left ( \mathbb { E } \| Z - \hat { Z } \| _ { L _ { q } ( [ 0, 1 ] ) } ^ { p } \right ) ^ { 1 / p }.$$

According to Theorem 4 in (Tuo and Wang, 2020), ∥ Z -ˆ Z ∥ L q ([0 1]) , has the order of magnitude O n ( -ν ) and sub-Gaussian tails, which leads to the desired result.

Lemma A.2. [Adapted from Proposition 1 and Corollary 2 in (Rieger and Wendland, 2017)] Let Φ : R d → R be a reproducing kernel of W r ; ⊗ d 2 ( R d ) , where Φ is the tensor product with Φ( x ) := ∏ d j =1 ϕ x ( j ) and x = ( x , . . . , x 1 d ) ⊤ , ϕ : R → R is a reproducing kernel of Sobolev space W r 2 ( R 1 ) of smoothness r &gt; 1 / 2 measured in L 2 ( R 1 ) -norm. Suppose we are given the sparse grid U ( η, d ) defined in eq. (10) over I d = [0 1] , d with n = card( U ( η, d )) points and an unknown function f ∈ W r ; ⊗ d 2 ( I d ) . Let s 0 be the solution to the norm-minimal interpolant on the sparse grid U ( η, d ) as follows:

$$\min \| s \| _ { w _ { 2 } ^ { r ; \otimes ^ { d } } } \, \colon s \in W _ { 2 } ^ { r ; \otimes ^ { d } } ( \mathbb { R } ^ { d } ) \quad \text{with } \, s ( \xi ) = f ( \xi ), \quad \xi \in \mathcal { U } ( \eta, d ).$$

Then, we have

$$\| f - s _ { 0 } \| _ { L _ { \infty } ( I ^ { d } ) } \leq C n ^ { - r + 1 / 2 } ( \log n ) ^ { ( r + 3 / 2 ) ( d - 1 ) + d + 1 } \| f \| _ { W _ { 2 } ^ { r ; \otimes ^ { d } ( \mathbb { R } ^ { d } ) } },$$

where C &gt; 0 is a constant.

Proof. See proofs of Proposition 1 and Corollary 2 in (Rieger and Wendland, 2017).

Theorem A.3. Let M = L q ([0 , 1] d ) for any q ∈ [1 , + ∞ ) , d &gt; 1 is the dimension of the space M . Let Z be a GP on [0 , 1] d with a product Matérn kernel K ( x x , ′ ) = ∏ d j =1 K x 0 ( ( j ) , x ′ ( j ) ) defined in eq. (5) with variance σ 2 = 1 and base kernel K 0 ( · , · ) of smoothness ν , and let ˆ Z be an inducing points approximation of Z under the inducing points on a sparse grid U ( η, d ) defined in eq. (10) . We denote the cardinality of the sparse grid U ( η, d ) by n . Then, as n → + ∞ , the order of magnitude of the approximation error is

$$W _ { 2 } ( Z, \hat { Z } ) = \mathcal { O } ( n ^ { - \nu } ( \log n ) ^ { ( \nu + 2 ) ( d - 1 ) + d + 1 } ).$$

Proof. If ϕ in Lemma A.2 is a Matérn kernel in one dimension, we can infer that r = ν +1 2 / by Lemma 15 in (Tuo and Wang, 2020), which states that the reproducing kernel Hilbert space (RKHS) N K ( R 1 ) is a Sobolev space W ν +1 2 / 2 ( R 1 ) with equivalent norms. We denote RKHS associated with product Matérn kernel K , ( · · ) by N ⊗ d K ( R d ) and its norm by ∥ · ∥ K .

Let ˆ f be the interpolation of the function f ∈ N ⊗ d K ( R d ) with respect to the approximation ˆ Z . Clearly, N ⊗ d K ( R d ) = W ν +1 2; / ⊗ d 2 ( R d ) . Proposition 2 in (Rieger and Wendland, 2017) asserts that ˆ = f s 0 where s 0 is norm-minimal interpolant defined in Lemma A.2. Therefore, for any f ∈ N ⊗ d K ( R d ) , we have

$$\| f - \hat { f } \| _ { L _ { \infty } ( [ 0, 1 ] ^ { d } ) } \leq C n ^ { - \nu } ( \log n ) ^ { ( \nu + 2 ) ( d - 1 ) + d + 1 } \| f \| _ { W _ { 2 } ^ { \nu + 1 / 2 ; \otimes ^ { d } ( \mathbb { R } ^ { d } ) } }.$$

Additionally, Lemma 15 in (Tuo and Wang, 2020) implies that

$$\mathbb { E } \| Z - \hat { Z } \| _ { L _ { q } [ 0, 1 ] ^ { d } } ^ { 2 } = \sup _ { \| f \| _ { K } \leq 1 } | f - \hat { f } | ^ { 2 }.$$

Note that Z and ˆ Z already live in the same probability space, that is ( Z, Z ˆ ) ∈ Γ( Z, Z ˆ ) . Therefore, by combining with eq. (33) and eq. (34), we can obtain the following results:

$$\mathfrak { s } \text{ and eq. (34), we can obtain the following results: } \\ W _ { 2 } ( Z, \hat { Z } ) & \leq \left ( \mathbb { E } \| Z - \hat { Z } \| _ { L _ { q } [ 0, 1 ] ^ { d } } ^ { 2 } \right ) ^ { 1 / 2 } \\ & = \sup _ { \| f \| _ { K } \leq 1 } | f - \hat { f } | \\ & \leq \sup _ { \| f \| _ { K } \leq 1 } C n ^ { - \nu } ( \log n ) ^ { ( \nu + 2 ) ( d - 1 ) + d + 1 } \| f \| _ { K } \\ & \leq C n ^ { - \nu } ( \log n ) ^ { ( \nu + 2 ) ( d - 1 ) + d + 1 }.$$

Thus, the order of the magnitude of the inducing points approximation on the sparse grid error is W Z,Z 2 ( ˆ ) = O ( n -ν (log n ) ( ν +2)( d -1)+ +1 d ) .

## A.2 Smolyak algorithm for GP sampling

Theorem A.4. Let K , ( · · ) be a kernel function on R d × R d , and U = U ( η, d ) a sparse grid design with level η and dimension d defined in eq. (10) . U ( η, d ) = ⋃ # » t ∈ G ( η ) U 1 ,t 1 ×U 2 ,t 2 ×··· × U d,t d where # » t = ( t 1 , . . . , t d ) , G ( η ) = { # » t ∈ N d | ∑ d j =1 t j = η } , η ≥ d ≥ 2 , η, d ∈ N + . Let f ( ) · ∼ GP ( 0 , K ( · , · )) be a GP, then sampling from f ( ) · on a sparse grid U can be done via the following formula:

$$f ( \mathcal { U } ) \leftarrow \sum _ { \vec { t } \in \mathcal { P } ( \eta, d ) } ( - 1 ) ^ { \eta - | \vec { t } | } \binom { d - 1 } { \eta - | \vec { t } | } \bigotimes _ { j = 1 } ^ { d } \text{chol} \left ( K ( \mathcal { U } _ { j, t _ { j } }, \mathcal { U } _ { j, t _ { j } } ) \right ) \mathbf f _ { n } \, _ { \vec { t } } \, \delta _ { x } ( \mathcal { U } _ { \vec { t } } ),$$

$$\delta _ { x } ( \mathcal { U } _ { \vec { t } } ) \coloneqq \begin{cases} 0 & \text{if $x\in U\,\mathcal{ }U_{\vec{t} }$} \\ 1 & \text{if $x\in U_{\vec{t} }$} \end{cases}, \quad \mathcal { U } _ { \vec { t } } \coloneqq \mathcal { U } _ { 1, t _ { 1 } } \times \mathcal { U } _ { 2, t _ { 2 } } \times \cdots \times \mathcal { U } _ { d, t _ { d } }, \quad \quad ( 3 7 )$$

where P ( η, d ) = { # » t ∈ N d | max( d, η -d +1) ≤ | # » t | ≤ η } , chol ( K ( U j,t j , U j,t j ) ) is the Cholesky decomposition of K ( U j,t j , U j,t j ) , f I n t # » ∼ N ( 0 I , n # » t ) is standard multivariate normal distributed, n # » t = ∏ d j =1 t j is the size of the full grid design U # » t = U 1 ,t 1 ×U 2 ,t 2 ×··· × U d,t d associated with # » t = ( t 1 , . . . , t d ) , δ x ( U # » t ) is a Dirac measure on the set U .

Proof. We know that sampling from f ( U ) ∼ N ( 0 , K ( U U , )) can be done via Cholesky decomposition, i.e., by computing chol ( K ( U U , ) ) f n sg where f n sg ∼ N ( 0 I , n sg ) is multivariate normal distributed on the sparse grid U of size n sg = |U| . Note that chol ( K ( U U , ) ) f n sg can be regarded as a function g evaluated on U such that g ( U ) = chol ( K ( U U , ) ) f n sg . According to eq. (12), we can obtain the following equation by applying Smolyak's algorithm P ( η, d ) to g ,

$$S m o l y a k ^ { \prime } s \, a l g o r i n \, \mathcal { P } ( \eta, d ) \, to \, \bar { g }, \\ \mathcal { P } ( \eta, d ) ( g ) & = \sum _ { \vec { t } \in \mathcal { P } ( \eta, d ) } ( - 1 ) ^ { \eta - | \vec { t } | } \binom { d - 1 } { \eta - | \vec { t } | } \left ( \bigotimes _ { j = 1 } ^ { d } \mathcal { U } _ { j, t _ { j } } \right ) ( g ) \\ & = \sum _ { \vec { t } \in \mathcal { P } ( \eta, d ) } ( - 1 ) ^ { \eta - | \vec { t } | } \binom { d - 1 } { \eta - | \vec { t } | } g ( \times _ { j = 1 } ^ { d } \mathcal { U } _ { j, t _ { j } } ) \delta _ { x } ( \mathcal { U } _ { \vec { t } } ) \\ & = \sum _ { \vec { t } \in \mathcal { P } ( \eta, d ) } ( - 1 ) ^ { \eta - | \vec { t } | } \binom { d - 1 } { \eta - | \vec { t } | } \bigotimes _ { j = 1 } ^ { d } \text{chol} \left ( K ( \mathcal { U } _ { j, t _ { j } }, \mathcal { U } _ { j, t _ { j } } ) \right ) f _ { n } \cdot \delta _ { x } ( \mathcal { U } _ { \vec { t } } ), \\ \text{where} \, \delta _ { x } ( \mathcal { U } _ { \vec { t } } ) \, \text{and} \, \mathcal { U } _ { \vec { t } } \, \text{are defined in eq. (37).  Note that $g$ is defined on $\mathcal{U}$, so the tensor product operation}$$

where δ x ( U # » t ) and U # » t are defined in eq. (37). Note that g is defined on U , so the tensor product operation ⊗ d j =1 U j,t j on g is actually evaluating g on the corresponding points, which is equivalent to applying Dirac delta measure defined on set × d j =1 U j,t j to g ( ⊗ d j =1 U j,t j ) .

## B Plots

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-2321bcf3-712c-4116-b696-728713ded882.jpg)

Figure 9: Sparse grids U ( η, d ) over [0 , 1] d of level η = 4 5 6 7 , , , and dimension d = 3 .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-80316a0c-99f1-4b51-b031-f7013d49e4f7.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-e977c8fe-15a2-45a0-bb99-931f8324536b.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-f553f3ba-4ce3-4a2f-ac93-351111e16e04.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-70a6a9e0-dbfd-404d-9506-79a8d9151838.jpg)

Figure 10: Residuals of PCG for Σ -1 U v for tolerance τ = 10 -3 with different preconditioners where Σ U = [ K U U , + σ -2 ϵ K U , X K X , U ] , U = U ( η, d ) is the sparse grid defined in eq. (10), vector v is randomly generated. Left: The sparse grid U ( η = 5 , d = 2) . Right: The sparse grid U ( η = 6 , d = 4) .

Figure 11: Time and accuracy of different algorithms for sampling from GP priors with Matérn 3/2 in dimension d = 4 . x -axis is the number of grid points n s . The sparse grid level of the proposed method is set as η = 6 . Left: Logarithm of time taken to generate a draw from GP priors. Right: Logarithm of 2-Wasserstein distances between priors and true distributions.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-72c6b79f-993b-44e3-9549-ca3db852a620.jpg)

Figure 12: Time and accuracy of different algorithms for sampling from GP posteriors with Matérn 3/2 and dimension d = 4 . x -axis is the number of observations n . The sparse grid level of the proposed method is set as η = 6 . The decoupled method is configured to use RFF priors with 2 8 = 64 features and exact Matheron's update. Left: Logarithm of time taken to generate a draw from GP posteriors over m = 1000 points. Right: Logarithm of 2-Wasserstein distances between posterior samplings and ground truth GP posteriors over m = 1000 points.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00206v1-figure-449c0732-5bf2-45fa-bd38-e0d7b0d2120b.jpg)

## References

| Ackley, D. (2012). A connectionist machine for genetic hillclimbing , volume 28. Springer Science &Business Media.                                                                                                            |
|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| Al Daas, H. and Grigori, L. (2019). A class of efficient locally constructed preconditioners based on coarse spaces. SIAM Journal on Matrix Analysis and Applications , 40(1):66-91.                                          |
| Al Daas, H., Jolivet, P., and Rees, T. (2023). Efficient Algebraic Two-Level Schwarz Preconditioner for Sparse Matrices. SIAM Journal on Scientific Computing , 45(3):A1199-A1213.                                            |
| Banerjee, S., Gelfand, A. E., Finley, A. O., and Sang, H. (2008). Gaussian predictive process models for large spatial data sets. Journal of the Royal Statistical Society Series B: Statistical Methodology , 70(4):825-848. |
| Bishop, C. M. et al. (1995). Neural networks for pattern recognition . Oxford university press.                                                                                                                               |
| Cai, X.-C. and Sarkis, M. (1999). A restricted additive Schwarz preconditioner for general sparse linear systems. Siam journal on scientific computing , 21(2):792-797.                                                       |
| Chen, H., Ding, L., and Tuo, R. (2022). Kernel packet: An Exact and Scalable Algorithm for Gaussian Process Regression with Matérn Correlations. Journal of Machine Learning Research , 23(127):1-32.                         |
| Chen, H. and Tuo, R. (2022). A Scalable and Exact Gaussian Process Sampler via Kernel Packets.                                                                                                                                |
| Cole, D. A., Christianson, R. B., and Gramacy, R. B. (2021). Locally induced Gaussian processes for large-scale simulation experiments. Statistics and Computing , 31:1-21.                                                   |
| Constantinou, P., Kokoszka, P., and Reimherr, M. (2017). Testing separability of space-time functional processes. Biometrika , 104(2):425-437.                                                                                |
| Cressie, N. (2015). Statistics for spatial data . John Wiley &Sons.                                                                                                                                                           |
| Curriero, F. C. (2006). On the use of non-Euclidean distance measures in geostatistics. Mathematical Geology 38:907-926.                                                                                                      |
| Cutajar, K., Osborne, M., Cunningham, J., and Filippone, M. (2016). Preconditioning kernel matrices. In International conference on machine learning , pages 2529-2538. PMLR.                                                 |
| Demmel, J. W. (1997). Applied numerical linear algebra . SIAM.                                                                                                                                                                |
| Diggle, P. J., Ribeiro, P. J., and Christensen, O. F. (2003). An introduction to model-based geostatistics. Spatial statistics and computational methods , pages 43-86.                                                       |
| Ding, L., Tuo, R., and Shahrampour, S. (2021). A sparse expansion for deep gaussian processes. arXiv preprint arXiv:2112.05888 .                                                                                              |
| Dolean, V., Jolivet, P., and Nataf, F. (2015). An introduction to domain decomposition methods: algorithms, theory, and parallel implementation . SIAM.                                                                       |
| Dowson, D. and Landau, B. (1982). The Fréchet distance between multivariate normal distributions. Journal of multivariate analysis , 12(3):450-455.                                                                           |
| Dryja, M. and Widlund, O. B. (1994). Domain decomposition algorithms with small overlap. SIAM Journal on Scientific Computing , 15(3):604-620.                                                                                |
| Engel, Y., Mannor, S., and Meir, R. (2005). Reinforcement learning with Gaussian processes. In Proceedings of the 22nd international conference on Machine learning , pages 201-208.                                          |
| FitzHugh, R. (1961). Impulses and physiological states in theoretical models of nerve membrane. Biophysical journal , 1(6):445-466.                                                                                           |
| Frazier, P. I. (2018a). A tutorial on Bayesian optimization. arXiv preprint arXiv:1807.02811 .                                                                                                                                |
| Frazier, P. I. (2018b). Bayesian optimization. In Recent advances in optimization and modeling of contemporary problems , pages 255-278. Informs.                                                                             |
| Garcke, J. (2013). Sparse grids in a nutshell. In Sparse grids and applications , pages 57-80. Springer.                                                                                                                      |
| Gardner, J., Pleiss, G., Weinberger, K. Q., Bindel, D., and Wilson, A. G. (2018). Gpytorch: Blackbox matrix- matrix gaussian process inference with gpu acceleration. Advances in neural information processing systems       |

| Gelbrich, M. (1990). On a formula for the L2 Wasserstein metric between measures on Euclidean and Hilbert spaces. Mathematische Nachrichten , 147(1):185-203.                                                                               |
|---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| Genton, M. G. (2001). Classes of kernels for machine learning: a statistics perspective. Journal of machine learning research , 2(Dec):299-312.                                                                                             |
| Genton, M. G. (2007). Separable approximations of space-time covariance matrices. Environmetrics: The official journal of the International Environmetrics Society , 18(7):681-695.                                                         |
| Girard, A., Rasmussen, C., Candela, J. Q., and Murray-Smith, R. (2002). Gaussian process priors with uncertain inputs application to multiple-step ahead time series forecasting. Advances in neural information processing systems , 15.   |
| Gneiting, T., Genton, M. G., and Guttorp, P. (2006). Geostatistical space-time models, stationarity, separability, and full symmetry. Monographs On Statistics and Applied Probability , 107:151.                                           |
| Golub, G. H. and Van Loan, C. F. (2013). Matrix Computations . JHU press.                                                                                                                                                                   |
| Gramacy, R. B. and Apley, D. W. (2015). Local Gaussian process approximation for large computer experiments. Journal of Computational and Graphical Statistics , 24(2):561-578.                                                             |
| Grande, R., Walsh, T., and How, J. (2014). Sample efficient reinforcement learning with Gaussian processes. In International Conference on Machine Learning , pages 1332-1340. PMLR.                                                        |
| Griewank, A. O. (1981). Generalized descent for global optimization. Journal of optimization theory and applications , 34(1):11-39.                                                                                                         |
| Grigorievskiy, A., Lawrence, N., and Särkkä, S. (2017). Parallelizable sparse inverse formulation Gaussian processes (SpInGP). In 2017 IEEE 27th International Workshop on Machine Learning for Signal Processing (MLSP) , pages 1-6. IEEE. |
| Hartikainen, J. and Särkkä, S. (2010). Kalman filtering and smoothing solutions to temporal Gaussian process regression models. In 2010 IEEE international workshop on machine learning for signal processing , pages 379-384. IEEE.        |
| Henderson, H. V., Pukelsheim, F., and Searle, S. R. (1983). On the history of the Kronecker product. Linear and Multilinear Algebra , 14(2):113-120.                                                                                        |
| Hensman, J., Durrande, N., and Solin, A. (2018). Variational Fourier features for Gaussian processes. Journal of Machine Learning Research , 18(151):1-52.                                                                                  |
| Hensman, J., Fusi, N., and Lawrence, N. D. (2013). Gaussian processes for big data. arXiv preprint arXiv:1309.6835 .                                                                                                                        |
| Hensman, J., Matthews, A., and Ghahramani, Z. (2015). Scalable variational Gaussian process classification. In Artificial Intelligence and Statistics , pages 351-360. PMLR.                                                                |
| Hestenes, M. R., Stiefel, E., et al. (1952). Methods of conjugate gradients for solving linear systems. Journal of research of the National Bureau of Standards , 49(6):409-436.                                                            |
| Journel, A. G. and Huijbregts, C. J. (1976). Mining geostatistics.                                                                                                                                                                          |
| Katzfuss, M. and Guinness, J. (2021). A general framework for Vecchia approximations of Gaussian processes. Statistical Science , 36(1):124-141.                                                                                            |
| Katzfuss, M., Guinness, J., and Lawrence, E. (2022). Scaled Vecchia approximation for fast computer-model emulation. SIAM/ASA Journal on Uncertainty Quantification , 10(2):537-554.                                                        |
| Kennedy, M. C. and O'Hagan, A. (2001). Bayesian calibration of computer models. Journal of the Royal Statistical Society: Series B (Statistical Methodology) , 63(3):425-464.                                                               |
| M. and Rasmussen, C. (2003). Gaussian processes in reinforcement learning. Advances in neural                                                                                                                                               |
| Kuss, information processing systems , 16.                                                                                                                                                                                                  |
| Kuss, M., Rasmussen, C. E., and Herbrich, R. (2005). Assessing Approximate Inference for Binary Gaussian Process Classification. Journal of machine learning research , 6(10).                                                              |

| Lin, J. A., Antorán, J., Padhy, S., Janz, D., Hernández-Lobato, J. M., and Terenin, A. (2023). Sampling from Gaussian Process Posteriors using Stochastic Gradient Descent. arXiv preprint arXiv:2306.11589 .                                              |
|------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| MacKay, D. J., Mac Kay, D. J., et al. (2003). Information theory, inference and learning algorithms . Cambridge university press.                                                                                                                          |
| Maddox, W. J., Balandat, M., Wilson, A. G., and Bakshy, E. (2021). Bayesian optimization with high-dimensional outputs. Advances in Neural Information Processing Systems , 34:19274-19287.                                                                |
| Mallasto, A. and Feragen, A. (2017). Learning from uncertain curves: The 2-Wasserstein metric for Gaussian processes. Advances in Neural Information Processing Systems , 30.                                                                              |
| Maruyama, G. (1955). Continuous Markov processes and stochastic equations. Rend. Circ. Mat. Palermo , 4:48-90.                                                                                                                                             |
| Marzouk, Y. M. and Najm, H. N. (2009). Dimensionality reduction and polynomial chaos acceleration of Bayesian inference in inverse problems. Journal of Computational Physics , 228(6):1862-1902.                                                          |
| Murray-Smith, R. and Pearlmutter, B. A. (2004). Transformations of Gaussian process priors. In International Workshop on Deterministic and Statistical Methods in Machine Learning , pages 110-123. Springer.                                              |
| Nagumo, J., Arimoto, S., and Yoshizawa, S. (1962). An active pulse transmission line simulating nerve axon. Proceedings of the IRE , 50(10):2061-2070.                                                                                                     |
| Neal, R. M. (2012). Bayesian learning for neural networks , volume 118. Springer Science &Business Media.                                                                                                                                                  |
| Nguyen, V., Deisenroth, M. P., and Osborne, M. A. (2021). Gaussian Process Sampling and Optimization with Approximate Upper and Lower Bounds. arXiv preprint arXiv:2110.12087 .                                                                            |
| Nickisch, H. and Rasmussen, C. E. (2008). Approximations for binary Gaussian process classification. Journal of Machine Learning Research , 9(Oct):2035-2078.                                                                                              |
| Nickisch, H., Solin, A., and Grigorevskiy, A. (2018). State space Gaussian processes with non-Gaussian likelihood. In International Conference on Machine Learning , pages 3789-3798. PMLR.                                                                |
| Nychka, D., Bandyopadhyay, S., Hammerling, D., Lindgren, F., and Sain, S. (2015). A multiresolution Gaussian process model for the analysis of large spatial datasets. Journal of Computational and Graphical Statistics , 24(2):579-599.                  |
| O'Hagan, A. (1978). Curve fitting and optimal design for prediction. Journal of the Royal Statistical Society: Series B (Methodological) , 40(1):1-24.                                                                                                     |
| Okhrin, Y., Schmid, W., and Semeniuk, I. (2020). New approaches for monitoring image data. IEEE Transactions on Image Processing , 30:921-933.                                                                                                             |
| Pardo-Iguzquiza, E. and Chica-Olmo, M. (2008). Geostatistics with the Matern semivariogram model: A library of computer programs for inference, kriging and simulation. Computers &Geosciences , 34(9):1073-1079.                                          |
| Plumlee, M. (2014). Fast prediction of deterministic functions using sparse grid experimental designs. Journal of the American Statistical Association , 109(508):1581-1591.                                                                               |
| Plumlee, M. (2017). Bayesian calibration of inexact computer models. Journal of the American Statistical Association , 112(519):1274-1285.                                                                                                                 |
| Quinonero-Candela, J. and Rasmussen, C. E. (2005). A unifying view of sparse approximate Gaussian process regression. The Journal of Machine Learning Research , 6:1939-1959.                                                                              |
| Rahimi, A. and Recht, B. (2007). Random features for large-scale kernel machines. Advances in neural information processing systems , 20.                                                                                                                  |
| Rasmussen, C. E. (2003). Gaussian processes in machine learning. In Summer school on machine learning , pages 63-71. Springer.                                                                                                                             |
| Rieger, C. and Wendland, H. (2017). Sampling inequalities for sparse grids. Numerische Mathematik , 136:439- 466.                                                                                                                                          |
| Roberts, S., Osborne, M., Ebden, M., Reece, S., Gibson, N., and Aigrain, S. (2013). Gaussian processes for time-series modelling. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences , 371(1984):20110550. |

| Saad, Y. (2003). Iterative methods for sparse linear systems . SIAM.                                                                                                                                                               |
|------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| Saatçi, Y. (2012). Scalable inference for structured Gaussian process models . PhD thesis, Citeseer.                                                                                                                               |
| Santner, T. J., Williams, B. J., Notz, W. I., and Williams, B. J. (2003). The design and analysis of computer experiments , volume 1. Springer.                                                                                    |
| Shewchuk, J. R. et al. (1994). An introduction to the conjugate gradient method without the agonizing pain.                                                                                                                        |
| Silverman, B. W. (1985). Some aspects of the spline smoothing approach to non-parametric regression curve fitting. Journal of the Royal Statistical Society: Series B (Methodological) , 47(1):1-21.                               |
| Smolyak, S. A. (1963). Quadrature and interpolation formulas for tensor products of certain classes of functions. In Doklady Akademii Nauk , volume 148, pages 1042-1045. Russian Academy of Sciences.                             |
| Snoek, J., Larochelle, H., and Adams, R. P. (2012). Practical bayesian optimization of machine learning algorithms. Advances in neural information processing systems , 25.                                                        |
| Spillane, N., Dolean, V., Hauret, P., Nataf, F., Pechstein, C., and Scheichl, R. (2014). Abstract robust coarse spaces for systems of PDEs via generalized eigenproblems in the overlaps. Numerische Mathematik , 126:741- 770.    |
| Srinivas, N., Krause, A., Kakade, S. M., and Seeger, M. (2009). Gaussian process optimization in the bandit setting: No regret and experimental design. arXiv preprint arXiv:0912.3995 .                                           |
| Stanton, S., Maddox, W., Delbridge, I., and Wilson, A. G. (2021). Kernel interpolation for scalable online Gaussian processes. In International Conference on Artificial Intelligence and Statistics , pages 3133-3141. PMLR.      |
| Stein, M. L. (1999). Interpolation of spatial data: some theory for kriging . Springer Science &Business Media.                                                                                                                    |
| Teckentrup, A. L. (2020). Convergence of Gaussian process regression with estimated hyper-parameters and applications in Bayesian inverse problems. SIAM/ASA Journal on Uncertainty Quantification , 8(4):1310- 1337.              |
| Thompson, W. R. (1933). On the likelihood that one unknown probability exceeds another in view of the evidence of two samples. Biometrika , 25(3-4):285-294.                                                                       |
| Titsias, M. (2009). Variational learning of inducing variables in sparse Gaussian processes. In Artificial intelligence and statistics , pages 567-574. PMLR.                                                                      |
| Toselli, A. and Widlund, O. (2004). Domain decomposition methods-algorithms and theory , volume 34. Springer Science &Business Media.                                                                                              |
| Trefethen, L. N. and Bau III, D. (1997). Numerical linear algebra, vol. 50.                                                                                                                                                        |
| Tuo, R. and Wang, W. (2020). Kriging prediction with isotropic Matérn correlations: Robustness and experimen- tal designs. The Journal of Machine Learning Research , 21(1):7604-7641.                                             |
| Tuo, R. and Wu, C. J. (2015). Efficient calibration for imperfect computer models. The Annals of Statistics , 43(6):2331 - 2352.                                                                                                   |
| Ullrich, T. (2008). Smolyak's algorithm, sampling on sparse grids and Sobolev spaces of dominating mixed smoothness. East Journal of Approximations , 14(1):1.                                                                     |
| Van der Vorst, H. A. (2003). Iterative Krylov methods for large linear systems . Number 13. Cambridge University Press.                                                                                                            |
| Wang, K., Pleiss, G., Gardner, J., Tyree, S., Weinberger, K. Q., and Wilson, A. G. (2019). Exact Gaussian processes on a million data points. Advances in neural information processing systems , 32.                              |
| Wang, W., Tuo, R., and Jeff Wu, C. (2020). On prediction properties of kriging: Uniform error bounds and robustness. Journal of the American Statistical Association , 115(530):920-930.                                           |
| Wenger, J., Pleiss, G., Hennig, P., Cunningham, J., and Gardner, J. (2022a). Preconditioning for scalable Gaussian process hyperparameter optimization. In International Conference on Machine Learning , pages 23751-23780. PMLR. |
| Wenger, J., Pleiss, G., Pförtner, M., Hennig, P., and Cunningham, J. P. (2022b). Posterior and computational uncertainty in gaussian processes. Advances in Neural Information Processing Systems , 35:10876-10890.                |

| Williams, C. and Seeger, M. (2000). Using the Nyström method to speed up kernel machines. Advances in neural information processing systems , 13.                                                                                                                       |
|-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------|
| Williams, C. K. and Rasmussen, C. E. (2006). Gaussian processes for machine learning , volume 2. MIT press Cambridge, MA.                                                                                                                                               |
| Wilson, A. and Nickisch, H. (2015). Kernel interpolation for scalable structured Gaussian processes (KISS-GP). In International conference on machine learning , pages 1775-1784. PMLR.                                                                                 |
| Wilson, J., Borovitskiy, V., Terenin, A., Mostowsky, P., and Deisenroth, M. (2020). Efficiently sampling functions from Gaussian process posteriors. In International Conference on Machine Learning , pages 10292-10302. PMLR.                                         |
| Wilson, J. T., Borovitskiy, V., Terenin, A., Mostowsky, P., and Deisenroth, M. P. (2021). Pathwise Conditioning of Gaussian Processes. J. Mach. Learn. Res. , 22:105-1.                                                                                                 |
| Xie, F. and Xu, Y. (2021). Bayesian projected calibration of computer models. Journal of the American Statistical Association , 116(536):1965-1982.                                                                                                                     |
| Yadav, M., Sheldon, D. R., and Musco, C. (2022). Kernel Interpolation with Sparse Grids. Advances in Neural Information Processing Systems , 35:22883-22894.                                                                                                            |
| Zafari, S., Murashkina, M., Eerola, T., Sampo, J., Kälviäinen, H., and Haario, H. (2020). Resolving overlap- ping convex objects in silhouette images by concavity analysis and Gaussian process. Journal of Visual Communication and Image Representation , 73:102962. |
Haoyuan Chen Rui Tuo
Gaussian Processes Sampling with Sparse Grids under Additive Schwarz Preconditioner
Gaussian processes (GPs) are widely used in non-parametric Bayesian modeling,
and play an important role in various statistical and machine learning
applications. In a variety tasks of uncertainty quantification, generating
random sample paths of GPs is of interest. As GP sampling requires generating
high-dimensional Gaussian random vectors, it is computationally challenging if
a direct method, such as the Cholesky decomposition, is used. In this paper, we
propose a scalable algorithm for sampling random realizations of the prior and
posterior of GP models. The proposed algorithm leverages inducing points
approximation with sparse grids, as well as additive Schwarz preconditioners,
which reduce computational complexity, and ensure fast convergence. We
demonstrate the efficacy and accuracy of the proposed method through a series
of experiments and comparisons with other recent works.