Kernel Interpolation with Sparse Grids
Authors: Mohit Yadav
Daniel Sheldon
Cameron N Musco
Venue: NeurIPS 2022 Accept (2022)
Keywords: Gaussian processes

Abstract: Structured kernel interpolation (SKI) accelerates Gaussian processes (GP) inference by interpolating the kernel covariance function using a dense grid of inducing points, whose corresponding kernel matrix is highly structured and thus amenable to fast linear algebra. Unfortunately, SKI scales poorly in the dimension of the input points, since the dense grid size grows exponentially with the dimension. To mitigate this issue, we propose the use of sparse grids within the SKI framework. These grids enable accurate interpolation, but with a number of points growing more slowly with dimension. We contribute a novel nearly linear time matrix-vector multiplication algorithm for the sparse grid kernel matrix. We also describe how sparse grids can be combined with an efficient interpolation scheme based on simplicial complexes. With these modifications, we demonstrate that SKI can be scaled to higher dimensions while maintaining accuracy, for both synthetic and real datasets. 

Kernel Interpolation with Sparse Grids
Mohit Yadav
Daniel Sheldon
University of Massachusetts Amherst ymohit@cs ceaas.eduu
University of Massachusetts Amherst sheldon&e umass. edu
Cameron Musco
University of Massachusetts Amherst cnuscoêcs. umass edu
Abstract
Structured kernel interpolation (SKI) accelerates Gaussian process (GP) inference by interpolating the kernel covariance function using dense grid of inducing points, whose corresponding kernel matrix is highly structured and thus amenable fast linear algebra. Unfortunately, SKI scales poorly in the dimension of the input points, since the dense grid size grows exponentially with the dimension. To mitigate this issue, we propose the use of sparse grids within the SKI framework. These grids enable accurate interpolation, but with number of points growing more slowly with dimension. We contribute a novel nearly linear time matrix- vector multiplication algorithm for the sparse grid kernel matrix. We also describe how sparse grids can be combined with an efficient interpolation scheme based simplices. With these modifications, demonstrate that SKI can be scaled to higher dimensions while maintaining accuracy.
Introduction
Gaussian processes (GPs) are popular prior distributions over continuous functions for use in Bayesian inference [14]. Due to their simple mathematical structure, closed form expressions can be given for posterior inference [20]. Unfortunately, well-established limitation of GPs is that they are difficult to scale to large datasets. In particular, for both exact posterior inference and the exact lig-likelihood computation for hyperparameter learning, one must invert a dense kernel covariance matrix K "x'," where n is the number of training points. Naively, this operation requires (n))) time and (n²²) memory.
Structured Kernel Interpolation. Many techniques have been proposed mitigate this scalability issue [23, 19, 24, 9]. Recently, structured kernel interpolation (SKI) has emerged as promising approach [9]. SKI, the kernel matrix is approximated via interpolation onto a dense rectilinear grid of m inducing points. In particular, K is approximated as W K'GW, where Ke E Rm is the kernel matrix on the inducing points and W R=xm an interpolation weight matrix mapping training points to nearby grid points. Typically, W is sparse, and KG is highly structured e.g. for shift invariant kernels, KG is multi-level Toeplitz. Thus, W, KG. and in turn the approximate kernel matrix W KGW admit fast matrix-vecton multiplication. This allows fast approximate inference and log-likeli ligeliieldhood computation via the use of iterative methods, e.g., the conjugate gradient algorithm.
SKI's Curse of Dimensionality. Unfortunately. SKI does not scale well to high-dimensional input data: the number of points the dense grid, and hence the size of KG. grows exponentially in the dimension d. Moreover, SKI typically employs local cubic interpolation, which leads to an interpolation weight matrix W with row sparsity that also scales exponentially d. This curse of
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

