The Price of Linear Time: Error Analysis of Structured Kernel Interpolation
Authors: Alexander Moreno
Justin Xiao
Jonathan Mei
Venue: ICML 2025 poster (2025)
Keywords: Gaussian processes
kernel methods
Primary area: theory->everything_else
TL;DR: We establish the first rigorous error bounds for Structured Kernel Interpolation (SKI), revealing dimension-dependent trade-offs between computational complexity and approximation accuracy.

Abstract: Structured Kernel Interpolation (SKI) scales Gaussian Processes (GPs) by approximating the kernel matrix via inducing point interpolation, achieving linear computational complexity. However, it lacks rigorous theoretical error analysis. This paper bridges this gap by proving error bounds for the SKI Gram matrix and examining their effect on hyperparameter estimation and posterior inference. We further provide a practical guide to selecting the number of inducing points under convolutional cubic interpolation: they should grow as \(n^{d/3}\) for error control. Crucially, we identify two dimensionality regimes for the SKI Gram matrix spectral norm error vs. complexity trade-off. For \(d<3\), \textit{any} error tolerance can achieve linear time for sufficiently large sample size. For \(d\geq 3\), the error must \textit{increase} with sample size for our guarantees to hold. Our analysis provides key insights into SKI's scalability-accuracy trade-offs, establishing precise conditions for achieving linear-time GP inference with controlled error.

The Price of Linear Time: Error Analysis of Structured Kernel Interpolation
Alexander Moreno
Justin Xiao
Jonathan Mei
Abstract
inference require inverse Gram matrix computations and log- determinant calculations. These issues necessitate scalable approximations for larger datasets.
Structured Kernel Interpolation (SKI) (Wilson Nickisch, 2015) scales Gaussian Processes (GPs) by approximating the kernel matrix via induc- ing point interpolation, achieving linear compu- tational complexity. However, it lacks rigorous theoretical eтrоr analysis. This paper bridges this gap by proving error bounds for the SKI Gram ma- trix and examining their effect on hyperparameter estimation and posterior inference. We further provide a practical guide to selecting the num- ber of inducing points under convolutional cubic interpolation: they should grow as n/33 for spec- tral norm error control. Crucially, we identify two dimensionality regimes for the SKI Gram ma- trix spectral norm error vs. complexity trade-off. For < 3, any error tolerance can achieve lin- ear time for sufficiently large sample size. For d 2 3, the error must increase with sample size for our guarantees to hold. Our analysis provides key insights into SKI's scalability aalability-accuracy accuracy trade- offs, establishing precise conditions for achieving linear-time GP inference with controlled error.
Structured Kernel Interpolation (SKI) (Wilson Nickisch, 2015) scales GPs by approximating the kernel matrix via inducing point interpolation. For stationary kernels, it achieves On + m logm)
m) (m inducing points) complex- ity by expressing the kernel via interpolation functions and an inducing point kernel matrix. Despite its effectiveness and popularity (> 600 citations: (Gardner et al., 2018) has 3.5k GitHub stars), SKI lacks theoretical analysis. Key ques- tions include: for a fixed SKI Gram matrix error using cubic convolutional interpolation, how many inducing points needed? If m is a function of when does O + m log m) complexity remain linear? What are the implications for hyperparameter estimation and posterior inference?
This paper bridges SKI's practice- practice-theory gap with: 1) The first error analysis for the SKI kernel, Gram matrix (spec- tral norm), and related quantities, yielding guidelines for selecting inducing points (m & n²/3 for spectral norm error control). 2) SKI hyperparameter estimation analysis. 3) SKI inference analysis. Key findings are: a) Two dimension- ality regimes link SKI Gram matrix error and complexity: for d < 3, any fixed spectral error is achievable in lin- ear time (sufficient n); for d 2 3, the error must increase with n for our linear time guarantee. b) For ju-smooth log- likelihoods, SKI-based gradient ascent approaches true stationary point's neighborhood at an O(1/K) rate; neigh- borhood size depends on SKI score function error (ignoring responses, it scales linearly with n if m xx d/33). This lever- ages Stonyakin et al. (2023)'s ' results on inexact gradient descent (d' Aspremont, 2008: Devolder et al., 2014).
1. Introduction
Gaussian Processes (GPs) (Rasmussen & Williams, 2006) are vital stochastic processes in machine learning and statis- tics, with applications including spatial data analysis (Liu & Onnela, 2021), time series forecasting (Girard al., 2002), bioinformatics (Huang et al., 2023), and Bayesian optimiza- tion (Frazier, 2018). GPs provide non-parametric frame- work for modeling distributions over functions, offering flexibility and uncertainty quantification. Their ability to incorporate prior knowledge via kernel choice makes GPs effective for regression and classification.
Sec. reviews related work and Sec. 3 gives SKI back- ground. Sec. 4 bounds key SKI errors (kernel, Gram, cross- kernel matrices). Sec. 5 analyzes SKI MLE and posterior errors. Sec. 6 presents empirical validation. Sec. 7 con- cludes with results, limitations, and future work.
However, GPs face substantial Onn²) (sample size n) com- putational and O(n²) memory bottlenecks, training and
2. Related Work
'MBZUAI Darwin Al IIonQ. Correspondence to: Alexander Moreno <alexander.f @ gmail.com>>
Related works fall into three main groups: theoretical analy- ses of Gaussian process regression or kernel methods using approximate kernels, SKI and its extensions, and papers
Proceedings of the ndnd International Conference Machine Learning. Vancouver, Canada. PMLR 267, 2025. Copyright 2025 by the author(s).

---

Error Bounds for Structured Kernel Interpolation
Quantity
Bound
SKI kernel error
37d
SKI Gram matrix error
0
ne
maya
SKI cross-kernel matrix error
0 mmax
c2d v/Pn?
SKI n -normalized score function error
ma/a
SKI posterior mean error
V VnO
max (T, n.) + VTnMa
omazza
SKI posterior covariance error
0
Tn" me" + TTmee" max(T,n)
m3ra
Table 1. Summaay of SKI error bounds with convolutional cubic interpolation. Variables: n. (urain/test sizes), (dimensionalityy. m (inducing points), (hyperparameters), c. > 0 (constant). Key: Gram matrix error is Onn m -33d))
developing techniques we leverage. In the first group, Burt et al. (2019; 2020) analyzed the sparse variational GP frame- work (Titsias, 2009; Hensman et al., 2013), deriving KL divergence bounds between true and variational approxi- mate posteriors. Moreno et al. (2023) gave bounds on SKI Gram matrix approximation error for univariate features, comparing to Nyström, but did not analyze downstream ef- feets approximate MLE or GP posteriors. Also, Wynne & Wild (2022); Wild et al. (2021) respectively gave a Banach space view of sparse variational GPs and Nyström connec- tions. Finally, Modell provide entry-wise error bounds for low-rank kernel matrix approximations; our approach also uses entry-wise bounds, but theirs target the best low-rank approximation, while ours are for the SKI Gram matrix. Only Moreno al. (2023) treated SKI specifically, and only in very special case setting. Our work gives end-to- end theory of SKI from elementwise to spectral error to estimation and posterior inference.
3. Gaussian Processes, Structured Kernel Interpolation and Convolutional Cubic Interpolation
This section provides background on Gaussian Processes (GPs) and two key techniques for enabling scalable infer- ence: Structured Kernel Interpolation (SKI) and Convolu- tional Cubic Interpolation. SKI (Wilson & Nickisch, 2015) addresses GPs scalability issue by approximating the kernel matrix through interpolation on a set of inducing points, leveraging the efficiency of convolutional kernels. partic- ular, cubic convolutional kernels, as detailed in Keys (1981), provide smooth and accurate interpolation scheme that forms the foundation of the SKI framework. In this paper, we focus on this cubic case as it is used by SKI. Future work may extend this to study higher-order interpolation methods. Here, we formally define these concepts and lay the groundwork for the subsequent error analysis.
In the second group, Wilson & Nickisch (2015)'s founda- tional work, which we analyze, introduced SKI for scalable large-scale GP inference. Kapoor et al. (2021) extended SKI to high dimensions via the permutohedral lattice. Ya- dav et al. (2022) developed a sparse grid kernel interpolation approach to address dimensionality. Most recently, Ban et al. (2024) proposed flexible SKI adjusting grid points based on kernel byperparameters. We focus our analysis on the orig- inal (Wilson& & Nickisch, 2015) technique, though future work could extend these latter papers' settings.
3.1. Gaussian Processes
A Gaussian process ~ GP(v, ki)) is a stochastic process (5(x))xex such that any finite subcollection [xx,))r: is multivariate Gaussian distributed. We assume that we have index locations x; E and observations yi R for set of training points = 1,..., n such that
= {xx) + E, ~(0
where v x kp X R are the prior mean and covariance functions, respectively. with positive semi-definite (PSD) kernel with hyperparameters O. Given (x,, 14 3 or alternatively X E ixd ,y E R". we wish to: 1) estimate hyperparameters 0 E O C RD of kernel ko (c.g. RBF kernel) 2) do Bayesian inference for the posterior mean (4()) E IR. and covariance Σ(.) E T at a set of test points {xe)²= Assuming v = (a mean-zero GP prior), for 1), one maximizes the logelikel
Relevant too are papers whose results or proof techniques we leverage or extend. We derive required multivariate extension to Keys (1981)'s error analysis for convolutional cubic interpolation. We also use a recent result from the inexact gradient descent literature (Stonyakin et al., 2023) to analyze gradient ascent on the SKI liglliieidood versus the true loglliiel Finally, we adapt proof technique (Bach, 2013; Musco & Musco, 2017), common for approxi- mate kernel ridge regression in-sam error, to bound test SKI mean function error.
1 " T(K +021)
C(0: x)
log K #
2'|
" log
(1)
2

---

to find 0 E D C θ where K R with entries ij ko (x,, x,) is the Gram matrix for the training dataset. For 2), given the kernel function, known observation variance o² > O and matrix of kermel evaluations between test and training points x E TRTxx xn the posterior mean and covariance are Error Bounds for Structured Kernel Interpolation
To learn kernel hyperparameters, one maximizes the SKI log-likelihood approximation (henceforth the SKI log- likelihood)
1 [(0: - 2 (K+
(4()) = x (K (K+ + +00))= (() K. + 1 - KK -o
(2)
log K o
log (2T)
(3)
Given the SKI kernel 1. X x X + R with learned hy- perparameters, one can do posterior inference of the SKI approximations to the mean (.) and covariance E() at a set of T test points as
Intuitively, the GP prior represents our belief about all pos- sible functions before seeing any data. When we observe data points, the posterior represents our updated belief it gives higher probability to functions that fit our observations while maintaining the smoothness properties encoded in the kernel. The posterior mean can be viewed as weighted average of these functions, where the weights depend on how well each function fits the data and satisfies the prior assumptions. The posterior variance indicates our remaining uncertainty it is smaller near observed points where we have more confidence, and larger in regions far from our data.
ia()) K, [K + 1091 y Σ(.) K.K+t + K + + 'Kx,
where K, x R is the matrix of SKI kernels between test points and training points and K, E RRTX xT is the SKI Gram matrix for the test points. Going forward, we may write (()) and Ž(0), dropping the explicit dependence on the data but implying it.
challenge that, between the log-likelihood and the pos- teriors, one first needs to compute the inverse regularized Gram matrix times the response vector, (K 12) y. Sec- ond, one needs to compute the rog-determinant log | I. These are both (n²²) computationally and Onn²) memory.
3.3. Convolutional Cubic Interpolation
Convolutional cubic interpolation (Keys, 1981) gives a con- tinuously differentiable interpolation of a function given its values on a regular grid, where its cubic convolutional kernel is piecewise polynomial function that ensures continuous differentiability. We formalize this using the definitions of the cubic convolutional interpolation kernel and tensor- product cubic convolutional function below. We also define an upper bound for the sum of weights for each dimension, which will be a useful constant going forward.
3.2. Structured Kernel Interpolation
Structured kernel interpolation (Wilson & Nickisch, 2015) or (SKI) addresses these computational and memory bot- tlenecks by approximating the original kernel function kg X x X + R, X S Rd by interpolating kernel val- ues at chosen set of inducing points
Definition 3.1. The cubic convolutional interpolation kernel u : R + R is given by
u
U
E
R md
|s|3
|11/ + 1, 3
0 18
<1, <2
um
u(s))
|s|3
|3/2 -4/s|
+
1 1 |s
The approximate kernel 7: X x X- - IR then:
0,
|8| 2 2.
k(x, = w(x) (x'(x')
Definition 3.2. Let x = (I), F2, x2) E Rd be a dimensional point and ʃ Rd - a function defined on a regular grid with spacing h in each dimension. Let ex denote the grid point closest to x. The tensor-product cubic convolutional interpolation function g Ridd + R is:
where Ky E Rmxm the inducing point kernel matrix, and w(x), w(x') R" are vectors interpolation weights using (usually cubic) convolutional kernel RR IR for the points and x', respectively. One then forms the SKI Gram matrix K- = W with W sparse matrix of L interpolation weights per row for a polynomial of degree / - 1. By exploiting the sparsity of each row, for stationary kernels this leads to a computational complexity of OnL + m log m) and a memory complexity of (nI.).
g(x)
Σ kE(-1,0,1,34
f(ex hik)
II j=1
[x (cx); loki u
where u is the cubic convolutional interpolation kernel and k = (kI,..., ka) is a vector of integer indices.
3

---

Error Bounds for Structured Kernel Interpolation
Definition 3.3 (Bound on interpolation weights). Let U C R be the (uniform) inducing grid with spacing h > 0, and let IR + I be the I-D interpolation kernel. Define
with the number of inducing points m. but may exhibit exponential error growth with increasing dimensions.
Its proof proceeds by induction on the dimension d. The base case for d = relies on the known univariate cubic convolution interpolation error bound from Keys (1981). The inductive step then demonstrates how this error bound accumulates, approximately multiplicatively with the con- stant с for each added dimension, when extending to the tensor-produs case.
Σ
I
c == sup
< 00.
h
ZER
ZEU N TT ) 40
For cubic interpolation u has compact support (only the four nearest grid points per dimension have nonzero weights), so c < DO. In d dimensions with tensor-product weights w(x) over Ud, note that each row has at most L = 44 nonzeros and
Lemma 4.1. The uniform error bound over compact do- main X' C RR for tensor- massorprodduc cubic convolutional in- terpolation of a thrice continuously differentiable function f :R + IR by its interpolant g .R. - R is
[|w(x)]l S da for all Rd
Going forward, we always assume that we use tensor prod- uct convolutional cubic polynomial interpolation, so that = as in (Wilson & Nickisch, 2015), but that we may vary the number of inducing points m. In particular, we will analyze how the number of inducing points affects error for different terms of interest, and how to choose the number of inducing points.
lf(x)) - -=(x) = ((C
sup =EX
equivalently O
cd
m37a
Proof. See Appendix B.I.I.
The following Lemma allows us to bound the absolute dif- ference between the true and SKI kernels uniformly with the same big-O error as for the underlying interpolation itself. The proof uses the the triangle inequality to decompose the error into two parts: the first is from a single interpolation, while the second from nested interpolations.
4. Important Quantities
This section derives bounds for key quantities in Structured Kernel Interpolation (SKT). Section 4.1.1 provides bound on the elementwise error between the true kernel and its SKI approximation. In Section 4.1.2, we extend this to the spectral norm error of the SKI approximation for the training Gram matrix and train-test kernel matrix. Finally, in section 4.2 we present conditions on the number of inducing points for achieving specific error tolerance e >> 0 and error needed to guarantee linear time complexity, noting linear time always holds for d with sufficiently large samples.
Lemma 4.2. Let 8m, denote the uniform error bound of tensor-product cubic convolutional interpolation (Lemma 4. I) for m inducing points and interpolation degree T_1/d (hence L = 4 for cubic interpolation, as used in SKI). The SKI kernel k. : X x X -> R approximating a thrice continuously differentiable (not necessarily PSD) kernel dimension has error
k X x X - R with uniform grid spacing h in each
4.1. Error Bounds for the SKI Kernel
[x(x,xx) - x(xx,x)) S 58,.L + A6m.L
This subsection analyzes the error introduced by the SKI approximation of any symmetric kernel function. We start by extending the analysis of Keys (1981) to the multivariate setting, deriving error bounds for multivariate cubic convolu- tional polynomial interpolation. We then use these to derive the elementwise error for the SKI approximation k((x, x') to a symmetric kernel, which may not be positive (semi)- definite. We next apply these elementwise bounds to derive spectral norm error bounds for SKI kernel matrices, which will be crucial for understanding the downstream effects of the SKI approximation Gaussian process hyperparameter estimation and posterior inference.
( c2d
= 0
m3/a
Proof. See Appendix 1.1.
4.1.2. SPECTRAL NORM ERROR
We now transition from elementwise error bounds to spec- tral norm bounds for the SKI gram matrix's approximation error, finding that it grows linearly with the sample size, ex- ponentially with the dimension, and decays as m -33/ with the number of inducing points. This is not only of indepen- dent interest but also important for nearly all downstream estimation and inference analysis. We also provide a bound on the spectral norms of the SKI train/test kernel matrix's 's approximation error. This is useful when analyzing the GP posterior parameter error.
4. 1.1.1 ELEMENTWISE
Our first lemma shows that multivariate tensor consor-productee cubic convolutional interpolation retains error cubic in the grid spacing of Keys (1981), which is equivalent to m -3/d decay
4

---

Error Bounds for Structured Kernel Interpolation
For this next lemma will express it both in the general interpolation setting and again give the specific big-O for convolutional cubic interpolation, but going forward we sometimes only show the latter setting in the main paper and derive the general settings in the proof. In particular, whenever we use big O-notation we are assuming convolu- tional cubic interpolation.
norm error derived in Proposition 4.3: this is sufficient condition for the desired spectral norm error to hold. It then relates the number of inducing points to the grid spacing in the SKI approximation, assuming regular grid with equal spacing in each dimension. By substituting this relationship into the sufficient condition, the proof derives the sufficient number of inducing points to control error.
Proposition 4.3. For the SKI approximation K of the true Gram matrix K, we have
Theorem 4.5. If the domain is [-D,D], then to achieve a spectral norm error of | K||2 S c, it is sufficient to choose the number of inducing points m such that:
|K-
K1 S n (6m,L + ct L.)
mi = (K)) d/3
=
ne?d
0
m//d
for some constant K' that depends only on the kernel func- tion and the interpolation scheme.
Proof. See Appendix B.1.4
Proof. See Appendix B.2.1.
Its proof leverages the previous elementwise error bounds. Since the error K - K is symmetric, its spectral norm is bounded by its maximum absolute row sum (I) ol norm). Each row sum involves n terms, each bounded by the ele- mentwise error: this gives n (elementwise error) structure. The next bound uses the property |All2 I ||A||||| which simplifies to ||A||2
The number of inducing points should thus grow:
Sub-linearly with the sample size and decrease in error for d < 3, linearly for d 3 and super-linearly for d > 3. Thus, as we want a tighter error tolerance or have more observations we need more inducing points, but at very different rates depending on the dimensionality.
2 S max([|A|]: ||A|| ox) for poten-
tially non-square matrices. The maximum absolute col- umn sum (II (]1) and row sum (II >)>) of the enror matrix K.x - k_x are bounded using the elementwise error from Lemma 4.2, scaled by T and n respectively.
Linearly with the volume of the domain (2D) Thus, if our observations are concentrated in a small region and we select an appropriately sized domain to cover it need fewer inducing points.
Lemma 4.4. Let Kx E Txx be the matrix of kernel evaluations between T test points and n training points, and let xxxxx x E Then
T' be the corresponding SKI approximation.
Exponentially with d², as we have a 2d term taken the power d/3. However, our empirical results suggest that this is quite pessimistic.
max(n, cd
|K - K. O
3/d
Proof.
See Appendix B.1.5.
The next Corollary establishes a condition on the spectral norm error, e, that ensures linear-ting O(n) computational complexity for SKI. The core idea that should be such that if we choose m based on the previous Theorem, (n log n) and thus log = O(n).
4.2. Achieving Errors in Linear Time
Here, show how many inducing points m sufficient to achieve desired error tolerance >> 0 for the SKI Gram matrix when using cubic convolutional interpolation. Based on the Theorem, we should grow the number of inducing points at an 4/'3 rate. We then show corollaries describing 1) how e and m must grow to maintain linear time 2) how the dimension affects whether the error must grow with the sample size ensure linear time SKI.
Corollary 4.6. If
e> +1+ e K' 2d D3
n(log 3/dd
(4)
C3/d
n3/d
for some constants K' K'C that depend on the kernel function and the interpolation scheme and we choose m. >> based on the previous theorem, then we have both ||
The following theorem shows the number of inducing points that will control Gram matrix error. It says that the number of inducing points should grow as d/3 to achieve a fixed error. The proof starts by lower bounding the desired spec- tral norm error with the upper bound on the actual spectral
K 2 e and SKI computational complexity of O(n).
Proof: See Appendix B.2.2.
5

---

Error Bounds for Structured Kernel Interpolation
Interestingly, the previous Theorem and Corollary im- ply fundamental difference between two dimensionality regimes. For d < 3, the choice of m required for a fixed grows more slowly than n/ log n. This means that for any fixed > 0, SKI with cubic interpolation is guaranteed to be a linear-time algorithm for sufficiently large n. In con- trast, for d 2 3, the choice of m required for a fixed >> eventually grows faster than n/ log n. Thus, to maintain linear-time complexity for d 2 3 and the guarantees from Theorem 4.5, we must allow the error e to increase with n. This demonstrates that the curse of dimensionality impacts the scalability of SKI, making it challenging to guarantee both high accuracy and linear-time complexity in higher dimensions. The next corollary formalizes this,
convergence to neighborhood of a critical point, only that at least one iterate will approach it. Given the challenges of non-concave optimization and the fact that we leverage fairly recent result, we leave stronger results to future work.
Let D C θ be a compact subset that we wish to optimize over. In the most precise setting we would analyze projected gradient ascent, but for simplicity we analyze gradient as- cent. Let kg X x X + E. be the SKI approximation of ko X' x X - R using m inducing points and interpolation degree L - 1. We are interested in the convergence proper- ties of inexact gradient ascent using the SKI leeliielioood, e.g.
0 0--1= 0x + nVL((u),
Corollary 4.7. For d < 3. for any e > O, Corollary 4.6 holds for any n sufficiently large, so that choosing m based on Theorem 4.5 is sufficient to achieve linear complexity. For d 2 3, c must grow with the sample size to guarantee linear complexity using Theorem 4.5.
where 7) E R is the learning rate and V (0x))) is the SKI score function (gradient of its log-iikelihodd. We assume: 1) a p-smooth log- oogllikelhhod. If we optimize on a bounded domain, then for infinitely differentiable kernels (c.g. RBF) this will immediately hold. 2) that the kernel's partial deriva- tives are themselves symmmetric (not necessarily PSD) ker- nels.
Proof. For d < 3, the RHS of Eqn. 4 decreases with with limit 0 and thus for sufficiently large sample size will be S satisfying the conditions to guarantee small error and linear time. For d 2 3, the RHS of Eqn. 4 grows with n,, so that must grow to satisfy the conditions for the guarantee. N
The true log- likelihood is p-smooth over D. That is, for all 0, 0' E D,
Assumption 5.1 (p-smooth-log--ikellinnn
V V s ul10 0'll
Assumption 5.2. (Kernel Partials) For each E (1,....p) p)) the partial derivative of ky with respect hyperparam- eter 02 E R. denoted as K/, (x, = 20,
5. Gaussian Processes Applications
@Ra(a.n") is also a
In this section, we address how SKI affects Gaussian Pro- cesses Applications. In Section 5.1 we address how using the SKI kernel and log- likeliiel affect hyperparameter estimation, showing that gradient ascent on the SKI log- likelihood approaches a ball around a stationary point of the true liglliillhood. In section 5.2 we describe how using SKI affects the accuracy of posterior inference.
symmetric (not necessarily positive semi-definite) kernel.
We next state results leading to our SKI score funetion error bound. We argue that we can apply the same elementwise error we derived previously to the SKI partial derivatives.
Lemma 5.3. [Bound on Derivative of SKI Kernel Error using Kernel Property of Derivative] Let Kl, (z, x²) be the SKI approximation of ko,, (I,x²), using the same inducing points and interpolation scheme as kg. Then, for all x, x' X' and all 0 E O, the following inequality holds:
5.1. Kernel Hyperparamete Estimation
Here we show that, for y-smooth log- likelihood, an iterate of gradient ascent on the SKI log-likelihood approaches neighborhood of stationary point of the true log-likelihood at O (+) rate, with k = 1, I,...,A K the iterations and with the neighborhood size determined by the SKI score function's error. To show this, we leverage a recent result for non-convex inexact gradient ascent (Stonyakin et al., 2023), which requires an upper bound on the SKI score function's is error. This requires bounding the spectral norm error of the SKI Gram matrix's partial derivatives. In order to obtain this, we note that for many symmetric kernels, under weak assumptions, the partial derivatives are also (not necessarily PSD) symmetric kernels, and thus we can reuse the previous results directly on the partial derivatives.
Oko(2, 00,
ako(22 x²)
= 1/,, (2, - (x()
00,
S 8'm, + +c &.m.I
2d
m33dd
where 8m.L is an upper bound the error of the SKI ap- proximation of the kernel ko, (x,,²)) with inducing points and interpolation degree L - 1, as defined in Lemma 4.2.
Note that Stonyakin et al. (2023) does not actually imply
Proof. See Appendix C.I.1
6

---

Error Bounds for Structured Kernel Interpolation
We then use the elementwise bound to bound the spectral norm of the SKI gram matrix's partial derivative error. This again leverages Proposition 4.3, noting these partial deriva- tives of the Gram matrices are themselves Gram matrices.
Proof. See Section C.I.3.
We apply Stonyakin et al. (2023) below: the result is the same in their paper (and assumes y-smoothness as we did C). but using gradient ascent instead of descent and using the score function error above. It says that at an O (K rate, at least one iterate of gradient ascent has its squared gradient norm approach a neighborhood proportional to the squared SKI score function's 's spectral norm error.
Lemma 5.4. [Purtial Derivative Gram Matrix Difference Bound] For any 1 E (I....,P).
0K
OK
S L.
Theorem 5.6. (Stonyakin et al., 2023) For inexact gradient ascent on L with additively inexact gradients satisfying [V Vā(0))l S Eg. we have:
00,
00:
2
nc2d
=0
m3/d
where L. J. is the bound on the spectral norm difference between the kernel matrices corresponding to kg, and its SKI approximation kep, (analogous to Proposition 4.3, but for the kernel k-)))
VL(00)) 2p(c" -[(0)))
2
min
+
(5)
K
2/
where L is the value at stationary point, L(0o) is the initial, function value, K is the number of iterations and Eg is the gradient error bound in the previous Lemma.
Proof. See Section C.1.2.
We now bound the SKI score function. The key insight to the proof is that the partial derivatives of the difference between regularized gram matrix inverses is in fact a dif- ference between two quadratic forms. We can then use standard techniques (Horn & Johnson, 2012) for bounding the difference between quadratic forms to obtain our result. The result says that, aside from the response vector's norm, the error grows quadratically in the sample size, at a square root rate in the number of hyperparameters and exponen- tially in the dimensionality. It further decays an må rate in the number of inducing points. Noting that maintain Gram matrix error, should grow at an d/3 rate, have that if lyyle O()n), the n normalized score function error can grow linearly with the sample size when choosing the number of inducing points based on Thcorem 4.5. To control it, we actually want m nd2//3
5.2. Posterior Inference
Finally, we treat posterior inference. As our hyperparameter optimization results only say that some iterate approaches a stationary point, we will focus on the error when the SKI and true kernel hyperparameter match. We add an assumption
Assumption 5.7. (Bounded Kernel) Assume that the true kernel satisfies the condition that (x)xx')) s -M for all x, x' E X.
Now we bound the spectral error for the SKI mean func- tion evaluated set of test points. The proof follows standard strategy commonly used for approximate kernel ridge regression. See Bach (2013); Musco & Musco (2017) for examples. The result says that the 12 error may grow exponentially in the dimensionality, super-linearly in the training sample size and at worst linearly in the test sample size. It decays at m rate in the number of inducing points. For controlled error want 2d/3
Lemma 5.5. [Score Function Bound (Final Version)] Let voL(0) and voL(0) be the true and SKI tog-likelihood gradients. Assume the kernel is sufficiently regular such that its partial derivatives are bounded by || I1 for a constant C. (n)).
Lemma 5.8. (SKI Posterior Mean Error) Let p(()) be the GP posterior mean at set of test points E Txd and ja(.) the SKI posterior mean at those points. Then the SKI posterior mean 12 error is bounded by:
Then the spectral norm of the score function error is bounded as:
[1ī((-) - -10())12
()Vo(()) voL(0)12 (
m3/d (1wl12n² n²)
VP
Jmmajd Vno maxx(T, m) + VTnMM o²
2d
Furthermore, under the common assumption that |1y113 Onn). the n "nnormalized error is bounded as:
Proof. See Appendix C.2.1.
We now derive the spectral error bound for the test SKI covariance matrix. The proof involves noticing that a key term a difference between two quadratic forms, and using
1 (oC(O) -o2(0)2==
-020V n
m3/dd

---

Error Bounds for Structured Kernel Interpolation
standard techniques for bounding such difference. The result shows that the error grows at worst super-linearly but subquadratically in the number of test points, quadratically in the training sample size and exponentially in the dimen- sion. The error scales with the number of inducing points at an m 1-3dd rate, so that it decays if < 3. If we select the number of inducing points to be proportional to nd/3,, then the error grows at rate n +d/3 for d 3, so that we do nor have error control for the covariance, despite having it for the Gram matrix. A future question is whether this bound can be improved.
inducing points Miotal (with Mpd per dimension, so Miotal = -mpa) and measured || || - K||2.
2. Error vs. n (Scaled myotal): We varied n (50 to 1000) and scaled Msotal k f--3/3 /3/3 (k- = 1.0), setting mpd max(4, ronnd (knddkk//////33)) We then measured
|| K|12.
6.2. Results and Discussion
The results are presented in Figure 1.
Lemma 5.9. [SKI Posterior Covariance Error] Let E(.) be the GP posterior covariance matrix at set of test points E RR xd and E() be its SKI approximation. Then
6.2.1. ERROR SCALING WITH NUMBER OF INDUCING POINTS (Mootal)
Figure I(a) displays spectral norm error vs. Mittall (log-log scale, fixed n).
20) 201l
VTMM
S m. +
o²
max m.
For d = 1, 1, error shows clear power-law decay with Mootal (linear trends in log-log plot).
+ In.mm. Tnme 2
TTnmc2dM M
Observed decay rates (slopes) are consistently steeper than the predicted - ): for d 1, slope -5.11 (expected -3.00); for d 2, 2288 (ex- pected 1.500); for d = 3, ~ 2205 (expected -1.00). 3/d
max(TT,m... m.).
0²
+ Timmal max(T,
=0
m3/d
This suggests the || K|
O /modd//mocd /m
bound, while holding, may be pessimistic for smooth RBF kernels and uniform data, where actual error de- cay is faster.
where PT.,mL is defined as in Proposition 4.3.
Proof. See Appendix C.2.2
Error magnitude clearly increases with d for fixed Miotal (e.g., at Mtotal & 64, error d = 3 > = 2 > d = 1). This aligns with the 2dd term in the bound, highlighting increased approximation challenge higher for constant meotal.
6. Empirical Analysis
To empirically investigate derived theoretical error bounds for the Structured Kernel Interpolation (SKI) Gram ma- trix approximation, we conducted numerical experiments. Experiments focused spectral norm error || K|12 be- havior function of data points (n), total inducing points (Mosta)) and dimensionality (d).
6.2.2. ERROR SCALING WITH SAMPLE SIZE (n) AND Mettal X nd/3
Figure 1(b) shows spectral norm error vs. n with Mootal R n//3
6.1. Experimental Setup
For d 1,2,3, eгror decreases with increasing n when Mostal is scaled thusly-a favorable finding. The- ory (substituting midel x n into the O(ne2d/n //m/d 3/d bound) suggests error stabilization around ((22d) (con- stant w.r.t. n).. The observed decrease indicates better practical performance of this scaling.
Experiments used synthetic datasets with input points z; E Rd drawn uniformly from [0, 10". We employed a stan- dard Radial Basis Function (RBF) kernel for all tests. SKI used cubic convolutional interpolation, consistent with our theoretical focus. All computations were performed using PyTorch and GPyTorch.
Though higher d (e.g., d 3) show larger error for smaller errors for d = 2, converge to similar low values as n + 1000. This noteworthy convergence implies the recommended mista scaling effectively mitigates initial error penalties from higher d, enabling comparable accuracy for larger datasets.
Two primary sets of experiments were run for dimensionali- ties d E {1,2,33
1. Error vs. Mistal (Fixed n): For fixed n (1000 for d 1, 500 for = 2, 250 for d = 3), we varied total
8

---

Error Bounds for Structured Kernel Interpolation
S0 Grom Motrll Erorvs, n 16tt thes (7")
56 Cram Mirri Erte vi Plg If and ri
or 10
dna Sod Fior 2280 Do---- 1.50
SerTart Fr -211
Hinde HOlDE
(a) Error vs. Meotal (Fixed m)
(b) Error vs. n (with Mastal OC n 2//3)
Figure 1. Empiriaal SKI Gram matrix spectral norm error K R2 (a) Error vs. Mestal (Fixed n): Error decay as Meotal increases
(fixed 1000 for d = 1, 500 for d = 2. 250 for d = 3). Log-log slopes (e.g.. = -5.11 for d = 1. 22.8 for d 2, ~ -2.0 for d = 3) are steeper than theoretical -3/d (-3.00, 30, -1.50, -.00 resp.). suggesting faster practical RBF error decay, Error magnitude increases with for given Mectat- (b) Error vs. (Mestat x n) Error as varies, with Mootal n//3 (mpd R n/3/3, so Meotal & n n For d = 1,2,3. 3, error decreases as n increases. Error levels for different d converge for larger n (e.g., n R 1000), showing scaling effectiveness against dimensionality impact.
6.3. Summary of Empirical Findings
with < 3 somewhat freely, but to use more care in the accuracy/speed tradeoff for d 3.
Empirical results largely support our theoretical framework. SKI Gram matrix approximation error is consistent with identified dependencies on n, Miotaal, and d. Scaling induc- ing points Mosial & 4/33 is highly effective, controlling and actively reducing approximation error as grows across tested dimensions. Observed Mittal еπтог decay rates for RBF kernels (which are entire, the highest level of smooth- ness) are faster than theoretical bounds suggest, indicat- ing conservative bounds and SKI's strong performance for smooth kernels. Error level convergence across dimensions for meotal n/ scaling larger is particularly encour- aging for practical SKI application.
This work could be extended by analyzing the error of SKI with other interpolation schemes e.g. Lagrange interpo- lation (Stoer et al., 1980), using potentially higher order polynomials. We could then analyze how to vary m and L jointly. Further, we could extend to the setting where hy- perparameters (e-g. lengthscale) vary with the sample size n. Additionally, one could analyze the error of SKI in more complex settings, such as when the inducing points are not placed on regular grid (Snelson & Ghahramani, 2006) or for non-stationary kernel functions, in which case the com- putational complexity would no longer be 0nL+ log m)) Further, we analyze the optimization properties under gra- dient ascent: would be interesting to analyze it under stochastic gradient ascent, analogous to Lin et al. (2024), but now using inexact noisy SKI gradients. Finally, one could analyze the methods for extending SKI to higher di- mensions (Kapoor et al., 2021; Yadav et al., 2022) and for faster SKI inference (Yadav et al., 2021).
7. Discussion
In this paper, we provided the first rigorous theoretical anal- ysis for structured kernel interpolation. takeaway is that control the SKI Gram matrix's spec- tral norm error, the number of inducing points should grow as n 4/3, Additionally, we showed the spectral norm error of the SKI gram and cross-kernel matrices, and how this impacts achieving a specific error in linear time. We then analyzed kernel hyperparameter estimation, showing that gradient ascent has an iterate approach ball around sta- tionary point, where the ball's radius depends the spectral error of the SKI score function. We showed that m = 24/3 sufficies to control the 1 normalized score function error. We concluded with analysis of the error of the SKI posterior mean and variance, where m = 222//3 controls error for the mean function. For practitioners, the takeaway is to use SKI
key practical
This paper used reasoning LLMs particularly Gemini Pro 2.0. The authors contributed the paper idea and early kernel and spectral norm eгror analysis. LLMs helped outline the statements to be made, turn initial rough descriptions into more formal language, and make proof attempts. In general, LLM proof attempts were wrong, but could drive insights into a working proof strategy. We also used versions with in- ternet access to help bring up relevant papers: hallucination rates were moderate.
9

---

Error Bounds for Structured Kernel Interpolation
Impact Statement
trol. IEEE transactions on pallern analysis and machine intelligence, 37(2):408-44 423, 2015.
This work contributes to a deeper theoretical understanding Structured Kernel Interpolation (SKI) (Wilson & Nickisch, 2015) for Gaussian Processes (GPs). By establishing error bounds and analyzing the impact of SKI on hyperparameter estimation and posterior inference, this research can lead to more confident use of approximate Gaussian Processes. These models have broad applications in various domains, including those mentioned in the introduction as well as robotics (Deisenroth et al., 2015), environmental model- ing (Desai al., 2023), and healtheare (Alaa & van der Schaar, 2017). Improved Gaussian Process models can hance prediction accuracy and decision-making potentially leading to advancements in robotics, more accurate envi- rommental predictions, and better healthcare outcomes. It is important to acknowledge that the application of Gaussian Process models also carries potential risks. For instance, in healthcare, inaccurate predictions or biased models can lead to misdiagnosis inappropriate treatment (Morley al., 2020). Therefore, understanding potential sources of eиrоr when using approximations can be crucial to understanding how reliable we can expect them to be.
Desai, A., Gujarathi. E., Parikh, s., Yadav, S., Patel, Z., and Batra, N. Deep Gaussian processes for air quality inference. In Proceedings of the 6th Joint International Conference on Data Science & Management of Data (10th ACM IKDD CODS and 28th COMAD), pp. 278- 279, 2023.
Devolder, O., Glineur, F., and Nesterov, Y. First-order methods of smooth convex optimization with inexact oracle. Mathematical Programming, 146(1):37-75, 2014.
Frazier, P. I. A tutorial on Bayesian optimization. arXiv preprint 1807.0:8 2018.
Gardner, J., Pleiss, G., Weinberger, K. Q., Bindel, D., and
Wilson, A. G. GPytorch: Blackbox matrix-mmtrix Gaus- sian process inference with gpu acceleration. Advances in neural information processing systems, 31, 2018.
Girard, A., Rasmussen, C., Candela, Q., and Murray- Smith, R. Gaussian process priors with uncertain inputs application to multiple-ste steipppeeetepppeee ahead time series forecasting. Advances in neural information processing systems, 15, 2002.
References
Alaa, A. M. and van der Schaar, M. Bayesian inference of in- dividualized treatment effects using multi-task Gaussian processes, Advances in Neural Information Processing Systems, 30, 2017.
Hensman, J., Fusi, N., and Lawrence, N. D. Gaussian processes for big data. In Uncertainty in Artificial Intelli- gence, 2013.
Bach, F. Sharp analysis of low-rank kernel matrix approxi- mations. In Conference on learning theory, pp. 185-209. PMLR, 2013.
Horn, R. A. and Johnson, C. R. Matrix analysis. Cambridge university press, 2012.
Huang, D., Jiang. J., Zhao, T., Wu, S., Lyu, Y., Feng, J., Wei, M., Zhu, Z., Gu, J., et al. DiseaseGPS: auxiliary diagnostic system for genetic disorders based on genotype and phenotype. Bioinformatics, 39(9)bbtad517, 2023.
Ban, H., Riemens, E. H., and Rajan, R. T. Malleable kernel interpolation for scalable structured Gaussian process. In 2024 32nd European Signal Processing Conference (EUSIPCO), pp. 1001 IEEE, 2024.
Burt, D., Rasmussen, C. E., and Van Der Wilk, M. Rates of convergence for sparse variational Gaussian process re- gression. International Conference on Machine Learn- ing, pp. 862-871. PMLR, 2019.
Kapoor, S., Finzi, M., Wang. K. A., and Wilson, A. G. Skiing on simplices: Kernel interpolation on the permuto- hedral lattice for scalable Gaussian processes. In Interna- tional Conference Machine Learning, pp. 5279-5289. PMLR, 2021.
Burt, D. R., Rasmussen, C. E., and Van Der Wilk, M. Con- vergence of sparse variational inference in Gaussian pro- cesses regression. Journal of Machine Learning Research, 21(131):1-63, 2020.
Keys, R. Cubic convolution interpolation for digital image processing. IEEE transactions acoustics, speech, and signal processing, 293))1153-1116,, 1981
d' ''spremont, A. Smooth optimization with approximate gradient. SIAM Journal on Optimization, 19(3):1171- 1183, 2008
Lin, J. A., Padhy, S., Antoran, J., Tripp, A., Terenin, A., Szepesvari, C., Hernández-LLoaall J. M., and Janz, D. Stochastic gradient descent for gaussian processes done right. In The Twelfth International Conference on Learning Representations, 2024. URL https: openreview net forum? (aaemmnddddfffffemm.m
Deisenroth, M. P., Fox, D., and Rasmussen, C. E. Gaussian processes for data-efficient learning in roboties and con-
10

---

Error Bounds for Structured Kernel Interpolation
Liu, G. and Onnela, J.-P. Bidirectional imputation of spatial GPS trajectories with missingness using sparse online Gaussian process. Journal of the American Medical In- formatics Association, 28(8):1777-1784 2021.
Yadav, M., Pleiss, G., Gardner, J., Weinberger, K. Q.. and Wilson, A. G. Faster kernel interpolation for Gaussian processes. In Proceedings of the 38th International Con- ference on Machine Learning, pp. 11279-112211 1228.. PMLR, 2021.
Modell, A. Entrywise error bounds for low-rank approxi- mations of kernel matrices. In The Thirty-eighth Annual Conference on Neural Information Processing Systems.
Yaday, M., Sheldon, D., and Musco, C. Kernel interpolation with sparse grids. In Advances in Neural Information Processing Systems, 2022.
Moreno, A., Mei, J., and Walters, SKI to go faster: Accelerating toeplitz neural networks via asymmetric kernels. arXiv preprint arXiv:2305.00 000028, 2023.
Morley, J., Machado, C. C., Burr, C., Cowls, J., Joshi, I., Taddeo, M., and Floridi, L. The ethics of Al in health care: mapping review. Social Science & Medicine, 260: 113172,
Musco, C. and Musco, C. Recursive sampling for the Nys- trom method. Advances in neural information processing systems, 30, 2017.
Petersen, K. B., Pedersen, M. S., al. The matrix cookbook.
Technical University of Denmark, 7(15):510, 2008
Rasmussen, C. and Williams, C. K. Gaussian processes for machine learning. MIT press, 2006.
Snelson, E. and Ghahramani, Z. Sparse Gaussian processes using pseudo-inputs. In Advances in neural information processing systems, pp. 1257-1264, 2006.
Stoer, J., Bulirsch, R., Bartels, R., Gautschi, W., and Witz- gall, C. Introduction to numerical analysis, volume 1993.
Springer, 1980.
Stonyakin, F., Kuruzov, I., and Polyak, B. Stopping rules for gradient methods for non-convex problems with additive noise in gradient. Journal of Optimization Theory and Applications, 198(2):531-551, 2023.
Titsias, M. K. Variational model selection for sparse Gaus- sian process regression. Report, University of Manchester, UK, 2009.
Wild, V., Kanagawa, M., and Sejdinovic, D. Connec- tions and equivalences between the nyström method and sparse variational Gaussian processes. arXiv preprint arXiv:2106.01 121, 2021.
Wilson, A. and Nickisch, H. Kernel interpolation for scal- able structured Gaussian processes (KISS-GP). In Inter- national conference on machine learning, pp. 1775-1784. PMLR, 2015.
Wynne, G. and Wild, V. Variational Gaussian processes: A functional analysis view. In International Conference on Artificial Intelligence and Statistics, pp. 4955 4971. PMLR, 2022.

---

Error Bounds for Structured Kernel Interpolation
A. Auxiliary Technical Results
Lemma A.1. Given a function J IRd - R of the form [ 22, ... ma) I j=1 Ji(tz), where each Jj B - R. Let
G G ( x G(d) be fixed d- d-dimensional grid, where each G {pf () p2 ..... ...,Pi;; is finite set of n
(3)
grid points along the j-th dimension for j = 1,2,.... d. Then the following equality holds:
d S
ΣΣ k==1j=1
I
Σ 50lk)
Proof. This is essentially repeated application of the distributive property.
By Induction on d (the number of dimensions):
Base Case (d = 1): When d = 1, the statement becomes:
11
1 ( = Σ
===1
This is trivially true.
Inductive Hypothesis:
Assume the statement holds for d m, i.e.
m1
122
Pmm
222
172
Σ
II
Σ 1
k1=1k== 11
=m==1==1
j=1
Inductive Step:
We need to show that the statement holds for d = m + 1. Consider the left-hand side for d = m + 1:
n2
Σ
Σ ...
Σ
II
1)(P))) (1)
km+1=1 j=1
We can rewrite this as:
fg
Aon 11
Σ
Σ
II (
Sm+ (olm+t))
km+1=1
Notice that the inner sum (over 1m4 does not depend ki, h2 ... ki. Thus, for any fixed values of k1, ka, kmm
can treat I f3 (t)) as a constant. Let C(ki,..., km)) = 11 (ok)) Then we have:
C(ki,
kim)
Σ
Em+i(pemtt (mn+1)
(Pkmtt
R1 =1 1==1
Rm =1
Rm+1=1
Now, the inner sum have:
Σ Rm
=mmii=1 fm+i (
(m++1)
is a constant with respect to ki...., km. Let's call this constant Sm+1. So we
12

---

Error Bounds for Structured Kernel Interpolation
Σ
(mm Sm+1
(
k-= k1
kg=1
Km
=1
By the inductive hypothesis, we can replace the nested sums with product:
Sm+1
II
Σ k==1
Σ
fm+
(Pim
(m+
II ==1
Σ Ry=1
j=l
Amm+1= =1
Rearranging the terms, get:
Am++
m++1
II Σ fiP
Σ
fm+1 1 (Dmmtt))
= II j=l
Σ ty
1=111=
k==
This is the right-han gentthhnddai side of the statement for m Thus, the statement holds for d + 1.
Conclusion:
By induction, the statement holds for all d 2 1. Therefore,
d
n Σ
= II ==1 R1=1
12=1
ka=1j=1
Claim Given convex combination C C= aA (1 B, where [0, 1). and A and B symmetric matrices, the eigenvalues of C lie in the interval [min (An (A), A, (B)) max (A) A1 (B))].
Proof. First, recall that for symmetric matrix A, the Rayleigh quotient R(A,x) x Ax x
is bounded by the smallest and
largest eigenvalues of A:
AA S &
Consider the Rayleigh quotient for the matrix C:
R(C,x) = (aA B)x aR(A (1-
x x
Since R(A, x) and R(B, x) are bounded by their respective eigenvalues, we have:
R(C,x) S aAi (B))))(BB
which implies:
R(C,,) S max (aax(11(a) A
Similarly,
R(C, x) 2 (Ai A, (B))
Thus, the eigenvalues of C aA (1 -a) bounded by:
min(A, A, (B) A(CC)= < (A),
13

---

Error Bounds for Structured Kernel Interpolation
B. Proofs Related to Important Quantities
B.1. Proofs Related to Ski Kernel Error Bounds
B.I.1. PROOF OF LEMMA 4.1
Lemma 4.1. The uniform error bound over a compact domain S Rd for tensor-product cubic convolutional interpolation of thrice continuously differentiable function f RRd -> by its interpolant g TR - R
sup (f() g(x) 0 =Ex
or equivalently O
od m372
Proof. We define a sequence intermediate interpolation functions. Let (o) f(x) be the original function. For
1,..., d, we recursively define (x) as the function obtained by interpolating 9:-1 along the i-th dimension using the cubic convolution kernel
Iz (cx); kh
(xxx) Σ Ji-1 (x ((cx); 24 k,) u
h
Here, Cx is the grid point closest to x, and e; is the i-th standard basis vector. Thus, 91 (x) interpolates f along the first dimension, (x) interpolates 9I along the second dimension, and so on, until ga(x) gxx) is the final tensor-product interpolated function.
We analyze the eиror accumulation across multiple dimensions using induction. Using (Keys, 1981), the error introduced by interpolating thrice continuous differentiable function along single dimension with the cubic convolution kernel is uniformly bounded over the interval domain by K Kh3 for some constant K > 0, provided the grid spacing is sufficiently small. This gives us the base case:
191(x) 90(()] Kh3
For the inductive step, assume that for some = the error is uniformly bounded by
[g Dk (x-1(x)]== S + kh.
14

---

Error Bounds for Structured Kernel Interpolation
We want show that this bound also holds for i k + 1. We can express the difference 9h+1 (x) Ih (x) follows:
yh (x + ((xx)x+1 16+1 + 1 [R+1 (cx)x+1 ki++ 1h
9k- 11-1(x)
-x(xx=
Σ
h
1
9k (x)
Σ
Σ 9k-1 (x + ((xx)) - Js + kehjoz + ((xx))+1 - Jx+1 + kiti
1-+1=-1
(cx)* kuhh ]
-
xR+1 (cx) ki+1h
u
h
[kk (cx)x kyh
Σ -k-1 (x+ + (((x) 26 + kxh))z)
kk
ck (cx)* - kahh
h
--=-1
x
9h-1 (x + ((cx)
Is + kxhjey + ((cx)++1 - 0641 + hi
h141=-1
Tk+1
(cx)x=11 ki--th) gk-1 +x+ ((cx) Is + kkh)oz)
h
Denote by 4 the one mim-emmesioonal cubic interpolation operator acting in coordinate i:
(7,h)(x) Σ 1(x ((cx): 24 + ee)) u [i (cx): th
With this notation 9i Lig-1 and 90 =
The inner term in the last expression is the interpolation error of 9k-1 along the (k+1)-st variable, evaluated at Yk. + ((cx)::
kxh)ez:
2
Ik- (yka + ((cx)++1 Pkt1 + (e+th))o+t) u Tk-1 (cx)x+1 kx+2/1) -k- (ykk)) ((I++I-Id)ge 1 (y
Σ
#441=-1
Hence,
Ik (cx)x kyh)
(x)+ 9x) 5
Σ
sup ((4h-) - i))
u
By the definition of e (uniform & bound on the weights), the sum of absolute weights S c, so
(x)+1()) (x)) s e. [[(T+1 Td)gk-1llss
We now apply the 1-D cubic interpolation error in the (k+1)-st variable to gh-1- Since previous interpolations act on coordinates 1, 1,...,k only, 9k-1 remains 3 in Ø441 with
||03 182_,9h-111 9h-1 In...---1002 (12,, A S 11o1 ... .......---11 11 110 oo ck = S f|
Thus there exists a constant K > (depending on ||03__ foloo but not on h) such that
[1(T++1 Id)gk-1lle S K x
15

---

Error Bounds for Structured Kernel Interpolation
and therefore
9k+1(x) - gx(x)) & c. =
This completes the inductive step.
Finally, we bound the total error [g(x) f(x) = -lga(x)= gu(x)) by summing the errors introduced each interpolation step:
d
d-1 K c 1=0
1g(x) - f(x)] - Σ 19:(x) 91 1x-1 - Σ c- Kh
The last sum is a geometric series, which evaluates to Kh Kha1- For a fixed c (independent of d), this expression (Oe)) when d is large. Therefore, tensor-prose productt cubic convolutional interpolation has te) error. Finally, noticing that h O mi/a gives us the desired result.
1- I-c"
B. 1.2. CURSE OF DIMENSIONALITY FOR KERNEL REGRESSION
The next lemma shows that when using product kernel for d-dimensional kernel regression (where cubic convolutional interpolation is a special case), the sum of weights suffers from the curse of dimensionality. The proof strategy involves expressing the multi-dimensional sum as a product sums over each individual dimension, leveraging the initial condition on the mie-dimensional bound for each dimension, and taking advantage of the structure of the Cartesian grid.
Lemma B.1. Let u R + be a one one-dimensoonal kernel function with constant be a d-dimensional product kernel defined as:
0 defined as 3.3. Let ua : R + R
d
x
1
Tj
x
Ud
h
h
==1
where x = (xx(1), _2(21 ,.., (()) &Rd and x;
(x. 3 x E R&d are d-dimensional points. Assume the data points
(1)
(2)
{xi)==1 (n may differ from the univariate case) lie on a fixed d-dimensional grid G G x G( x G(d), where each GU3 {p(") (5) (j)
P2 : ..... Pa; } is finite set of ny grid points along the j-th dimension for j = 1, ..., d. Then, for any x E IRd, the sum of weights in the d-dimensional kernel regression is bounded by cd:
(1)
a
x
Di
S
Ud
Proof. Let the fixed d-dimensional grid be defined by the Cartesian product of d sets of 1-dimensional grid points: (j)
G= G() ( x G(d), where GO) {p("), .p......ppp ....,pp?? is the set of grid points along the j-th dimension.
We start with the of weights in the d-dimensional case:
1
2
z(
ΣI 1=
Ud
u
h
h
Since the data points lie on the fixed grid G, we can rewrite the outer sum as a nested sum over the grid points in each dimension:
10)) - x u
d
((0)
Pki
(1)
=
"
h
h
12=1
ka=1j=1
Now we can change the order of summation and product, as proven Lemma A.1:
16

---

Error Bounds for Structured Kernel Interpolation
d
200 -
Pkj (3)
n2
x10)
(3) Pk;
ΣI u 1==1k==1 k==1z=1
II
Σ
=
"
h
h
By the assumption of the lemma, we know that for each dimension j. the sum of weights bounded by c. Note that (n ny kk== is simply set of points in R, thus:
x((i)
(
Pks
S
Therefore, we have:
x1133
())
)
Pky
Σ
ca
h
k
Thus, we have shown that:
x
S se
h
1=1
1.1.3. PROOF OF LEMMA 4.2
Lemma 4.2. Let 8, denote the uniform error bound of tensor-product cubic convolutional interpolation (Lemma 4. f. m inducing points and interpolation degree L,1d (hence 1,1/d = for cubic interpolation, as used in SKI). The SKI
kernel k; X x X + IR approximating a thrice continuously differentiable (not necessarily PSD) kernel k X X + R with uniform grid spacing in each dimension has error
k(x,x') - x((xxx)) S 1, + Adm.L
cdd
0
m33d
Proof. Recall that SKI approximates the kernel as
k(x, x') ~ f(xxxxxx x') = w (x)((x)),
Let KU,'' n be the vector of kernels between the inducing points and the vector x'
1k x') - k((x x')l = |k(xx x') - wx(x) KU,x' + wx)) KUxx' - w(xx Kuw((x))
1x(xx x') w(x) Ku,x + /w(x) Kuz"' w(x) krw(('))
S 8,, + |x(x) T KU, un(x) T Kuw(x')]
since |k(x, x') wx Kux is single polynomial interpolation
(6)
Now note that w(x) e Rm a sparse matrix with at most non-zero entries. Thus, letting w(x) E R be the non-zero entries of w(x) and similarly Kux R4 be the entries of Ku, in the dimensions corresponding to non-zero entries of
17

---

Error Bounds for Structured Kernel Interpolation
w(x) E IRm while Ky LR is the analogous matrix for Ku, we have
[w(x) Ku, u0(x) (x)((x))) = ww(x) T KUX w(x) Kuw(x')
I (x)(x) KU, - Kuwx(x')ll
ea) Kux Kvw((x)) (oo Lemma B.1
cd 8
where the last line follows as each element of (x)((x') a polynomial interpolation approximating each element of Kux* Plugging Eqn. into Eqn. 6 gives us the desired initial result of
kxx, x') - k(x, x')l S ōm, + c 6m.
and Lemma 4.1 gives us the result when the convolutional kernel is cubic.
1.4. PROOF OF PROPOSITION 4.3
Proposition 4.3. For the SKI approximation K of the true Gram matrix K, we have
| K||2 S n (5m,L + aat (mmLL
Mn.mI.
ne2d
m33/d
Proof. Recall that for any matrix |A||2 S || Since K is symmetric, we have
K S K K K
K
K oc
Furthermore, K K] is the maximum absolute row sum K K.
Since there are rows and, by Lemma 4.2, each
element of K K is bounded by 6, + d em5m, in absolute value, have
KK K I n (5m,1 (mme) -n
Therefore,
K
K|12 Ln
B. PROOF OF LEMMA 4.4
Lemma 4.4. Let X Rxx" be the matrix of kernel evaluations between T test points and training points, and let K. x E Txxn be the corresponding SKI approximation. Then
max(nn 2d
K K. x =
m3/d
Proof. Using the same reasoning as in Proposition 4.3, we have
K x K xl122 S V /IIK K. xlll K. x.llla xl1oc
I max (||K,x K, | K xl
Now, K. K. xl11 is the maximum absolute column sum, which is less than or equal to T(6m + (mmmm) = mm Similarly, K K xlloo the maximum absolute row sum, which is upper bounded by n n65, + Amm,))) = -n.m.
Therefore,
K KxX12 S max(1T, I.
18

---

Error Bounds for Structured Kernel Interpolation
B.1.6. ADDITIONAL SPECTRAL NORM BOUNDS
Lemma B.2. Let x 1T be cross kernel matrix between T test points and n. training points, where the SKI
approximation uses inducing points. If the kernel function & is bounded such that [l<(xxx²)] M for all x,x' E X, then:
K x1122 S VTnM
Proof.
K x/12 S V X
VTnM
Lemma B.3. Let K. E RTxx be the matrix of SKI kernel evaluations between T test points and training points, where the SKI approximation uses inducing points. Let ((.) E TRR and W(X) Raxm be the mutrices of interpolation weights for the test points and training points, respectively. Assume that the interpolation scheme is such that the sum of the absolute values of the interpolation weights for any point is bounded by ch, where c > is a constant. Let Ku E Rmxm be the kernel matrix evaluated at the inducing points. If the kernel function & is bounded such that [x(xxx')| S M for all x,x' E X, then:
K. x1l22 S VTnme M
Proof. By the definition of the SKI approximation and the submultiplicativityy of the spectral norm, we have:
K x212 | (())k s |/llll]]K W(X)
We now bound each term.
1. Bounding [/(2)]l2 and [|W(X) |l: Since the spectral norm is induced by the Euclidean norm, and using the assumption that the sum of absolute values of interpolation weights for any point is bounded by cd, we have
II(-))22 < VIWOlIOWWL s VTd. d VTc.
Similarly, [|(X)] S Vic"
2. Bounding Ku 1l: Since Ky is symmetric, ||Kull of k),, and each row has m entries, so |Kull> S mM. Thus, |uulla S M.M.
Kullcc Each entry of Ky is bounded by M (by the boundedness
Combining these bounds, we get:
K x S (MM))) 2d M
as required.
Lemma B.4. Let K be the SKI approximation of the kernel matrix K, and o² > 0 be the observation variance. Assume that K, K are PSD. The spectral error of the regularized inverse can be bounded as follows:
(kk (K+O')) (K+0' 0²1) o2) -1
Tn,m
<
Proof. Note that
(K +ko0)) (K
[K + +k'0 IK [K- K)) (K + -21)
19

---

Error Bounds for Structured Kernel Interpolation
Taking the spectral norm, we have
|| K + +to1)
(K (K+
1 K + of
K k (K
(K+ + -
by Proposition 4.3 (K (K+02)) -111 (K + 022) - -
S m.
In,ma
B.2. Proofs Related to Linear Time Analysis
B.2.1. PROOF OF THEOREM 4.5
Theorem 4.5. If the domain is [- DJd, then to achieve a spectral norm error of | Kl S is sufficient choose the number of inducing points such that:
Th (1) 4/3
M
for some constant K' that depends only on the kernel function and the interpolation scheme.
Proof. We want to choose m such that the spectral norm error
K Kl12 S. From Proposition 4.3, we have the bound:
K
K| S n(1 + e² )m)LL
For cubic interpolation, from Lemma 4.1, we have the scaling for the elementwise error:
8, K'
where K' a constant that depends only on the kernel function and the interpolation scheme. Therefore, a sufficient
condition to ensure K
K
|| S
is:
n(1 + KK Se
(8)
The inducing points are placed on a regular grid with spacing in each dimension over the domain [-D, DJd. The number of inducing points m is related to h by:
h 2D m 1/d
Substituting this into the sufficient condition (8), we get:
n(1 + K'
2D mll/
3
S e
Rearranging to isolate m, we obtain:
m3/d (1 +
( + cd K' (8) )3) 4/3
(-1
20

---

Error Bounds for Structured Kernel Interpolation
B.2.2. PROOF OF COROLLARY 4.6 Corollary 4.6. If
(1 (1+c n(logn)) 3)/d
E2
(4)
C3/d
3/d
for some constants K' C > O that depend on the kernel function and the interpolation scheme and we choose m > based on the previous theorem, then we have both |
K1[2 S t and SKI computational complexity of (n))
Proof. Assume that
(1 + 2c D'
n(log n)3/dd n//d
C3/d
Rearranging this obtain
"h + K
d/3
C
log nl
0 n
log
Now taking
( (2)) d/3 m
n
we have that
0
and by Theorem 4.5, K
K||2 S
Now plugging in
log into log we obtain
log
n1
n
O (m log m) =
log
logn
log m1
n
n
0
log
log log
log n
log
(n)
desired.
C. Proofs Related to Gaussian Process Applications
C.1. Proofs Related to Hyperparameter Estimation
C.I.1. PROOF OF LEMMA 5.3
Lemma 5.3. [Bound on Derivative of SKI Kernel Error using Kernel Property of Derivative] Let Wlg (2. x') be the SKI
approximation of Kg, x.), using the same inducing points and interpolation scheme ko- Then, for all x. x! E X and all O E O the following inequality holds:
Öke (z, x²)
(ko x²)
Ke, (xx Ro a0),
09:
00,
$ 8m,, + c. c.
2d
m//
where &m.L is an upper bound on the error of the SKI approximation of the kernel K/, (x, x') with m inducing points and interpolation degree 1.2/d_ 1, as defined in Lemma 4.2.
21

---

Error Bounds for Structured Kernel Interpolation
Proof. By assumption, K/ x²) (ha(a.a)) is a symmetric kernel. The SKI approximation of Ro, (ar, a²) using the same 00;
inducing points and interpolation scheme 2x", is given by (1r, x'). For the kernel kg, (1r, ('),)),, have:
Kg, (a, x2) Rig (x,))")) S 8,,
where m, is the upper bound on the error of the SKI approximation of Ko, (a, x') as defined in Lemma 4.2.
Now, we need to show that 06a((x.x) (1, x²). Recall that the SKI approximation kg (xx, (") is a linear combination of 3e, kernel evaluations inducing points, with weights that depend on and r':
6o x²)
W W x' (o u)
where Wz((1, x²) are the interpolation weights. Taking the partial derivative with respect to 8,, we get:
Oko x¹)
(kg w)
W x²)
00,
00,
j=1
172
Wy((1, xf Dklg (ug. un)
This is precisely the SKI approximation of the kernel k%g, (x,,x") using the same inducing points and weights:
(x( = Σ>
Wala x 3ko, (uj. us).
1m1
Therefore, Okulz.a)) (x, 09,
Substituting this into our inequality, we get:
Okip Ōkg x²)
(ko(x,, x¹)
N/. (xx, Ko x²))
09,
00:
S 8.,,1. + I'
C.1.2. PROOF OF LEMMA 5.4
Lemma 5.4. [Partial Derivative Gram Matrix Difference Bound] For any I e {L,...,p]
OK
Ok
00,
00,
2
2d
0
33dd
where LL,, is the bound on the spectral norm difference between the kernel matrices corresponding to Ko, and its SKI approximation k/, (analogous to Proposition 4.3, but for the kernel W/
22

---

Error Bounds for Structured Kernel Interpolation
Proof. Let 2 be the kernel matrix corresponding to the kernel Ko,(x, x') corresponding to its SKI approximation k/ (2,,
(he(m..')
and let K2 be the kernel matrix
00,
From Lemma 5.3, we have:
Oko(x,x))
ko, (2,
(9)
08:
Therefore:
0K
0K
Ko.,
(10)
00:
00,
By Proposition 4.3, have bound on the spectral norm difference between kernel matrix and its SKI approximation. Let Net% N...,. be the corresponding bound for the kernel Ka,, and its SKI approximation kg Then:
[ - Ko,ll2 S Inm.
(11)
Thus,
OK 0k 00,
=|[K,, - K2,el22 S m.
This completes the proof.
C.1.3. PROOF OF LEMMA 5.5
Lemma 5.5. [Score Function Bound (Final Version)] Let V (0C and V vo((0) be the true and SKI log- log-likelihood gradients.
Assume the kernel sufficiently regular such that its partial derivatives are bounded by ||
0K
06, ll2 S C, for constant
C. = (n...
Then the spectral norm of the score function error bounded as:
(|Vo(()) - vož(0)12= (0 VP
m3/d (llvli2n2 + n²))
Furthermore, under the common assumption that =1 = (n)) the 2/2 norrmllized error is bounded as:
-Z -0 ( VD-VV--V vp-"² n²
m3/d
Proof. We start with the expressions for the gradients:
VC(0) v ( 2
(-22 (K -
log K + o²Il log(2rr)
"
2
Vε(0) T (k -
log KK + 0²1|
log log(2rr)
2
Thus, the difference is:
23

---

Error Bounds for Structured Kernel Interpolation
1
1
V
2 (K + +
2 log K + 0411
2 (K + -2-1--yy log K + +021
2
V
28 +k (K + 0 +
y
v
log KK 0²1| - log K + 2
T2
We will bound Ti and T2 separately.
Bounding Ti:
1
Vo (xT
Ti 2
1
30, yT
(K -K+
2
1
a (20 (IK+ ( (K+ + 02I)-1) -1) y)
P max
V
2
ISISP
yT (IK+ -(K
(K + o²I)-11
max 1s
2
00,
We will then bound 2,y T (IK + ----0-1----0 (K +22)))))
y Using the following equality 70, we can express this derivative as quadratic form as a difference between two quadratic forms and apply standard techniques
8
- (
24

---

Error Bounds for Structured Kernel Interpolation
for bounding differences between quadratic forms.
θ 30,
+ +2+1--1-11- (K - + oa1) -01)-) y)
y
(IK -21) 1 (K + ofo)) -11
$ 1
CS inequality
00,
1
(K + d21)-1 (K + +2²1 -1)
08,
(00
0 K
(00,k)
a K
= 1
(K + oorr--
(K 1k+o 1-+0o1)---
- + (K+ + +o-1)-1
(K (K
00,
000
à K
i K (K + -02))
= 11x11
- -K + o 2) -1
K
+
00,
001
00,
+(K +021)-1
K(K + +0²1)-1
00,
2
0
= =11v112
(K + 1T -1
K
K
(K + 0²1)
00,
00,
-KK -KKoo*1)- -1
0 K
(K
+
1 + (K + 0 (
3 20,
[ (K + -
k)
001
< ||
[K -
K
K (K + -
000
00:
(a)
(0. k) (K (K+
(KK + - - (K + 1
00,
(K< (bx
0 K)
+ K + -2))
(K + -
08,
We now explicitly bound (a), (b), and (c).
11.
θ K
0
(a) s
(K + 22)
K
(K
+ o 21
00,
00,
2
K + 0²T)
-T-11112
0 K
0
K
00
00,
12
K
K
04
00,
00,
(Using Lemma 5.4)
(b) < |K + +o21)==00000-0 o2 - (K + 021)-1112
K
|[(K +
00,
12
< 1
Ô
1 (K + +0 - (K + -000)~1111 112
K
a0
7nmmL g6
0 K
(Using Lemma B.4)
00
2
25

---

Error Bounds for Structured Kernel Interpolation
where we able apply the last line since K, K are GP kernels and thus PSD. Then
(b) < Tn.m.L
i
K
o6
00,
1
< Cn mmm.
by assumption
06
and finally
a
(c) - (K + o21)
K
1(K + 0-1)-1 -or---K--00 + -
00:
1
1 (K + 0 0²1)-1
0
(K + 0²1) 1
K
30
mn.mI
0 K
(Using Lemma B.4)
σ6
00:
M..m.L σ6
Cn
Combining these, we obtain
Ti < 20
1
|1yl1?
vD
max IstSp
(hhme o² mmmmm.Cn
2
Bounding the Trace Term Error (T2):
To bound the error in the trace term of the gradient, we begin with the standard identity for the gradient of log-determinant with respect to parameter 0; (see e.g. (Petersen et al., 2008)):
i a0,
log K(0)) Tr K(0)-1 (K)
09,
The Tr v log JK + - log K + o²1l
represents the L2 norm of the gradient error vector for the log-determinant. We first analyze the magnitude of its 1-th component, which denote (T2)r:
(T2)
2/2 Tr (K (K
021+
oK 00,
Tr
(K OK)) 00:
OK
By applying the triangle inequality to the trace, can decompose the error:
((T)): S 2 Tr ((K+ o'1) ---K-------00 (K+ -K+021)
0K 00,
1
Tr
(K (K+
0K aK
2
00,
aA
We bound each of the two resulting parts using the property [T((A)] I n A, (K + I))
A ||2
B|22-
Let A,
(K o²I) and
1. The first part is bounded using ||- - A= < In from Lemma B.4:
Tr (A A
oK
< n2 A A
UK
< n 0
Vn,m. Cn
a0,
ae,
2
2. For the second part, since K positive semi-definite by construction, we have
A
1|2 S 1/02.
This yields:
Tr A 0k
ok
OK
ak
1
< nl/
A
112
< n o2 In,m,L.
00,
ae,
09,
a0,
26

---

Error Bounds for Structured Kernel Interpolation
Combining these gives the bound for single component: [(72))l O n 7m.mm.L
an am..m.
+
Combining the Bounds:
The total error is the L2 norm of the error vector, which can bounded by vp times maximum component- error.
(V voC(O)]22 & vD max
( ((T2)))
IsiSP
We combine our bound for (T2); with the bound for the quadratic term error, ((T))) =0 to arrive at the final expression:
[1y|
(1im,, & + 22 (a.
o²
1y11%%
n² To
n n
VoC(O)- -
04
Assuming |1yl2 Onn) and substituting the known bounds for y and y', the dominant terms yield the simplified scaling relationship:
Vvc( -000l (224VP n
m33/d
C.2. Proofs Related to Posterior Inference
C.2.1. PROOF OF LEMMA 5.8
Lemma 5.8. (SKI Posterior Mean Error) Let pu()) be the GP posterior mean at set of test points E T and ii(:) be the SKI posterior mean those points. Then the SKI posterior mean 12 error bounded by:
[p()) m(-112 2d
om Vno max((T + VTnMn
Proof. We start by expressing the difference between the true and SKI posterior means:
K. x (K (K+ y k. (K (k+ + +
= (K K,)) (K (K+o') -1 y k K. [ (K+O')) - (K -(K+o-)) x//
27

---

Error Bounds for Structured Kernel Interpolation
Applying the triangle inequality and submultiplicative property gives:
x² | K - K. xxl12 12 +
[ + 02) - + 02T)
-1
||y/|
1
max
m,, )
1b/12 + || ||
(K< + o21)
- -KK+ -22) -1
| Lemma 4.4
o²
max mL LL)
1yl12 + VTnM
(K
(K -K -1
||ll| Lemma B.2
max
mmT I...m.
1yl12 +
vTnM
fmm,,1lyll2 Lemma B.4
o²
0²
max (MT, m.), + VTn
11yll
02
In,m,
2dd
11y122 max + VTnMn o²
= o2m3/a llylall
C.2.2. PROOF OF LEMMA 5.9
Proof.
First, note that
|Σ()) - 2()1k2 S |K..- - K
+ K x(K + o21)- - K x(K + 1211)11ks 1Kx 1|2
S m. + || (K + -221)-2-- - x + 0²1) --- _Kx 11,
where we used Proposition 4.3 and the fact that JK,- - KK,,Il S I. Now, we bound the second term, which is a different between two quadratic forms:
K (K + 12 K x ok Kx,
< K [K + +221)= -Kx, Kxx K + +021) - Kx,112
+ KK ((K + +02) -xx. - K + +021)= Kx, ll2
K. x (K. + 0- 22
(x,)]l2 + ||K x(K 0²T) x -02)) |xx, x 1|2
K. [K + 1011~~111 |x. - -Kx.112 + | (K + - x + 0 1) - Kx.
Xl Kx.112 + + -
K.
+ +0²1) 1|2 ||2||| Kx, 11.
where we used the fact that (K + (K+021)- 1 I.
Next, we bound the term
K xK + 0²1)-1 1 x + x221
K x((K + 22) - K, (K + 00 -01)
(K + 0 12 1 K x (K + =2
K. x KK + -21)-1 x 021)-111 1/2
+
< (K + +0) 1 - K. K ( + 021)
+ K. x (K- + +021)=1 - ((K + 0² -211-1112
K [ + 02 (K + o2 -11 1112
K x KK +0²1)
K x1xx1211( (K + 0212 -1 (K + +022) 2|2 + K. x Kx11211(K (K + +
K x||2
In
+ K kx 2
28

---

Error Bounds for Structured Kernel Interpolation
where we used Lemma B.4 in the last inequality. Substituting this back into the main inequality, we get:
x 1P1))KKx K. x(K ( 1xx
xxx 1x. 1/2 +
||, m.mm + + K. x112
1) ||xx, ll2
xllall - -xx, 1k,1112- +
In.ma
||K x [x,,xxx 1 oz K K /alll//Kxx 112.
σ1
Using Lemma 4.4 and the fact that |[x,, - Rx,112 S max( 7.m,..)) and that K,x = KX we have ||| = |Kx,-ll2.
112. Also, by assumption, |Kxx,-1a VTnM. Using Lemma B.3, have ||Kx, 2l S VTnme M. Substituting these bounds, get:
x O2 Kx - KK + 0022)=
VTnM
max Imm, +
nnmmL
(VTTM) (V Tnme M) + max m, L. mn,,, ( ( M)
1
o2
04
o²
VTnM
VTnme M
max Lmm.
mnm Tnma M
max (T.x L)
0²
04
Finally, substituting this back into the original inequality, obtain the desired bound:
VTnM
E(- E(-1l2 I mI +
max M, mn.
o²
T M
"n.mm
Tnme 2d M
max In.m...)
04
o²
Tn-meddM2 M2
Vnma M max(T,n)
=0
m3/d
29


Reviews:
[{"review_id": "JOZpXs6wCN", "reviewer": "Reviewer_y1bt", "summary": "This paper presents a theoretical analysis for the structured kernel interpolation (SKI, Wilson and Nickisch, 2015), a popular method for scaling Gaussian processes (GPs). The paper proves the error bounds for the SKI gram matrix and examine its effect on GP hyperparameter estimation and posterior inference. Crucially, they identify the number of inducing points needed for error control is $O(n^d/3)$. When $d\\leq 3$, SKI can achieve linear complexity for any error tolerance, but for $d>3$ the error must increase with sample size to maintain linear time. This analysis provides key insights into SKI's accuracy-scalability trade-offs.\n\n## Updated after rebuttal:\nThe paper provides a novel theoretical contribution to SKI analysis. The biggest limitation of the current submission is the lack of practical implication. Furthermore, according to authors rebuttal, the experiment results suggest that the bound can be pessimistic and the practical performance of SKI seems to be better than what the worst-case theory indicates. Based on these points, I would suggest authors to incorporate a more comprehensive empirical validation in the next draft to bridge the theory and practice. I will adjust my score from 4 to 3.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your positive feedback and insightful suggestions! We are encouraged that you found the theoretical claims rigorous and the proofs correct.\n\nWe will address the clarity issues raised:\n* We will explicitly define $k=1,\\ldots,K$ as the iterations when we first use it in Section 5.1.\n* We will clarify the definition of interpolation error in Lemma 4.1 (as the uniform error bound over the domain using functions) and $\\delta_{m,L}$ in Lemma 4.2 (as the tensor-product convolutional cubic interpolation error). We will also ensure the specific norm (spectral, elementwise, L2) is clear whenever \"error\" is discussed.\n* We will correct the noted typos (e.g., \"healthcar\").\n* We confirm Lemma 4.2 should use an inequality ($\\le$), and will correct this.\n* Regarding Lemma A.1, you’re right that it essentially shows the distributive property; we included the inductive proof for formal completeness but will add a few sentences that summarize it before the full proof.\n* Regarding LLM usage: LLMs assisted with brainstorming, literature search, and initial proof drafting. However, all final mathematical arguments, proofs, and derivations presented were significantly refined, validated, and verified by the authors.\n\nRegarding the request for empirical validation: We appreciate the suggestion and performed experiments simulating the spectral norm error $||K - \\tilde{K}||_2$ for varying $n$ and $d$ using an RBF kernel on synthetic data. Specifically, we set the number of inducing points according to our theoretical scaling, $m = \\lceil n^{d/3} \\rceil$, for dimensions $d=2, 3, 4$. We then plotted the log of the approximation error against the log of the sample size $n$.\n\nOur theoretical bound (Prop 4.3 combined with $m \\propto n^{d/3}$) predicts that the error should asymptotically approach a constant $O(1)$ independent of $n$. However, our experiments showed a more favorable outcome for all tested dimensions ($d=2, 3, 4$), the approximation error actually decreased as the sample size $n$ increased with $m$ scaled precisely as $n^{d/3}$. This suggests our theoretical upper bound, particularly the $O(n)$ factor in the numerator derived from the infinity norm bound, is likely pessimistic, and the practical performance of SKI with this scaling is better than the worst-case guarantee. We will add these figures along with experimental details to the next revision.\n\nFurthermore, upon reviewing our theoretical derivation (specifically Corollary 4.7 and its proof), there is a typo: the correct regimes are $d<3$ and $d \\ge 3$, but this was incorrectly stated as $d \\le 3$ vs $d > 3$. We will correct this inconsistency throughout the text. Thank you again for your valuable feedback which has helped improve our paper.", "claims_and_evidence": "This paper remains theoretical. The theoretical claims are rigorous. However, I'd appreciate some empirical simulation to validate the theoretical rates, in particular for the regime of $d>3$.", "methods_and_evaluation_criteria": "The paper does not involve proposing a new method.", "theoretical_claims": "I have checked the proofs of main lemmas and theorems up to section 4 and they seem to be correct.", "experimental_designs_or_analyses": "NA", "relation_to_broader_scientific_literature": "This paper fills in the gap of theoretical analysis of SKI. SKI has been a popular scalable GP method, however, there has not been analysis in its approximation error and how that impacts the downstream GP learning and inference.\n\nThe error analysis builds on many existing theoretical tools, such as  convolutional cubic interpolation error analysis and gradient ascent analysis. But the application of these analysis to GP is still novel.", "essential_references_not_discussed": "NA", "other_strengths_and_weaknesses": "Strength: the paper is well-organized and for most part the intuition behind the theorems and proofs are clear. \n\nWeakness: \n\n- Writing clarity. I'd appreciate the authors make explicit references of the definitions for various symbols. For example, the K in section 5.1 (line 307 on page 6) is used without reference. It only became clear to me in Theorem 5.7 that $K$ stands for the number of iterations. \n\n- Technical clarity. A lot of technical terms are used without clarification/definition. For example, in lemma 4.1 how is the error of tensor-product cubic convolution interpolation defined? And in lemma 4.2, what is $delta_{m,L}$ exactly? I only got a sense of these terms after going over the proofs. Also please clarify what norm you used when mentioning \"error\". \n\n- Missing empirical validation. While I understand the paper's main contribution is the theory, I'd still appreciate some (toy) empirical study that confirms the theory, in particular in  the two regimes $d\\leq 3$ and $d>3$ and see how their behaviors diverge. And I think highlighting the empirical behavior when $d>$ can really help practitioners understand the practical limitation of SKI.", "other_comments_or_suggestions": "- Sec 4.1 title, Ski -> SKI\n\n- Section 7: typo, healthcar -> healthcare", "questions_for_authors": "- In Lemma 4.2, should we have inequality $\\leq$ instead of $=$ in the main equation?\n\n- App. A: is lemma A.1 essentially the distributive property of multiplication over addition?\n\n- The authors stated that the paper has largely used LLM. Can the authors clarify if any part of the proof is LLM-written? If so, which part, and if the authors have carefully examined through these parts?", "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": "", "extra_scores": {}}, {"review_id": "A6XWArV5Q9", "reviewer": "Reviewer_AYYe", "summary": "In this paper, the authors prove several bounds for quantities of interest when using the structured kernel interpolation approximation of Gaussian processes. These bounds include element-wise errors in the kernel matrix, the spectral norm of the difference between the approximated and true Gram matrices, the log marginal likelihood (referred to as the score function in the paper), and the $L^2$/spectral norm error between the posterior mean and covariance for the same set of hyperparameters. These bounds are derived using results from optimization theory and kernel ridge regression.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thanks for your helpful review and for recognizing the paper's comprehensive scope. We acknowledge the concern regarding the density of the paper due to proofs being in the appendix. In the next version, we will incorporate additional proof intuitions or summaries of the key steps for our main theorems and lemmas into the main body of the paper to improve readability, as requested. Specifically, while the main text already outlines the core strategy for some results (like Theorem 4.5 relating error bounds to grid spacing and $m$, and Lemma 5.6 using quadratic form bounding techniques for the score function), we will add or enhance the intuitions for several key steps:\n* For **Lemma 4.1** (Interpolation Error), we will note the proof uses induction on dimension $d$, showing how the error bound accumulates multiplicatively based on the 1D case from Keys (1981).\n* For **Prop 4.3** (Gram Matrix Error) and **Lemma 4.4** (Cross-Kernel Error), we will clarify that the proofs rely on standard matrix norm inequalities ($\\Vert A\\Vert_2 \\leq \\Vert A\\Vert_\\infty$ for Prop 4.3 using symmetry, and $\\Vert A\\Vert_2 \\le \\sqrt{\\Vert A\\Vert_1 \\Vert A\\Vert_\\infty}$ for Lemma 4.4) combined with the elementwise error bounds derived in Lemma 4.2.\n* For **Lemma 5.9** (Posterior Mean Error), which currently mentions the standard strategy, we will elaborate on the specific steps: using algebraic manipulation (add and subtract terms) and the triangle inequality to decompose the error $\\Vert \\mu - \\tilde{\\mu}\\Vert_2$ into terms involving the error in the cross-kernel matrix ($K_{\\cdot,X} - \\tilde{K}_{\\cdot,X}$) and the error in the inverse Gram matrix ($(K+\\sigma^2I)^{-1} - (\\tilde{K}+\\sigma^2I)^{-1}$), then applying our previously derived bounds for these components (Lemma 4.4 and Lemma B.4).\n* For **Lemma 5.10** (Posterior Covariance Error), which currently highlights the quadratic form bounding, we will add detail on bounding the prior covariance term ($K_{\\cdot,\\cdot} - \\tilde{K}_{\\cdot,\\cdot}$) using Prop 4.3, and explicitly mention how the update term's error decomposition utilizes the bounds derived in previous lemmas.\n\nWe will address the formatting and writing suggestions:\n* We will carefully revise referencing for consistency (including theorems in Sec 4) and ensure consistent capitalization (e.g., Gram matrix, SKI).\n* We will use textual citations (\\citet/\\textcite) where appropriate for better flow.\n* We will define the SPD (symmetric positive definite) acronym upon first use (line 125), but change it as you suggested to PSD for consistency with common GP literature.\n* We will review the older citations (Lagrange, 1795; Kolmogorov, 1940). These were included to credit the original foundational works. We will check for standard modern references.\n\nThank you for your constructive feedback, which will help us enhance the clarity and presentation of our work.", "claims_and_evidence": "The authors' claims are substantiated with formal proofs provided in the appendix.", "methods_and_evaluation_criteria": "There are no evaluation criteria.", "theoretical_claims": "I have checked some of the initial proofs and the ideas behind them, and they sound reasonable. However, I am not familiar with the literature on approximate gradient descent/ascent.\n\nOverall, I found the paper overly dense. Although earlier results are used in the later theorems, the paper focuses on theoretical analysis, yet all proofs are relegated to the appendix. I have set my score to *weak accept* due to this. I believe the paper would be more valuable if the authors included at least a high-level description of their proofs in the main paper, as these are the paper’s main contributions.", "experimental_designs_or_analyses": "There are no experiments.", "relation_to_broader_scientific_literature": "The contributions of the paper provide an important theoretical lens on the properties of SKI, whose experimental validations show that its strength lies in lower-dimensional data; this property could be related to the authors' bounds, which change asymptotic behavior when $d \\\\geq 3$. As far as I know, such an analysis has not been done before and represents a good initial step toward analyzing methods based on local kernel interpolation.", "essential_references_not_discussed": "I believe no essential references are missing.", "other_strengths_and_weaknesses": "A strength of this paper is that the authors explore all facets of Gaussian process approximation, from the kernel to the final predictive distribution.", "other_comments_or_suggestions": "- On line 254, the theorem in question was not referenced.  \n- Please revise the referencing of equations and sections, as the same paragraph (lines 294–300) contains inconsistent capitalization.  \n- \"Gram matrix\" should always be capitalized, as it is derived from the surname of Jørgen Pedersen Gram. Similarly, \"SKI\" should always be capitalized.  \n- Some citations would be better suited as textual citations rather than parenthetical citations. For example: *“In the second group, the foundational work by (Wilson & Nickisch, 2015) that [..]”* would be better displayed as *“In the second group, the foundational work by Wilson & Nickisch (2015) that [..]”*. You can achieve this by replacing `\\cite` with `\\citet` (if using natbib) or `\\textcite` (if using biblatex).  \n- On line 125, left column, the acronym SPD (semi-positive definite?) is not introduced. In GP literature, the acronym PSD (positive semi-definite) is more commonly used.  \n- Please avoid citing very old or inaccessible papers. For example, Lagrange interpolation is referenced as *“Lagrange, J.-L. Leçons élémentaires sur les mathématiques. Imprimerie de la République, 1795.”* However, it is unclear whether this citation format is correct, as I could not locate this version under this specific publisher. Additionally, it is not specified where in the text the interpolation is defined. The same issue applies to Kolmogorov’s 1940 paper, which is cited in relation to Gaussian processes.", "questions_for_authors": "No additional concerns.", "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": "", "extra_scores": {}}, {"review_id": "ncw0mzFmLO", "reviewer": "Reviewer_2Sc8", "summary": "For a fixed kernel, non-zero noise-variance parameter, and spatial dimension, the authors give an asymptotic bound on the number of inducing points needed to reach some accuracy threshold with the SKI approximation as a function of the number of data points.  Upper bounds on error in the kernel matrix, the log-likelihood, and the posterior mean and variance are all considered.  Based on their arguments, the authors recommend a scaling of inducing points like n^(d/3), at least when cubic interpolation is used for mapping onto the structured grid in SKI.\n\n## Update after rebuttal\n\nI do not believe that these bounds are tight; from discussion with the authors during the rebuttal period, it seems as though they are in agreement.  Nonetheless, after thinking about it, I think this is an interesting enough step in the right direction that I would be fine with seeing it in print.  I have thus updated my recommendation to a weak accept.", "questions": "", "limitations": "", "rating": 3, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your thorough and critical reading of our manuscript. We appreciate the detailed feedback, particularly regarding the assumption about kernel derivatives.\n\n**Regarding the critique of Assumption 5.3:** You are absolutely correct. Our claim that the partial derivative of the kernel $k_\\theta$ with respect to lengthscale $\\ell$ for RBF is a valid SPD kernel (Assumption 5.3 and a short description about RBF kernel’s before the formal statement) is false. For stationary kernels $k(x,y) = \\phi(||x-y||^2/\\ell^2)$, the derivative with respect to $\\ell$ is symmetric but indefinite. Thanks for noticing this. Crucially, we do not actually need positive definiteness. We only rely on the spectral norm error we derived: this starts with tensor-product convolutional cubic interpolation error (not for kernels specifically). It then moves to elementwise kernel error, which relies on neither positive definiteness nor symmetry. After that the spectral norm bound relies on symmetry (so that the $l_1$ and $l_\\infty$ norms are equal), but not positive definiteness. We will thus adjust the assumption to not claim positive definiteness.\n\n**Regarding Bound Tightness/Experiments:** We appreciate the request for empirical validation. We performed experiments measuring $||K - \\tilde{K}||_2$ on synthetic data using an RBF kernel. Specifically, we set the number of inducing points according to our theoretical scaling, $m = \\lceil n^{d/3} \\rceil$, for dimensions $d=2, 3, 4$. We then plotted the log of the approximation error against the log of the sample size $n$.\n\nOur theoretical bounds predict that the error should asymptotically approach a constant $O(1)$ for this scaling ($m \\propto n^{d/3}$). However, the experiments revealed a more favorable trend: for all tested dimensions ($d=2, 3, 4$), the approximation error actually decreased as $n$ increased. This suggests our theoretical upper bounds are likely pessimistic as you suggested and that the method's practical scaling performance is better than the worst-case guarantee.\n\nFurthermore, a careful re-reading of our theory (Theorem 4.5 proof) confirmed the correct threshold condition is $d<3$ vs $d \\ge 3$; we had made a typo and written $d \\le 3$ vs $d > 3$ and will correct this throughout the paper. We will include plots and details of these experiments in the next revision of the paper.\n\n**Regarding Novelty (Curse of Dimensionality):** While the curse of dimensionality impacting SKI is known, our contribution lies in its theoretical characterization that we can control error by growing inducing points at an $n^{d/3}$ rate and that for $d<3$ this guarantees linear time controlled error, further supported by our empirical findings showing performance exceeding the pessimistic bounds.\n\n**Regarding Proof Style/LLMs:** Lemma A.1 was included for completeness. Regarding LLMs, they assisted in outlining and literature search, but all final proofs were a combination of human and AI derived and human verified.\n\n**Regarding Generalization Error:** We agree kernel matrix error isn't the full story. Our bounds on the posterior mean and covariance approximation errors (Lemmas 5.9, 5.10) provide a more direct theoretical link between the SKI approximation and the resulting predictive distribution's accuracy.\n\nThank you again for your detailed critique. Clarifying the assumption issue and incorporating the empirical results will significantly strengthen the paper.", "claims_and_evidence": "The authors analyze error as a function of the number of inducing points.  This is not the same as the *time* to do almost any calculation involved.  Even a matrix-vector product with SKI involves some FFTs, adding a log term into the mix.  Moreover, the SKI construction just gives a matrix-vector product; computing likelihoods and linear system solves involves some additional iterative method.  The rates of convergence can indeed be bounded by the condition number, and this is controlled via the bound on the kernel (times the number of data points) divided by the noise variance.  But these are all potentially rather pessimistic bounds with very large constants.", "methods_and_evaluation_criteria": "This paper is purely theoretical.  No empirical evaluation is given, though I would have appreciated it in some places (it would be nice to have a sense of how tight the bounds are).", "theoretical_claims": "I checked the first proofs in the supplement, though I started skimming toward the end.  There are some erroneous statements, though I am not sure how significant they are.  For hyper-parameter estimation, the authors claim that dk/dtheta is itself a kernel; but this clearly cannot be true for any family k(x,y) = phi(norm(x-y), theta) s.t. phi(0,theta) is constant, since in that case the diagonal of the kernel is zero.  Hence, for example, this statement is invalid for most length-scale parameters for standard kernels -- the derivative matrices will be symmetric, but indefinite.", "experimental_designs_or_analyses": "There were no experiments.  I would have appreciated some, as the bounds given here are mostly upper bounds for which it is hard to judge tightness.", "relation_to_broader_scientific_literature": "The authors claim there is not much other error analysis in the literature of this type.", "essential_references_not_discussed": "I did not see any essential references not discussed.", "other_strengths_and_weaknesses": "Some of the proofs in the supplementary material were a bit of a slog.  The first proof, in particular, is essentially an application of distributivity; writing down an induction is fine, but seems like more than really needed.\n\nAs indicated above, I am not convinced of the tightness of most of these bounds.  Some computational experiments (or matching lower bounds, but that seems harder) would have assuaged my skepticism.\n\nThe fact that SKI suffers a curse of dimensionality is fairly well-known.  I also observe that the SKI kernel is a kernel, even if it is only an approximation to the original kernel.  However, the original kernel is usually chosen by a modeler based on a rough assessment of the regularity of the function at hand.  So altogether, it's not so straightforward to go from statements about error in the kernel matrix (for example) to statements about generalization error.\n\nThe comment about the use of LLMs in writing the model was interesting, though perhaps left me a little less trusting of the results than I otherwise would have been.", "other_comments_or_suggestions": "My comments all fit into the boxes above.", "questions_for_authors": "No important questions.", "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": "", "extra_scores": {}}, {"review_id": "8eDLTH7aYD", "reviewer": "Reviewer_vBju", "summary": "The authors provide a theoretical treatment of structured kernel interpolation in Gaussian processes, particularly focusing on cubic convolutional interpolation methods. They provide a comprehensive characterization of error and complexity. Their main results focus on the interplay between the number of inducing points and the error bound (both on the “elementwise” kernel values themselves and on the corresponding Gram matrix spectral norm), finding that there is a distinct difference in this relationship as the dimension increases. Namely, they show that any error tolerance can be achieved in linear time (given enough samples) for dimension less than or equal 3. Moreover, they bound the error of kernel hyper-parameter inference in SKI and posterior parameter approximation error.\n\n## update after rebuttal: I believe my initial score holds after the rebuttal, as the paper remains a strong contribution, in my opinion.", "questions": "", "limitations": "", "rating": 4, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you for your support and positive assessment of our work as a valuable contribution! We are glad that you found the proofs straightforward and the work well-positioned. Regarding your points:\n* **Clarifying the SKI Implementation (Sec 3.2 / 3.3 Connection):** We will make the connection between the general SKI formulation and cubic interpolation clearer. We will explicitly state that the weight vectors $w(x)$ in $\\tilde{k}(x,x^{\\prime})=w(x)^{\\top}K_{U}w(x^{\\prime})$ are constructed using the tensor-product cubic interpolation function $g(x)$ from Definition 3.2. In this context, $g(x)$ uses the cubic kernel $u(s)$ from Definition 3.1, and the function values being interpolated, $f(c_x+hk)$, corresponding to kernel evaluations between inducing points, $k(u_i, u_j)$, which form the matrix $K_U$. We will add cross-references to make this link explicit.\n* **Writing Suggestions:**\n    * We will define the SPD (symmetric positive definite) acronym on first use (line 125) but instead replace it with the suggested PSD (positive semi-definite) nomenclature for consistency.\n    * We will rephrase \"compute the action of the inverse\" to simply \"compute the inverse times a vector\" or similar.\n    * We will format the definition of the inducing points matrix U on its own line for clarity.\n    * We will ensure consistent capitalization of SKI.\n    * Regarding the comment on Line 234 (\"linearly for $d=2$\" should be \"linearly for $d=3$\"): Thank you for prompting us to double-check. You are correct and we will fix this. We appreciate the careful reading.\n    * We will fix the typo \"healthcar\".\n\nThank you again for your encouraging review and constructive suggestions.", "claims_and_evidence": "The main claim is related to the choice of the number of inducing points to achieve (sub-)linear time complexity with respect to a fixed error tolerance. This is shown rigorously in Section 4, where it is noted that this is possible for any error tolerance for dimension at most 3, and that the error tolerance must grow in higher-dimensional settings to satisfy linear time complexity. Other key claims are related to the hyper-parameter estimation via inexact gradient ascent on the approximate log-likelihood and approximation error of the posterior parameters. In all cases, the authors leverage common techniques in the literature to prove their findings.", "methods_and_evaluation_criteria": "As this is a theoretical work, there are no technical evaluation criteria. The authors employ common and new methods from the literature to show their results, such as recent work related to inexact gradient descent and common techniques related to error characterization in kernel regression.", "theoretical_claims": "I have checked through the proofs in Appendix B, and they seem relatively straightforward.", "experimental_designs_or_analyses": "Not relevant, as there are no experiments included in this theoretical work.", "relation_to_broader_scientific_literature": "The authors position their work in the context of three scientific domains. (1) Within the domain of theoretical analysis of Gaussian process regression/kernel regression with approximated kernels, they have contributed an analysis of a method which has previously been rarely treated: SKI. It is noted that prior works analyzing SKI focused only on one-dimensional features and did not analyze the posterior approximation. (2) In the context of the SKI literature, this works provides a theoretical basis for the widely-used method. (3) The authors also relate to theoretical works from which they utilize results, such as the inexact gradient descent literature.", "essential_references_not_discussed": "There are no missing references, to my knowledge.", "other_strengths_and_weaknesses": "This is a strong paper and provides valuable theoretical insight into a widely-used method. I believe this will be a useful contribution to the community. The only weaknesses would be a bit of lack of clarity in Section 3, as the exact implementation of the SKI method (particularly, the cubic convolutional interpolation) is not clear to me.", "other_comments_or_suggestions": "Left Column Line 125: I assume SPD kernel means positive semidefinite here and throughout the paper. I believe PSD is the more common nomenclature, and this acronym should be fully spelled out when first introduced.\n\nLeft Column Line 160-161: “compute the action of the inverse of the regularized…” I don’t think the word “action” needs to be included here. Is it not correct to just say “compute the inverse of the regularized…”\n\nRight Column Line 116-119: The inducing points matrix U definition should not be part of the paragraph, but should be separated into its own line.\n\nLeft Column Line 216: SKI should be capitalized\nRight Column Line 234: “linearly for d=2” should be “linearly for d=3”\n\nRight Column Line 242: healthcar -> healthcare", "questions_for_authors": "Right Column Line 125-128: It is not clear to me how the weight vectors w(x) can be obtained from the convolutional interpolation kernel u(s) of definition 3.1-3.2. Can this be made more clear, or at least can a more direct connection between Section 3.2 and 3.3 be made?", "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": "", "extra_scores": {}}]

Author rebuttal:

Meta-review: This manuscript develops error bounds that pertain to SKI, a popular method leveraged for efficient computation with Gaussian Processes. The reviewers agree that the theoretical contributions are mostly sound and sensible. Moreover, given the broad use of SKI there is value in the contribution. The reviewers do feel that there is room for improvement in aspects of the presentation. The main limitation of the manuscript in its present form is the lack of experiments to illustrate the bounds (and their looseness or tightness); the inclusion of experiments as discussed during the rebuttal will strengthen the manuscript (even if they show there is potential for further improvement of the bounds in certain situations).
Decision: Accept (poster)