dimensionality is well-known issue with the use of dense grid interpolation. It has been studied extensively, e.g. in the context of nigh-dimensional interpolation and numerical integration [17, 5].
In the computational mathematics community, an important technique for interpolating functions in high dimensions is sparse grids [2]. Roughly sparse grid is union of rectilinear grids with different resolutions in each dimension. In particular, it is union of all 2 x 2 x X 2 sized grids, where Di=1 / S l, for some maximum total resolution L. This upper bound on the total resolution limits the number of points each grid while grid can be dense be dense in all dimensions. See Figure for an illustration. Sparse grids have interpolation accuracy comparable to dense grids under certain smoothness assumptions on the interpolated function [22], while using significantly fewer points. Concretely, for any function f IR + I with bounded mixed partial derivatives, a sparse grid containing ( pd))) points can interpolate as accurately as a dense grid with O points, where i the maximum grid resolution [21].
few dimensions, no grid can
Combining Sparse Grids with SKI. Our main contribution to demonstrate that sparse grids be used within the SKI framework significantly improve scaling with dimension. Doing so requires several algorithmic developments. When the inducing point grid is sparse, the kernel matrix on the grid, KG. longer has simple structure. E.g.. it is not Toeplitz when the kernel shift invariant. Thus, naive matrix- vector multiplications with Kg require time that scales quadratically, rather than near-linearly in the grid size. This would significantly limit the scope of performance improvement from using sparse grid. To handle this issue, we develop near-linear' time matrix- vector multiplication (MVM) algorithm for any sparse grid kernel matrix. Our algorithm recursive, and critically leverages the fact that sparse grids can be constructed from smaller dense grids and that they have tensor product structure across dimensions. For an illustration of our algorithm's complexity versus that of naive quadratic time MVMs, see Figure (c).
24 II
(b)
Asymptotic order of growth
100
Deuse-M MVM Sparse-MVS
10
MVM
4
=
4
1
10
12
1
(a)
Figure 1: Iuustration of sparse grid construction for d 2 and maximum resolution é 2. (a) A dense grid with resolution (4,,4) has 2 and 2 points in each dimension; the grids with square points resolution I + Σ 2. (b) The sparse grid G2.2 is the union of rectilinear grids with total resolution at most 2. The 17 points with red squares belong to the sparse grid. The dense grid includes the additional points with blue dots, for total of 49 points. (c) The asymptotic order of growth in performing singlc kernel MVM operation for both grids with 2 unique points in each dimension, ignoring constants. For sparse grids, Ours-Mi -MV (i.c., the proposed MVM algorithm) improves significantly over the naive implementation (Sparse-MVM),
total
second key challenge that, while sparse grids allow for grid size that grows as a much more mild exponential function of the dimension d, the bottleneck for applying SKI on large datasets can come computing MVMs with the interpolation weight matrix E Raxm For classic high-dimensional interpolation schemes, like cubic interpolation, each row of W has ((2d) non-zero entries, i.e., the kernel covariance for each training point is approximated by a weighted sum of the covariance ( grid points. When the number of training points is large, storing W in
"Near- 'Near-linear linear' here means running in time O(m log m) for sparse grid with points.

---

memory, and multiplying by it, can become prohibitively expensive. To handle this issue, we take an approach similar to that of Kapoor al. [13] and employ simplicial basis functions for interpolation. whose support grows linearly with d. Combined with our fast MVM algorithm for sparse grid kernel matrices, this interpolation scheme lets us scale SKI to higher dimensions.
In summary, we propose the use of sparse grids to improve the scalability of kernel interpolation for GP inference relative to the number of dimensions. To this end, we develop an efficient nearly linear time matrix-vector multiplication algorithm for the sparse grid kernel matrix. Furthermore, we also propose the use of simplicial interpolation to improve scalability of SKI for both dense and sparse grids. We show empirically that these ideas allow SKI to scale to at least 10 dimensions and perform competitively with state-of-the art GP regression methods. We provide an efficient GPU implementation of the proposed algorithm compatible with GPy [9], which is available at https:// /gg thub. com/ymohit //skisg and licensed under the MIT license.
2 Background
Notation. We let N and No denote the natural numbers and N U {0} respectively. Matrices are represented by capital letters, and vectors by bold letters. I denotes the identity matrix, with dimensions apparent from context. For a matrix M, mvm(M) denotes the number of operations required to multiply M by any admissible vector.
In GP regression, the training data are modeled as noisy measurements of a random function / drawn from GP prior, denoted f ~ -N (0, -)), where k Rd x Rd - R is covariance kernel. For an input X,, the observed value modeled as y; = (x;) + 6, with E; (O,, o²). Observed training pairs (x,, yi) are collected as X [x X , Xnl 1aad and y = [y1..... Un. E 1 The kernel matrix (on training data) is Kx = [k(x, xy)li E R The GP inference tasks are to compute the posterior distribution of f given (X, which itself is Gaussian process, and to compute the marginal log-likelihoon logp(y) p(y). ). Naive approaches rely on the Cholesky decomposition of the matrix Kx = Kx + I² which takes (n²²) time; see Rasmussen [20] for more details.
To avoid the (n²2) running time of naive inference, many modern methods use iterative algorithms such the conjugate gradient (CG) algorithm to perform GP inference in a way that accesses the kernel matrix only through matrix vector multiplication (MVM), i.c., the mapping v + K x [9]. These methods support highly accurate approximate solutions to GP posterior inference task as well as hyper-parameter optimization. The complexity of posterior inference and one step of hyper- ppprmmprrretter optimization is (pm²) for CG iterations, as mvm x) n². In practice, < suffices [27, 9].
2.1
SKI: Structured Kernel Interpolation
SKI further accelerates iterative GP inference by approximating the kernel matrix in way that makes matrix-vector multiplications faster [27]. Given a set of inducing points U 12d, SKI approximates the kernel function as k(x, x') w Ky Wx', where Ku e || is the kernel matrix for the set of inducing points U, and the vector Wx [l contains interpolation weights interpolate from U any x. The SKI approximate kernel matrix Kx W W'',, where W | is the matrix with 7th row equal to Wx.-
To accelerate matrix-vector multiplications with the approximate kernel matrix, SKI places inducing points on a regular grid and uses grid-based interpolation. This leads to a sparse interpolation weight matrix W for example, with cubic interpolation there are ((4)) entries per row, so that mvm(W) ((n4))
a kernel matrix Ky that is multi-level Toeplitz (if A is stationary) [27], so that mvm (K)) O log [U]). Overall, mvm(Xx) O(n4d + I0] log Ul). which is much faster than n² for small d. However, SKI quickly becomes infcasible in higher dimensions due to the 4d entries per row of W and curse of dimensionality for the number of points in high-dimensional grid: specifically, [ = md for a grid with m points in each dimension.
2.2 Sparse Grids
Rectilinear grids. We first give a formal construction of rectilinear grids, which will later be the
foundation for sparse grids [7]. For a resolution index / E No. define the 1-d grid Q4 as the centers
3

---

of 2' equal partitions of the interval [0, 1], which gives 94 = (1/2 1 i S 2 and is odd}. The fact that the position index must be odd implies that grids for any two different resolutions are disjoint, Moreover, resolution-position index pairs (1, i) uniquely specify grid points in the union
Ule Q of 1-d rectilinear grids.
To extend rectilinear grids to d dimensions, let 1 E NJ denote a resolution vector. The corresponding rectilinear grid is given by 9 == w - where & denotes the Cartesian product. A grid point in 92 is indexed by the pair i1) of a resolution vector and position vector i, where (1,,4) gives the position in the 1-d grid S2 for dimension j. This construction of rectilinear grids yields three essential properties that will facilitate formalizing sparse grids: (1) grid On is uniquely determined by its resolution vector 1, (2) grids Ω4 and St with different resolution vectors are disjoint, (3) the size [Sa] 2||||| of a grid is determined by the L1 norm of its resolution vector.
Construction of Sparse Grids. Sparse grids use rectilinear grids as their fundamental building block [22] and exploit the fact that resolution vectors uniquely identify different rectilinear grids. Larger grids are formed as the union of rectilinear grids with different resolution vectors. Sparse grids use all rectilinear grids with resolution vector having Ln norm below specified threshold. Formally, for a resolution index é E Ng, the sparse grid Ge.d dimensions is Ge, = ULiii,, <e Sq. Figure illustrates the construction of the sparse grid G2.2 from smaller 2-d rectilinear grids with maximum resolution 2, i.e., (Ω [|||] S 2). The figure also illustrates another important fact: the sparse grid Ge,a has total of ((2)) distinct and equally spaced coordinates in each dimension, but many fewer total points than dense d-fold Cartesian product of such 1-d grids, which would have ((2dd) points.
Sparse grids have a number formal properties that are useful in algorithms and applications [22, 7]. Proposition below summarizes the most relevant ones for our work. For completeness, proof appears in appendix A. For more details, see Valentin [25].
Proposition 1 (Properties of Sparse Grid). Let Ge [0, 1] be sparse grid with any resolution E E No and dimension d E N. Then the following properties hold:
(P1) | = ( ed-))
(P2) V E N, se = ge.s S Ge
(P3) Ge,d = 1 (2, s GE-1,a-1) and Ge Ui
Property P1 shows that the size of sparse grid with ((2²) points in each dimension grows more slowly than a dense grid with the same number of points each dimension, since ( < O(2 Properties P2 and P3 are consequences of the structure of the set {|||]| { and the sparse grid construction. Property P2 says that a sparse grid with smaller resolution is contained in one with higher resolution. Property P3 is crucial property, and says that d d-dimensional sparse grid can be constructed via Cartesian products of I-dimensional dense grids with sparse grids in dimensions.
3 Structured Kernel Interpolation on Sparse Grids
To scale kernel interpolation higher dimensions, we propose to select inducing points = Ge,d a on a sparse grid and approximate the kernel matrix as W KG W for a suitable interpolation matrix W adapted to sparse grids. This will require fast matrix-vector multiplications with the sparse grid kernel matrix Kg, a and the interpolation matrix W. We show how to accomplish these two tasks in Sections 3.1 and 3.2 for the important case of stationary product kernels [10].
3.1 Fast Multiplication with the Sparse Grid Kernel Matrix
Algorithm is an algorithm compute Kg v for any vector v. The algorithm uses the following definitions. For any finite set U, Kr = [ x') x')x The rows and columns of Kü are "U-indexed", meaning the entries correspond to elements of U under some arbitrary fixed ordering.
Suppose i/2/t1 +1 are both grid points and k L. Then i j2*=1 ! is even, contradiction.

---

For U S V. we introduce selection matrix Su.v to map between U-indexed and V- V-indexed vectors. has entries (S equal one if the ith element of equal the jth element of V. and zero otherwise. Also, Sv.u S For V Viindexed vector zy, the multiplication Sc.y Zy produces U-indexed vector by selecting entries corresponding to elements in U, and for a U-indexed vector zu. the multiplication Sv,uzu zu produces V-indexed vector by inserting zeros for elements
U.
Algorithm Sparse Grid Kernel-MVM Algorithm
Input: RR9 and K R 9r.al
Output:
mvm ( v), where, mvm (K, v) denotes K'v obtained using this algorithm. 1: Let V; be the result of reshaping Sp, 12, 1 v into a [Ω4] x Ge-i,d-1 matrix; this contains entries corresponding to the ith grid in the decomposition of P3.
2: ifd then
return K
Base case
4: end if
5: & Pre-computation for to / do
7: A; = Kg Sq,, S2 V
8: B, mvm (K ia 11 V,T)
Recursively multiply Kg. 5.a with columns of V/F.
end for
10: & Main loop
11: for i to do
)
12: AT mvm
Kg Σ Sq, AjSgr- 1.S9--500-1B1--
Recurse like line 8.
13:
B, Sp, Kon.
Σ S (.
14: u, = vec [Ai] + vec [Bi]
15: end for
Theorem 1. Let Ker,a be the kernel matrix for a d-dimensional sparse grid with resolution for stationary product kernel. For any E Rl Algorithm computes Kg in ((d22)) time.
The formal analysis and proof of Theorem appears in Appendix B. The rumning time of Ogd2))) is nearly linear in [Ge,al and much faster asymptotically than the naive MVM algorithm that matcrializes the full matrix and has running time quadratic in (Gra).
Algorithm is built on two high-level observations. First, the decomposition of Property P3 from Proposition and the fact that the kernel follows product structure across dimensions are used to decompose the MVM into blocks, each of which is between sub- sub-grids which are the product of I I-dimensional rectilinear grid and sparse grid in d - -1 dimensions. Therefore, the overall MVM computation can recursively decomposed into MVMs with sparse grid kernel matrices dimensions. This observation is in part inspired by Zeiser [28], which also decomposes computation with the matrix on sparse grids by the resolution of first dimension. The base case occurs when 1. We assume that Toeplitz structure, which arises due to the kernel being stationary, is leveraged to perform this base case MVM in O(E2f) time. Algorithm I can also be extended to non-stationary product kernels by using standard MVM routine for the base case, which changes the overall running-time analysis but is still more efficient than the naive algorithm.
Secondly, by Property P2, the kernel matrix multiplication for any grid of resolution é also includes the result of multiplication for grids of lower resolution and the same number of dimensions. Thus, the results of the multiplications for many individual blocks can be obtained by using the appropriate selection operators with the result of the multiplication with the kernel matrix Kg in Line 12 of the algorithm. Further intuition and explanation provided in the Appendix B.I.

---

Improving batching efficiency. The recursions in Lines 8 and 12 can be batched for efficiency, since both are multiplications with the same symmetric kernel matrix Ke Similarly, the recursion spawns many recursive multiplications with kernel matrices of the form Ge dd' for S l < l and 1 1 d' < d, and the calculation can be reorganized to batch all multiplications with each K This is significant savings, because there are only d({ ((E+1) distinct kernel matrices, but the recursion has a branching factor of (f + 1), so spawns many recursive calls with the same kernel matrices.
3.2 Sparse Interpolation For Sparse Grids
We now seek to construct the matrix W, which interpolates function values from the sparse grid Ge,d to training points x; E Rd, while ensuring that each row of W is sparse enough to preserve efficiency of metrix-vector multiplications with W. This requires sparse interpolation rule for sparse grids.
To set up the problem, we consider interpolating function f observed at points in generic set U. Let fo = [ff(x) A linear interpolation rule for U is a mapping x + Wx E R used to approximate fxx at an arbitrary point as f(x = w fc. The density of an interpolation rule is the maximum number of non-zeros in Wx for any x.
The combination technique for sparse grids constructs an interpolation rule by combining interpolation rules for the constituent rectilinear grids.
For each I, let w' be an interpolation rule for the rectilinear grid 9 with maximum density e. The combination technique gives an interpolation rule wx. for the sparse grid Ge,a with density most e (É+d-1). The factor (E+d-1) is the number of rectilinear grids included in Ge
Proposition 2.
(d=1)
We use the combination technique to construct the sparse-grid interpolation coefficients Wx. for each training point x; and stack them in the rows of W, which can be done in time proportional to the number of nonzeros. Details are given in Appendix A. The combination technique can use any base interpolation rule for the rectilinear grids, such multilinear, cubie, simplicial interpolation.
3.3 Simplicial Interpolation
The density of the interpolation rule is critical consideration for kernel interpolation techniques with or without sparse grids. Linear and cubic interpolation in d dimensions have density O(2") and ((4*), respectively. This "second curse of dimensionality" makes computations with W intractable in higher dimensions independently of operations with the grid kernel matrix. For sparse grids, the density of W increases by additional factor of 1==
We propose to use simplicial interpolation [11] for the underlying interpolation rule avoid expo- nential growth of the density. Simplicial interpolation refers to a scheme where Rd is partitioned into simplices and a point interpolated using only the d+ 1 extreme points of the enclosing simplex, so the density of the interpolation rule exactly d d Simplicial interpolation was previously proposed for sparse grid classifiers in [8]. In work closely related to ours, Kapoor et al. [13] used simplicial
interpolation for GP kernel interpolation, with the key difference that they use the permutohedral
lattice as the underlying grid, which has a number of nice properties but does not come equipped with fast specialized routines for kernel matrix multiplication.
In contrast to Kapoor et al. [13], we maintain rectilinear and/or sparse underlying grids, which preserve structure that enables fast kernel matrix multiplication. For rectilinear grids, this requires partitioning each hyper-rectangle into simplices, so the entire space is partitioned by simplices whose extreme points belong to the rectilinear grid. Then, within each simplex, the values at the extreme points are interpolated lincarly. In general, there are different ways to partition hyper-rectangles into simplices-we use the specific scheme detailed in [11]. For sparse grids, we then use the combination technique, leading to overall density of (d + 1) ( Table provides the MVM complexities of W and kernel matrices for different interpolation schemes and both grids. More details on how to perform simplicial interpolation with rectilinear and sparse grids are given in the Appendix B.2.

---

Table 1: The MVM complexities (i.e., mvm(()) of the interpolation matrix W and the kernel matrix Ke for different interpolation bases and d-dimensional grids with 2 unique points in each dimension and n. data points. A star indicates approaches proposed in this work. 'Dense' denotes rectilinear grid.
Grid
Basis
mv)((W)
mvm (c)
Cubic
O(n 44)
Ole d. d
Dense
Linear
O(n.27
Oit. d. d
Simplicial"
O(n.d*
O(t: d
Cubic*
Onn.4 44. (
d=1
O .2)
Sparse
Linear*
O(n.2 24.
[1+d
O 2
Simplicial
O(n
[tt
O 2²)
d-1
4 Experiments
In this section, we empirically evaluate the time and memory taken by Algorithm I for matrix- multiplication with the sparse grid kernel matrix, the accuracy of sparse grid interpolation and GP regression as the data dimension d increases, and the accuracy of sparse grid kernel interpolation for GP regression on real higher-dimensionnll datasets from UCL. Hyper-parameters, data processing steps, and optimization details are given in Appendix C.I.
MVM time [in secs)
10a
Build time (in secs)
10
109
10
10
Spone- XVM Chm-AVVM- CUM-MMMM-
10-4
10
10a
10a
10
10a
10a
10
|
|9]
10
10
Memory (in GB)
CG time (in secs)
10a
10-
10"
10-3
10
10
10
10
10F
10 |
10
G
and increasing resolution f. From top-left to bottom-right: time for one MVM operation, build (pre- preepprccessig) time for the kernel matrix, peak memory usage, and typical time to solve a linear system using CG (time for build plus 50 MVMss. Each plot shows the naive quadratic MVM algorithm (Sparse-MVN) the recursive implementation of Algorithm (Ours-Ay MVM-RV R), and the efficiently batched iterative implementation Algorithm (Ours- MVM-I). All measurements are averaged over trials; enror hars represent twice standard error.
Figure 2: Matrix-vector multiplication (MVM) resource usage relative to sparse grid size
Sparse grid kernel MVM complexity. First, we evaluate the efficiency of MVM algorithms. We compare the basic and efficient implementations of Algorithm the naive algorithm, which constructs the full kernel matrix and scales quadratically with the sparse grid size. Algorithm has significant theoretical advantage terms of both time and memory requirements as the grid size grows. Figure 2 illustrates this for d = 6 by comparing MVM time and memory requirements for sparse grids with resolutions E (2..... 9} (roughly 100 1M grid points). The MVM time, preprocessing time, and memory consumption of Algorithm all grow more slowly than the naive
7

---

algorithm, and the efficient implementation of Algorithm is faster for [g] larger than about 10t, after which the naive algorithm also exceeds the 10 GB memory limit. For comparison, at é (about 40K grid points), Algorithm uses only 0.05 GB of memory. These results indicate that the proposed algorithm is crucial for enabling sparse grid kernel interpolation in higher dimensions. Figure also depicts the typical time to run GP inference (i.e., preprocessing time plus 50 MVM operations).
Sparse grid interpolation and GP inference accuracy on synthetie data. A significant advantage of sparse grids over dense rectilinear grids is their ability to perform accurate interpolation in higher dimensions. We demonstrate this for d 6 by interpolating the function (xx cos [| ( from observation locations on sparse grids of increasing resolution onto 200 random points sampled uniformly from [0. 1] Figure 3, left, shows the interpolation error for dense and sparse grids with both cubic and simplicial interpolation. Sparse (cubic) is significantly more accurate than Dense (cubic), and Sparse (simplicial) more accurate than Dense (simplicial).
We next evaluate the accuracy of GP inference in increasing dimensions. We keep the same function fx) and generate observations as y f(x + (0, 0.05). For all d, we use é 4 for the sparse grid and compare to the dense grid with the closest possible total number of grid points (i.e., [g], points in each dimension) For d 2 8, performance is better with sparse grids than dense grids for both interpolation schemes. Remarkably, our proposal to use simplicial interpolation with dense grids allows SKI to scale typically infeasible for d 2 4.
10, which is significant improvement over prior work, in which SKI is
10
Interpolation Eroo
Dani: (cubk)
0.2
Sparee (culise
10-
Test RMSE
Dense (simpáeaa))
Spari (simplicial)
10-2
Dana: (subic) (cubed) Dene (simolittc)
0.1
10-3
Spare (сmppeaal)
10
10
10
10
0.0
2
8
10
|9]
Figure 3: Comparing dense and sparse grids with cubic and simplicial interpolation schemes on synthetic data. Dense (cubic/simplicia)) dense grid methods with cubic/simplicial interpolation schemes; similarly, sparse (cubic/simplicial) are for sparse grids. Left: Function interpolation error vs. the grid size d = 6. Right: Test root-mean-squaee error (RMSE) for GP regression for increasing dimensions. For both tasks, sparse grid methods outperform dense grids. For d = 10, both methods with cubic interpolation run out of GPU memory, which is 48 GB for this experiment.
GP regression performance UCI datasets. To evaluate the effectiveness of our proposed methods for scaling GP kernel interpolation to higher dimensions, we consider all UCI [4] data sets with dimension 8 < d < 10. We compare our proposed methods, Dense-grid (dense SKI with simplicial interpolation) and Sparse-grid (sparse-grid SKI with simplicial interpolation) to SGPR [24], SKIP [10], and Simplex-GP [13]. For SGPR, report the best results using 256 or 512 inducing points. For SKIP, 100 points per dimension are used. For Simplex-GP [13], the blur stencil order set to Table 3 shows the root mean squared error (RMSE) for all methods. Our methods have performance comparable to and often better than SGPR and prior methods for scaling kernel interpolation higher dimensions. This shows that sparse grids and simplicial interpolation can effectively scale SKI to higher dimensions to give a GP regression framework that is competitive with state-of ofte--of-th-attt art approaches. We report additional results and analysis in Appendix C.
'We tried to match the sizes of both dense and sparse grids while ensuring that the dense grid always had at least as many points as the sparse grid, to give a fair comparison. Precisely, (d, dense grid size, sparse grid size) tuples are (2, 144, 129). (4, 1296, 796). (6, 4096, 2561), (8, 6561, 6401). (10, 59049, 13441).

---

Table 2: Test rool-meen--suuare-rroo (RMSE) UCI regression datasets with dimensions (i.e., 8 S S 10).
See text for algorithm descriptions and settings. All mean and standard deviations are computed over trials. Pendulum is abbreviated as Pendu to accommodate table horizontally.
Datasets (d)
SGPR
SKIP
Simplex GP
Densc-gra grid
Sparse- grid
Energy (8)
1.509 0.003
5.762 + 0.000
3.076 + 0.012
1.333 0.009
0.715 1 0.004
Energy (8)
1,509 0.003
5.762 + 0.000
3,076 + 0.012
1.333
0.009
0.715 + 0.004
Concrete (8)
12.727 = 0.018
12.727 = 0.001
12.727 + 0.000
12.191 = 0.001
8.655 + 0.002
Kin40k (8)
0.168 + 0.009
0.174 +: 0.001
0.287 + 0.003
0.205 + 0.003
0.483 + 0.000
Fertility (9)
0.197 = 0.016
0.183 + 0.000
0.187 + 0.001
0.182 = 0.000
0.194 + 0.002
Pendu (9)
1.948 1 0.021
2.947 1: 0.000
2.577 1. 0.009
2.053 1: 0.010
2.103 1 0.015
Protein (9)
0.605 0.001
0.778 + 0.000
0.582 + 0.018
0.736 = 0.002
0.595 + 0.001
Solar (10)
0.790 = 0.026
0.780 = 0.002
0.792 + 0.000
0.775 : 0.002
0.748 + 0.006
Related Works
Beyond SKI and its variants, a number scalable GP approximations have been investigated. Most notable are the different variants of sparse GP approximations [26, 23, 18]. For inducing points. these methods require either ((nm²) time for direct solves or (2nm) time for approximate kernel MVMs iterative solvers. While these methods generally do not leverage structured matrix algebra like the SKI framework, and thus have worse scaling in terms of the number of inducing points m, they may achieve comparable accuracy with smaller m, especially in higher dimensions. Some work utilizes the SKI framework to further boost the performance of sparse GP approximations [12].
closely related work to ours is on improving the dimension scaling of SKI through the use of low-rank approximation and product structure [10]. In another very closely related work, Kapoor al. [13] recently proposed scale SKI to higher dimensions via interpolation on a permutohedral lattice. Like our work, they use simplicial interpolation. Unlike our work, the kernel matrix on the permutohedral lattice does not have special structure that admits fast exact multiplications: they instead use a locality-based approximation that takes into account the length scale of the kernel function by only considering pairs of grid points within a certain distance. Another related direction of research is the adaptation of sparse grid techniques for machine learning problems, e.g., classification and regression [16] and data mining [1]. These methods construct feature representations using sparse grid points [3], often by selecting a subset of grid points relevant for the dataset [6, 3].
Discussion
This work demonstrates that two classic numerical techniques, namely, sparse grids and simplicial interpolation, can be used scale GP kernel interpolation to higher dimensions. SKI with sparse grids and simplicial interpolation has better or competitive regression accuracy compared to state-of-the- GP regression approaches several UCI benchmarking datasets with 8 to 10 dimensions.
Limitations and future work. Sparse grids and simplicial interpolation address two important bottlenecks when scaling kernel interpolation to higher dimensions. Sparse grids allow scalable matrix-vector multiplications with the grid kernel matrix, and simplicial interpolation allows scalable multiplications with the interpolation matrix W. The relatively large number of rectilinear grids used to form sparse grid-i.e., the factor of ((+d-") in Proposition 2-is one limiting factor that makes multiplication by W more costly. Future research could investigate methods to mitigate this extra cost, and explore the limits of scaling even higher dimensions with sparse grids.
Acknowledgments
This material is based upon work supported by the National Science Foundation under Grant Nos. 1661259, 1749854, 2046235 and 1763618, along with Adobe and Google Research Grants. We thank Javier Burroni, Abhinay Agarwal, Shiv Shankar and Shreyas Chaudhari for their feedback.

---

References
[1] H-J Bungartz, Dirk Pflüger, and Stefan Zimmer. Adaptive sparse grid techniques for data mining. In Modeling, simulation and optimization of complex processes, pages 121-130. Springer, 2008.
[2] Hans-Joachim Bungartz and Michael Griebel. Sparse grids. Acta Numerica, 13:147-269, 2004.
[3] Tri Dao, Christopher De Sa, and Christopher Ré. Gaussian quadrature for kernel features. In Advances in Neural Information Processing Systems 30 (NeurIPS), page 6109-6119, 2017.
[4] Dheeru Dua and Casey Graff. Uci machine learning repository, 2017. URL http:////arcc ics. uci edu/ml.
[5] A. Dutt, Miao Gu, and Vladimir Rokhlin. Fast algorithms for polynomial interpolation, integra- tion, and differentiation. SIAM Journal on Numerical Analysis, 33:1689-1711, 1996.
[6] Jochen Garcke. A dimension adaptive sparse grid combination technique for machine learning. Anziam Journal, 48:C725-C777 2006.
[7] Jochen Garcke. Sparse grids in nutshell. In Sparse grids and applications, pages 57-80. Springer, 2012.
[8] Jochen Garcke and Michael Griebel. Classification with sparse grids using simplicial basis functions. Intelligent data analywis, 6(6) 2002.
(9] Jacob Gardner, Geoff Pleiss, David Bindel, Kilian Q. Weinherger, and Andrew Gordon Wilson. GPyTorch: Blackbox matrix-matrix Gaussian process inference with GPU acceleration. In Advances in Neural Information Processing Systems 31 (NeurIPS), pages 7576-7586, 2018.
[10] Jacob Gardner, Geoff Pleiss, Ruihan Wu, Kilian Q Weinberger, and Andrew Gordon Wilson. Product kernel interpolation for scalable Gaussian processes. arXiv: 00000000000... 2018.
[11] John H Halton. Simplicial multivariable linear interpolation. Tehenical Report 91-002, Univer- sity of North Carolina at Chapel Hill Department of Computer Science, 1991.
[12] Pavel Izmailov, Alexander Novikov, and Dmitry Kropotov. Scalable gaussian processes with billions of inducing inputs via tensor train decomposition. In In Proceedings of the 21st International Conference on Artificial Intelligence and Statistics (AISTATS), volume 84 of Proceedings of Machine Learning Research, pages 726-735. PMLR, 2018.
[13] Sanyam Kapoor, Mare Finzi, Ke Alexander Wang. and Andrew Gordon Gordon Wilson. Skiing on simplices: Kernel interpolation on the permutohedral lattice for scalable gaussian processes. In Proceedings of the 38th International Conference Machine Learning (ICML., pages 5279-5289, 2021.
[14] Radford M. Neal. Bayesian Learning for Neural Networks. Springer- Verlag, 1996.
[15] Michael Obersteiner and Hans-Joachim Bungartz. A generalized spatially adaptive sparse grid combination technique with dimension-wiss refinement. SIAM Journal Scientific Computing, 43(4):A2381-A2400 2021.
[16] Dirk Michael Pilüger. Spatially adaptive sparse grids for high-dimensional problems. PhD thesis, Technische Universität München, 2010.
[17] George M Phillips and Peter J Taylor. Theory and applications of numerical analysis. Elsevier, 1996.
[18] Joaquin Quiñooner--Cnnder and Carl Edward Rasmussen. A unifying view of sparse approxi- mate Gaussian process regression. Journal of Machine Learning Research, 6(Dec): 1959, 2005.
[19] Ali Rahimi and Benjamin Recht. Random features for scage kernel machines. In Advances in Neural Information Processing Systems 20 (NeurIPS), pages 1177-1184, 2007.
[20] Carl Edward Rasmussen. Gaussian Processes in Machine Learning. Springer, 2004.
[21] Winfried Sickel and Tino Ullrich. Spline interpolation on sparse grids. Applicuble Analysis, 90 (3-4):3773833, 201
[22] Sergei Abramovich Smolyak. Quadrature and interpolation formulas for tensor products of certain classes of functions. In Doklady Akademi Nauk, volume 148, pages 1042-1045. Russian Academy of Sciences, 1963.
10

---

[23] Edward Snelson and Zoubin Ghahramani. Sparse Gaussian processes using pseudo-inputs. In Advances Neural Information Processing Systems 18 (NeurIPS), pages 1257-1264, 2005.
[24] Michalis Titsias. Variational learning of inducing variables in sparse gaussian processes. In In Proceedings of the 12th International Conference on Artificial Intelligence and Statistics (AISTATS), pages 567-574. PMLR, 2009.
[25] Julian Valentin. B-splines for sparse grids: Algorithms and application to higher- diaher-dimennionall optimization. arXiv:1910.05379, 2019.
[26] Christopher KI Williams and Matthias Seeger. Using the Nyström method to speed up kernel machines. In Advances in Neural Information Processing Systems 14 (NeurIPS), pages 682-688, 2001.
[27] Andrew Wilson and Hannes Nickisch. Kernel interpolation for scalable structured Gaussian processes (KISS-GP). In Proceedings of the 32nd International Conference on Machine Learning (ICML), pages 1775-1784, 2015.
[28] Andreas Zeiser. Fast matrix-vecter multiplication in the sparse-guu gause--rrddee galerkin method. Journal of Scientific Computing, 47(3):328-346, 2011.

---

Checklist
1. For all authors...
(a) Do the main claims made in the abstract and introduction accurately reflect the paper's
contributions and scope? [Yes]
(b) Did you descrihe the limitations of your work? [Yes] See Section 6.
(c) Did you discuss any potential negative societal impacts of your work? [N/A] As our method is an improvement on existing methods therefore it doesn't have any direct negative consquences.
(d) Have you read the ethics review guidelines and ensured that your paper conforms to them? [Yes]
2. If you are including theoretical results...
(a) Did you state the full of assumptions all theoretical results? [Yes]
(b) Did you include complete proofs of all theoretical results? [Yes]
3. If you ran experiments...
(a) Did you include the code, data, and instructions needed to reproduce the main experi- mental results (either the supplemental material or as a URL)? [Yes] See Section C.
(b) Did you specify all the training details (e.g., data splits, hyperparameters, how they were chosen)? [Yes] See Section C.I.
(c) Did you report error bars (e-g., with respect the random seed after running experi- ments multiple times)? [Yes] See Figures 2 and 3.
(d) Did you include the total amount of compute and the type of resources used (e.g., type of GPUs, internal cluster, or cloud provider)? [Yes] See C.
4. If you are using existing assets (e.g., code, data, models) or curating/releasing new assets...
(a) your work uses existing assets, did you eite the creators? [Yes]
(b) Did you mention the license of the assets? [Yes]
(c) Did you include any new assets either in the supplemental material or as a URL? [Yes] See Section 1.
(d) Did you discuss whether and how consent obtained from people whose data you're you're using/curating? [N/A] We use only publicly available datasets.
(c) Did you discuss whether the data you are using/curating contains personally identifiable information or offensive content? [N/A] We use only publicly available datasets.
5. If you used crowdsourcing or conducted research with human subjects...
(a) Did you include the full text of instructions given participants and screenshots, if applicable? [N/A]
(b) Did you describe any potential participant risks, with links to Institutional Review Board (IRB) approvals, applicable? [N/A]
(c) Did you include the estimated hourly wage paid to participants and the total amount spent participant compensation? [N/A]
12


Reviews:
[{"review_id": "mYz-R3bSGSv", "reviewer": "Reviewer_9bwU", "summary": "This manuscript proposes using sparse grids within the structured kernel interpolation framework to accelerate matrix vector products necessary when working with a kernelized Gaussian Process. Of particular note, sparse grids are preferable relative to dense grids as the dimension of the problem increases. A fast algorithm is provided for computing matrix vector products with the sparse grid based inducing point matrix and appropriate interpolation schemes are discussed. Finally, numerical results show the efficacy of the method for both synthetic and real world benchmark problems.", "questions": "- Is the plot in figure 1 (c) meant to show arithmetic operations or is this an actual timing experiment? \"Complexity\" is not a particularly descriptive y-axis. The statement \"improves\" suggests this is an actual experiment, but \"complexity\" suggests this is simply flop counting.\n\n- Why are $\\bar{\\mathbf{a}}_i$ and $\\bar{\\mathbf{b}}_i$ lower case when they represent matrices? Also, maybe it would be clearer to make the recursive calls more closely match the algorithm input/output definition (technically, they can be fit into it by just calling the algorithm several times, but it does read a bit oddly given the input statement).\n\n- Is there an explanation for the high variance in certain data points of Figure 2?\n\n- Is Figure 2 showing the total MVM time or just the part with the kernel matrix (i.e., only Algorithm 1)? the text suggests this is an evaluation of Algorithm 1 and, therefore only $K$ is considered. However, by showing CG time the manuscript suggests it is the whole process including the interpolation matrices.\n\n- In figure 3 what is the dense grid size? is it made to match the total number of points in the sparse grid or does each dimension match the maximum resolution. The results suggest the former (since the sparse grid interpolation is doing better), but this should be made clear.\n\n- Is there a good explanation for why the sparse (simplicial) interpolation seems to be improving more slowly than simplicial interpolation using the dense grid?", "limitations": "- One limitation not discussed is the need to choose $\\ell.$ While not dissimilar to parameters that appear in related methods (so I don't think this is a problem, per se), the manuscript does not provide much guidance on how select this parameter based on, e.g., the desired accuracy and that is a current limitation.", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Overall, I think that this paper makes a good contribution to the SKI framework. While sparse grids are certainly a well-studied and developed area, to the best of my knowledge they have not found use in this setting. Moreover, the development of a fast algorithm for applying the resulting inducing point matrix is key for this application. Accordingly, the idea + algorithm development represents a clear strength of the manuscript. Moreover, the exposition is relatively easy to read and the key points are well explained.\n\nThe main weakness of the paper is in its experimental evaluation of the method. Some of this may be relatively easy to address by clarifying certain details (see some specific questions below). Nevertheless, there are a few larger issues outlined below: \n\n- The main focus of the paper seems to be on the efficiency of all parts of the SKI pipeline (i.e., both the inducing point matrix vector product and the interpolation). However, the most \"comparative\" experiment (table 1) simply looks at RMSE and comes off a bit as: \"all these methods do reasonable in most cases and no one dominates the rest.\" While this is fine for showing the final result of the proposed method is comparable with other work, it does little to highlight the advantages (and/or potential disadvantages) of the proposed method relative to prior work. While there is a very small comparison in the appendix, given the focus of the paper it seems that balanced the numerical experiences more towards exploring the computational efficiency of the method would be preferable. For example, is the proposed method uniformly preferable? or when in say $d=2$ or $d=3$ does the practical efficacy of the FFT help dense grid methods win out? Where might one see the crossover to sparse grids being preferable? and more. Ideally the numerical experiments could help practitioners understand which method should be used in which settings (at least to a first approximation). \n\n- The manuscript would benefit from a more clear breakdown between the interpolation time and the inducing point matrix multiplication. For example, which one is the dominant cost? This would help direct future work and also make it easier to evaluate the limitations where the manuscript suggests improvements to the interpolation are important (suggesting that this ends up being the dominant cost).\n\nIn summary, there is a clear algorithmic contribution in the manuscript that suggests the possibility of a practically efficient method in slightly larger dimensions than SKI is typically used (though there are still limits). Moreover, it seems likely that these idealized gains translate to practical improvements over prior methods. However, the experimental evaluation makes assessing the latter point somewhat challenging and the manuscript would benefit significantly by making a clearer case here.\n\nAssorted notes:\n\n- The definition of $V_i$ could be made a bit more clear; maybe a figure or more precise mathematical statement in the appendix?\n\n- The fact that certain base case operations in Algorithm 1 are Topelitz and, therefore, can be accelerated using the FFT should probably be mentioned in the main text discussion of Algorithm 1. It is clear in the appendix, but important enough to warrant inclusion above.\n\n- The authors note that SKI is typically infeasible for $d \\geq 4$ and then immediately proceed to run a variant of SKI in higher dimensions. I think it is better to stick to the more general \"scaling\" arguments rather than trying to draw somewhat arbitrary boundaries between methods based on the dimension.\n\n- Proposition 2 is styled different than the other theoretical statements (e.g., Proposition 1 and Theorem 1). While either is fine, the styling should be consistent.\n\n- The related work section focuses mainly on SKI type methods. While that being the focus is certainly appropriate given the contribution of this manuscript, given that the problem being solved is GPs and not \"improve SKI\" the manuscript may want to allude to the fact that there are other \"styles\" of algorithms that are also being developed to try and scale to larger dimensions (though they all certainly have their own limitations and pros/cons relative to SKI). For example, one common strategy is the use of rank-structured factorizations (as with SKI these are most effective for $d\\leq 3$ but there has been work on \"scaling up\" to more moderate dimensions:\n\n[Yang, Changjiang, Ramani Duraiswami, Nail A. Gumerov, and Larry Davis. \"Improved fast gauss transform and efficient kernel density estimation.\" In Computer Vision, IEEE International Conference on, vol. 2, pp. 464-464. IEEE Computer Society, 2003.]\n\n[Ambikasaran, Sivaram, Daniel Foreman-Mackey, Leslie Greengard, David W. Hogg, and Michael O’Neil. \"Fast direct methods for Gaussian processes.\" IEEE transactions on pattern analysis and machine intelligence 38, no. 2 (2015): 252-265.]\n\n[Ryan, John P., Sebastian E. Ament, Carla P. Gomes, and Anil Damle. \"The fast kernel transform.\" In International Conference on Artificial Intelligence and Statistics, pp. 11669-11690. PMLR, 2022.]\n\n- It seems ref 24 is a duplicate of 23? I believe this is just one paper from 2005 unless there is another reference I am unfamiliar with that is meant by the second appearance.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "SvAjOcJxGlD", "reviewer": "Reviewer_pbvp", "summary": "The manuscript proposes a methodology to speed up structured kernel interpolation Gaussian processes by using a grid with  a number of grid points that does not scale exponentially in the number of data dimension. The sparse grid approach is borrowed from computational mathematics and used to render Gaussian processes more scalable. A number of experiments on datasets with d>>4 illustrates the potential of the method.", "questions": "* The grid point density is not evenly distributed over the input space. Does this imply a preprocessing to the data to concentrate the relevant parts close to the coordinate axis where the density is higher?\n* In Figure 1, you could align the color of the complexity panel and the grid visualisation panels.\n* The term \"nearly linear\" should be properly defined or changed.", "limitations": "Yes.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "1) Clarity\nThe paper is well written. The notation language and figures are concise and well accessible.\nTwo typos: Line 85: scalibility, Line 287: Sentence needs to be revised.\n\n2) Originality\nThe proposed (thinner) grid seems to not have been applied to scalable Gaussian process inference.\n\n3) Significance\nThe method is a proper extension of SKI and is not specific to a particular covariance or likelihood function and can be used in a wide range of inference schemes. As demonstrated in the experiments, SKI is applicable to moderate dimensions up to 10 which was previously challenging.\n\n4) Reproducibility\nAs the code for the experiments is available on https://anonymous.4open.science/r/skisg-485E/README.md, the results should be simple to reproduce. I'm assuming that the code will be shared at a more accessible location upon acceptance.\n\n5) Empirical analysis\nThe experimental section reports a number of results on runtime and memory footprint and demonstrates the utility of the method.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "M8YpNpjhhH-", "reviewer": "Reviewer_vR8u", "summary": "This study presents a way to improve the scalability of Structured kernel interpolation (SKI) for high dimensional dataset.  Specifically,  authors introduce the sparse grid inputs for constructing the inducing kernel Gram matrix $K_{G}$  in approximate kernel $W^{T}K_{G}W$. To use the sparse grid for SKI, they propose the MVM algorithm for the introduced sparse grid to have the near-linear computation for MVM, and the interpolation scheme for computing the $W$ (sparse weight). In this procedure, authors describe the properties of sparse grid and explain how the computational cost induced from sparse grid inputs could be much smaller compared to that of regular grid for SKI. \nThrough experiments, authors validate that using the sparse grid with the proposed schemes is computationally efficient while obtaining the comparable prediction performance.", "questions": "$\\bullet $  Q1. For the scalable GP inference for high dimensional dataset, combining the random Fourier feature with fast food tricks enables the efficient inference in [1].  For the results in table 1, I am quite curious whether the proposed SKI framework could be comparable to other scalable GP inference scheme.\n\n[1] A la Carte – Learning Fast Kernels, AISTATS 15\n", "limitations": "I think that if the validation procedure in experimental section is improved, this work could get a better evaluation.\nThus, I write some concerns about experiment section and some suggestions.\n\n$\\bullet $ For the experiment on synthetic dataset of $\\cos({\\vert{x}\\vert}_{1})$, it seems that the more accurately the kernel is approximated in left figure ($d=6$), the more accurately it seems to predict in right figure ($d=6$). However, what if synthetic set are generated from different process, I am curious whether the effectiveness of sparse grid still remains (for generalization on a variety of datasets). If possible, what about applying the same experiment on some of the UCI datasets used in Table1?  If authors show the consistent results from the high dimensional UCI datasets, this method will be regarded as more powerful.\n\n$\\bullet $ For the UCI regression task, Table 1 does not report the standard deviation of the Test RMSE for 3 trials. Thus, I am not sure whether the proposed scheme results in stable training and prediction performance.  Also, for evaluating the GP regression task, it would be better to report the test negative log likelihood (nll), which is the widely used metric for GP regression task, as done [1].\n\n$\\bullet $ For regression task, comparing the training progress (xaxis: training time, yaxis: test rmse ) over varying number of sparse grid points $|\\mathcal{G}|$  would be better to make sure that the sparse grid results in fast training for high dimensional dataset due to the reduced computational cost as shown in Figure 2 of [2].\n\n\n[1] SKIing on Simplices: Kernel Interpolation on the Permutohedral Lattice for Scalable Gaussian Processes, ICML 21\n\n[2] Product Kernel Interpolation for Scalable Gaussian Processes, AISTAST18\n\n\n ", "rating": 5, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "$\\textbf{Strengths}$\n\nThe existing SKI framework has been the scalable issue for high dimensional dataset, because the computational cost of MVM for the regular grid is $O(N4^{d} + |U| \\log{|U|})$, and this implies that high dimensional data set (large $d$) increases the computational cost exponentially, and thus incurs the computational issue for employing the SKI. This work alleviates this scalable issue by inducing the sparse grid, and then proposing the corresponding MVM for sparse grid inputs. Since the introducing the sparse grids in the naïve manner does not result in the improvement of the computational cost as shown in Figure 1.(c), I think that making the near-linear computation for MVM computation with sparse grid is a quite novel part of this work.\n\n$\\textbf{Weaknesses}$\n\nAlthough this work proposes a novel inference scheme for overcoming the issue of SKI, I think that the validation procedures of the proposed method in the experiment section is not sophisticated, and seems to a little bit weak to validate the generalization of the proposed work on a variety of dataset, as comparing the experiments validating the existing variants of SKI.\nFor further details, see Section Limitations below.\n\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "JPUgdcpuvm", "reviewer": "Reviewer_GXVc", "summary": "This paper uses sparse grids in structured kernel interpolation (SKI), which speeds up Gaussian process inference by replacing the kernal matrix with a kernel matrix evaluated at a structured set of points. In sparse grids the number of points grows faster with the dimension than in full grids typically used in SKI.\n\nThe two main contributions of the paper are (1) linear time algorithm for multiplying sparse grid kernel matrices with vectors that can be used in solving linear systems arising in GP inference via the conjugate gradient algorithm and (2) a simplicial interpolation scheme for approximating the full kernel matrix.\n", "questions": "Which class of kernels does the method apply to? Scattered remarks in the main paper and supplement seem to indicate that the kernel should be a product of one-dimensional stationary kernels but this is not made clear anywhere. In either Section 2 or 3 the reader should be explicitly told which kernels are considered in the paper.", "limitations": "One of the main limitations of the work is that although sparse grids are an improvement over full grids in this regard, the dimension that one can work in still remains fairly low (at most ten in this paper, I think). The authors are sufficiently upfront about this issue.", "rating": 8, "confidence": 4, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Although the idea of replacing full grids with sparse grids is quite obvious and it is rather well known that, as opposed to orders of magnitudes more or something like that, one gains only a few extra dimensions with sparse grids (I appreciate that the paper does not attempt to hide this fact as is evident from line 80 and 272-273), I find this paper a very nice and useful contribution for the field. The paper is very well written and pleasant to read. ", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept