High-Dimensional Gaussian Process Regression with Soft Kernel Interpolation
Authors: Chris L Camaño
Daniel Huang
Venue: Accepted by TMLR ()

Abstract: We introduce Soft Kernel Interpolation (SoftKI), a method that combines aspects of Structured Kernel Interpolation (SKI) and variational inducing point methods, to achieve scalable Gaussian Process (GP) regression on high-dimensional datasets. SoftKI approximates a kernel via softmax interpolation from a smaller number of interpolation points learned by optimizing a combination of the SoftKI marginal log-likelihood (MLL), and when needed, an approximate MLL for improved numerical stability. Consequently, it can overcome the dimensionality scaling challenges that SKI faces when interpolating from a dense and static lattice while retaining the flexibility of variational methods to adapt inducing points to the dataset. We demonstrate the effectiveness of SoftKI across various examples and show that it is competitive with other approximated GP methods when the data dimensionality is modest (around $10$).

Published in Transactions on Machine Learning Research (08/2025)
High-Dimensions Gaussian Process Regression with Soft Kernel Interpolation
Chris Camaño
camano@caltech.oo edu
Deportment of Computing and Mathematical Sciences, California Institute of Technology, Pasadena, CA, 91125, USA
Daniel Huang Base26, CA, USA
dan@basc26labs. com
Reviewed on OpenReview: https: //tp://ppppe opwrr/veen...e...e net/ form?? id- (9b2FIjjWW
Abstract
We introduce Soft Kernel Interpolation (SoftKI), a method that combines aspects of Struc- tured Kernel Interpolation (SKT) and variational inducing point methods, to achieve scal- able Gaussian Process (GP) regression
high-dinensional datasets. SoftKI approximates a kernel via softmax interpolation from a smaller number of interpolation points learned by optimizing a combination of the SoftKT marginal log-likelihood (MLL), and when needed, an approximate MLL for improved numerical stability. Consequently, can overcome the dimensionality scaling challenges that SKI faces when interpolating from a dense and static lattice while retaining the flexibility of variational methods to adapt inducing points the datasct. We demonstrate the effectiveness of SoftKT across various examples and show that it is competitive with other approximated GP methods when the data dimensionality is modest (around 10).
1 Introduction
Gaussian procossos (GPs) flexible function approximators based Bayesian inference. However, there are scaling concerns. For dataset comprising data points, constructing the x kernel (covariance) matrix required for GP inference incurs a space complexity of ( Naively, posterior inference using direct solvers costs (n²²) time, which quickly becomes infcasible for moderate n. In recent years, this shortcoming has led to the widespread adoption conjugate gradient (CG) based GPs enabling "exact" GP inference
O(n²) time (Cunningham et al., 2008 Cutajar et al., 2016). Despite these developments, exact GP regression remains challenging in practice, often requiring bespoke engineering solutions (Wang et al. 2019 et al. 2
Gardner
A popular response to this scaling issue is based on nariational inference methods Titsias et al. 2013). These approaches build a variational approximation of the posterior GP by learning the locations m < n inducing points (Quinonero- Candela Rasmussen) 2005 Snelson & Ghahramani 2005
2009
Hensman
Inducing points and their corresponding inducing variables introduce latent variables with normal priors that form a low-rank approximation of the covariance structure in the original datasct. This approach improves the time complexity of posterior inference O(um²) for Sparse Gaussian Process Regression (SGPR) Titsias 2009) and O(m²) for Stochastic Variational Gaussian Process Regression (SVGP) (Hensman et al. 2013 with the latter introducing additional variational parameters.
Another successful approach is based on Struclured Kernel Interpolation (SKT) Wilson & Nickisch 2015) and its variants such as product kernel interpolation (SKIP) Gardner et al.) 2018b). Simplex-SKI (Kapoor et al. 2021), and Sparso-Gree Grid SKI (Yadav et al. 2022). SKI-based methods achieve scalability by constructing computationally tractable approximate kernels via interpolation from pre-computed and dense rectilinear grid of interpolation points. This structure enables fast matrix-vector multiplications (MVMs), which leads to

---

Published in Transactions on Machine Learning Research (08/2025)
downstream acceleration when paired with CG-based kernel inversion strategies. Superficially, this approach effectively responds to the aforementioned scaling issues, but unfortunately has the consequence of causing the complexity of posterior inference to become explicitly dependent on the dimensionality d of the data (e.g., O(14 dm" log m) for a MVM in SKI and (nd + m)) for a MVM in Simple--SKT). Morcover, the static grid used in SKT approximations does not have the flexibility of a variational method to adapt to the dataset at hand. This limitation motivates us to scarch for accurate and scalable GP regression algorithms that can attain the best of both approaches.
In this paper, introduce soft kernel interpolation (SoftKI), which combines aspects of inducing points and SKI to enable scalable GP regression on high-dimensional datasets. Our main observation is that while SKI can support large numbers of interpolation points, it uses them in a sparse manner only a few interpolation points contribute to the interpolated value of any single data point. Consequently, we should still be able to interpolate successfully provided that we place enough of them in good locations relative to the dataset, i.e., learn their locations as in variational approach. This change removes the explicit dependence of the cost of the method on the data dimensionality d. It. also opens the door to more directly optimize with an approximate GP marginal log- ligliikeihhood (MLL) to learn the locations of interpolations points opposed to variational approximation that introduces latent variables with normal priors over the values obscrved at inducing point locations.
Our method approximates a kernel by interpolation from a softmax of m < n learned interpolation points, hence the name soft kernel interpolation (Section 3.1) The interpolation points are learned by optimizing a combination of the SoftKI MLL, and when needed, an approximate MLL for improved numerical stability in single-precision floating point arithmetic (Section 3.2 It is able to leverage GPU acceleration and stochastic optimization for scalability. Since the kernel structure is dynamic, we rely on m << interpolation points to obtain a time complexity of O(nm²) and space complexity of ((m)) for posterior inference (Section 3.3 We evaluate SoftKI on a variety of datasets from the UCI repository (Kelly et al. 2017) and demonstrate that it achieves test root mean square error (RMSE) comparable to other inducing point methods for datasets with moderate dimensionality (approximately = 10) (Section [41]. To further explore the effectiveness of SoftKI in high dimensions, we apply SoftKI to molecule datasets (d in the hundreds thousands) and show that SoftKI is scalable and competitive in these settings well (Section 4.2 Lastly, we examine the numerical stability of SoftKI (Section
2 Background and Related Work
Notation. Matrices will be denoted by upper-case boldface letters (e.g. A E R_xm). When specifying a matrix entry-wise generated by an underlying function we write A = [g(i, meaning Aij = gi. At times, given points F1,.... Fn E Rdd, we stack them into E Rnd by x" = (IT...., ,. so that entries d(i - 1) +1 through di are the components of x,. Any function f : Rd + R extends column-wise via ((x)
2.1 Gaussian Processes
Let k Rd x Rd + R be a positive semi-definite kernel function. A (centered) Gaussion process (GP) with kernel k is a distribution over functions f: Rd + R such that for any finite collection of inputs {x:)a the vector of function values fxx is jointly Gaussian fxx) ~ NO, Kxx), where the kernel matrix Kxx E Ruxn has entries [Kxlij = k( 2z)) Given two collections {rr and (a/ by the cross-covariance matrix is denoted Kxx² E Raxx' is defined by Kaeli k(r,, 2²)).
To perform GP regression on the labeled dataset D = ((1,, : Dy E IRd, yi E RR) we assume the data is generated as follows:
f(x) N(0,
(GP)
y ~ [(fxx))((((x,)))) 381)
(likelihood)
where of function sampled from a GP and each observation y/ is evaluated x; perturbed by inde- pendent and identically distributed (ii.d.) Gaussian noise with zero mean and variance B², The posterior
2

---

Published in Transactions on Machine Learning Research (08/2025)
predictive distribution has the following closed-form solution [Rasmussen & Williams, 2005
x( = (K. + - K.. K XKx(Xx + -x=)
(1)
where = 3 Using direct methods, the time complexity of inference is Onn)) which is the complexity of solving the system of linear equations in variables (Kxx A)a = y for so that the posterior mean (Equation
1 is K.x
A GP's hyperperrmeeers 0 include the noise 8 and other parameters such those involved in the definition of kernel such as its lengthscale i and output scale Thus = (3, The hyperparameters can be learned by maximizing the MLL of a GP
logply x; 0) = Ny O, Kxx + A())
where (V Σ) is notation for the probability density function (PDF) of Gaussian distribution with mean f and covariance E and we have explicitly indicated the dependence of Kxx and A on A.
2.2 Sparse Gaussian Process Regression
Many scalable GP methods address the high computational cost of GP inference by approximating the kernel using Nyström method (Williams & Seeger 2000). This approach involves selecting smaller set of inducing points, z C x, to serve as representatives for the complete dataset. The original n x kernel Kxx is then approximated as
SGPR
KxnKn Ri Kxx:
xx
leading an overall inference complexity of O(mm²).
Variational inducing point methods such as SGPR Titsias 2009) combines the Nyström GP kernel approxi- mation with variational optimization procedure so that the positions of the inducing points can be learned. The variational approximation introduces latent inducing variables f(z) to model the values observed at inducing points z with joint distribution
p(f(x), ( . Kxx Kxz
)
Kzx Kmn
The inducing points are learned by maximizing the evidence lower bound (ELBO)
ELBO(4) Z2) [lo(p(( Z)] (L(g(Z)) (Z))
(2)
where KL is the KL divergence between two probability distributions, q is computationally tractable
variational distribution, Y is observed data (Y= =y for SGPR), and Z are latent variables (Z 1z=(6 (f(xx,x,, f(z))
for SGPR). In the case of SGPR regression, the variational family f()x p(f(x) ((()))((())(((()) is chosen to approximate
so that Equation 2 simplifies to
EL) = log N (y O, K
SGPR XX
A) tr (Kxx K SGPR
Since the ELBO is lower bound on the MLL logp(y x; 0), we can maximize the ELBO via gradient-based optimization as a proxy for maximizing logp(y p(y x; () to learn the location of the inducing points z by treating it GP hyperparamoter, i.e., 0 (3,6,a,2 The time complexity of computing the ELBO Onm?)).
The posterior predictive distribution has closed form
NUU)) (.
K,z K,z(K -C-))
where C = Kyn + Kxx - The time complexity of posterior inference is Onmm² + m²) since it requires solving Ca KzxA -'y for a, which is dominated by the cost of forming C.
Additional variational methods. Building on SGPR, SVGP (Hensman al. 2013) extends the optimiza- tion process used for SGPR to stochastic variational inference. This further reduces the computational
3

---

Published in Transactions on Machine Learning Research (08/2025)
cost of GP inference to Omm²) since optimization is now over variational distribution on 2, rather than the full posterior. Importantly SVGP is highly scalable to GPUs since its optimization strategy is amenable to minibatch optimization procedures. More recent research has introduced Variational Nearest Neighbor Gaussian Processes (VNNGP) (Wu et al. 2022). which replaces the low-rank prior of SVGP with a sparse approximation of the precision matrix by retaining only correlations among each point's 's K nearest neighbors. By constructing sparse Cholesky factor with at most K + nonzeros per row, VNNGP can reduce the per-iteration cost of evaluating the SVGP objective (Equation 2.2 to ((1s + m)) when minibatching over ny data points and my inducing points.
2.3 Structured Kernel Interpolation
SKT approaches scalable GP regression by approximating a large covariance matrix Kxx by interpolating from cleverly chosen interpolation points. In particular, SKI makes the approximation
K K Wxn Kzn W R Kxx.
where Wxe is sparse matrix of interpolation weights, Wex = W xx: and z are interpolation points. The interpolation points can be interpreted as quasi-inducing points since
W
K KSKI
Kyx K (Kux Wxx) Kxz K K,x
where Rxs WxrKex ax- However, interpolation points are not associated with inducing variables since they are used for kernel interpolation rather than defining a variational approximation. The posterior predictive distribution then the standard GP posterior predictive
p([() R N/K. (K ( + +A) -,) K.x (K + -K-)
where Kxx replaced with KSKT When Equation solved via lincar CG with structured interpolation (e.g., cubic), the resulting method is called KISS-GP (Wilson & Nickisch 2015 Since the time complexity of CG depends access to fast MVMs, KTSS-GP makes two design choices that are compatible with CG. First, KISS-GP uses cubic interpolation weights (Keys) 1981] that Wxx xz has only four non-zero entries in each row, i.c., Wx is sparse. Second, KISS-GP chooses the interpolation points z to be on fixed rectilinear grid in Rdd which induces multilevel Toeplitz structure in Kzz if k( y) is a stationary kernel.
Putting the two together, the resulting kernel KSKI is a highly structured structured matrix which enables fast MVMs. If we wish have m distinct interpolation points for each dimension, the grid will contain md interpolation points leading to a total complexity of GP inference of O(n4 + dmd log(m)) since there are O(4") nonzero entries per row of W and MVM with Tocplitz structured kernel like K can be done in O(1z) log([z))) time via a Fast Fourier Transform (FFT) (Wilson) 2014). Thus, for low dimensional problems, SKI can provide significant speedups compared to vanilla GP inference and supports a larger number of intorpolation points compared to SGPR.
xx
2.4 Other Variants of Structured Kernel Interpolation.
While SKI can provide significant speedups over vanilla GP inference, its scaling in d restricts its usage to low-dimensional settings where 4d < n. To address this, Gardner et al. (2018b) introduced SKIP, which optimizes SKT by expressing a d-dimensional kernel as a product of d one-dimensional kernels. This reduces the MVM cost with Kz from (omm)) kernel instead of md. However, SKIP limited to dimensions roughly in the range of 10 30 for large datasets, it requires substantial memory and may suffer from low-rank approximation errors.
Odm logm))) using only m grid points per component
Recent work on improving the scaling propertics of SKI has been focused on efficient incorporation of sim- plicial interpolation, which when implemented carefully in the case of Simplex-GP (Kapoor et al. 2021) and Sparse- Grid GP (Yadav et al. 2022) can reduce the cost of a single MVM with Wxz down to O (nd²). These methods differ in the grid construction used for interpolation with Simplex-She opting for permutohedral lattice (Adams et al. 2010), while Sparse- Grid GP uses sparse grids (Bungartz & Griebel 2004) and cus- tom recursive MVM algorithm. These modifications vastly improve KISS-GP leading to Kya MVM costs of
4

---

Published in Transactions on Machine Learning Research (08/2025)
Table 1: Cost of a single MVM with the approximate covariance for various GP approximations. For listed methods based SKI take points in each dimension. The variable represents the rank used in the SKIP approximation (r s 100), and the variable l is the grid resolution taken in Sparse-Grid Grid GP (typically x 11-5) 5)
Method
Time Complexity of one MVM
KISS-GP (cubic)
(n44 + dmm log(m))
SKIP
(d(rn+ m mlog m)) (d m))
Simplex-Gr
Sparso-Grid GP (simplicial)
0 (nd2+e 2222)
SGPR
O(nm² + m
SoftKI
0 Onnm² + m²)
(n m)) for Simplex-GP and O (nd=+++ + (d2t) for Sparse-Grid GP where is a small constant. While both of these methods greatly improve the scaling of SKI, they still feature a problematic dependency on d limiting their extension to arbitrarily high dimensional datasets. The complexity of SKI and related variants are tabulated in Table 1
3 Soft Kernel Interpolation
In this section, we introduce SoftKI (Algorithm 1 SoftKI takes the same starting point as SKI, namely that KSKI W xx Kz Was ax can be used as an approximation of Kxx- However, we will deviate that we will abandon the structure given by lattice and opt learn the locations of the interpolation points instead. This raiscs several issues. First, revisit the choice of interpolation scheme since static structure (Section 3.1. Second, we require efficient procedure for learning the interpolation points (Section 3.2 Third, we require an alternative route to recover a scalable GP since our approach removes structure in the kernel that was used in SKI for efficient inference (Section 3.3
longer assume
3.1 Soft Interpolation Strategies
Whereas cubic (or higher-order) interpolation scheme is natural in SKI when using static lattice structure, our situation is different since the locations of the interpolation points are now dynamic. In particular, we would ideally want our interpolation scheme to have the flexibility to adjust the contributions of interpolation points used to interpolate a kernel value in a data-dependent manner. To accomplish this, we propose softmar interpolation weights.
Definition 3.1 (Softmax interpolation weights). For a finite collection of inputs {a:) H; E Rd, and a set of interpolation points { 2z E IRd, define the softmar interpolation weight matrix as
exp (- || 13
Ex2
(4)
exp (- | 24
With this definition in hand we can now define the SoftKI kernel approximation
SoftKI xx
Exg Epx R Kxx-
Since a softmax can be interpreted as probability distribution, we can view SoftKI as interpolating from a probabilistic mixture of interpolation points. Each input point x; is interpolated from all points, with an exponential weighting favoring closer interpolation points (See Figure [ Each row Σr,z represents the softmax interpolation weights between a single point J; and all m interpolation points. Note that Exz is not strictly sparse since softmax interpolation continuousl assigns weights over each interpolation point for each Zt- However, d increases, the weights for distant points become negligible. This leads to
5

---

Published in Transactions on Machine Learning Research (08/2025)
(44 Structured interpolation Points
Cubic Interpolation Weights
Leumed interpolation Points
Softmax Interpolation Weights
00
Dats Points
Date Points
Interpolstion
Interpolation
Figure Structured & Soft Kernel Interpolation: Comparison of interpolation procedure for local cubic interpolation during KISS-GP (left) and softmax interpolation (right) during SoftKI. Here white points indicate interpolation points 2, and black points are data points 2; E In these diagrams magenta. line segments indicate which interpolation points are being used for a given data point during the interpolation procodure of cach method. KISS-GP, a small local subset of grid points is used, whereas in KI all interpolation points are included, weighted by a softmax distribution centered at given data point.
effective sparsity, with each 2, predominantly influenced only by its nearest interpolation points during the roconstruction of Kxx-
To add more flexibility to the model, the softmax interpolation scheme can be exlended to contain learnable temperature parameter I' as below
|. / 25|1) k= exp (- 11: 2 II)
exp
Exs
This additional hyperparameter is needed because unlike the SGPR kernel K SGPR KxeKKiKK K2x where the lengthscale influences each term, the lengthscale only affects Kz in the SoftKI kernel. The temperature acts as lengthscale-Hi hyperparameter on the interpolation scheme that controls the distance between the datapoints x and the interpolation points z. When = 1, we obtain the original scheme give in Equation Analogous to how automatic relevance detection (ARD) (MacKay et al., 1994) can be used set lengthscales for different dimensions, we can also use different temperature per dimension. When learning both temperatures and lengthscales, cap the lengthscale range for more stable hyperparameter optimization since they have push-and-pull effect.
3.2 Learning Interpolation Points
Similar SGPR, we can treat the interpolation points z hyperparameters of GP and learn them by optimizing appropriate objective. We propose optimizing the MLL log p(y () with method based on gradient descent for SoftKT which for D, (SfttKI(0) + A()) has closed the form solution
1 Do y log llgg(22)
logply x;
(5)
and derivative
i Olog x)
0D, D + tr D 2D, 00 )
D
(6)
00
a0
Each evaluation of logp(y x; 0) has time complexity O(nm²) and space complexity O(mm) where m is the number of interpolation points.
Hutchinson's pseudoloss. Occasionally, we observe that the MLL is numerically unstable in single- precision floating point arithmetic (see Section 4.3 for further details). To help stabilize the MLL in these

---

Published in Transactions on Machine Learning Research (08/202))
Leamned Induxing Fritt
Rezd inteppiation Roims
Leamnd Inte poracioion Peints
SoftKii
0.5
0.80
0.15
3.5
0.
-0.30
64
-0.5
22
SGPR Absolute trror
SKI Absolute Error
SoftKI Absolute tmrr
10
20
0.0
0.5
Abroc dee
1.0
0.5
10
10
20-20
2
Figure 2: Inducing vs. Interpolation Points: Comparison of inducing points learned by SGPR, static rectilinear lattice points in SKI, and adaptive interpolation points learned by SoftKl on the Ricker wavelet [Ricker] 1953) after 100 epochs of hyperparameter optimization. Contour plots of the posterior mean are paired with each method's absolute-error. Each method achieves comparable accuracy: SGPR's inducing points follow the sample distribution, and SoftKT's interpolation points adapt to the true function's local geometry. Additional experiment details can be found in Appendix C.1
situations, use approximate MLL based work in approximate GP theory (Gardner et al. 2018b; Maddox et al. 2021 Wenger et al. 2023) that identify efficiently computable estimates of the log determinant (Equation and trace term (Equation 6) in the GP MLL. In more detail, these methods combine stochastic trace estimation via the Hutchinson's trace estimator GGrard 1989 Hutchinson 1989) with blocked CG so that the overall cost remains quadratic in the size of the kernel. Maddox et al. (2022) term this approximate MLL the Hutchinson's pseudoloss and show that remains stable under low-precision conditions.
Definition 3.2 (Hutchinson Pseudoloss). Let un, u1,.... ur be the solutions obtained by using block conju- gate gradients for the system Do Un H₁ ..... u) (y w.... ww...wi), where each w, E R" is Gaussian random vector normalized to unit length. The Hutchinson pseudoloss approximation of the GP marginal log likeli- hood in Equation is given as
log ply xx
un D Ug +
Σ u (Dow,)
(7)
with derivative
i log ply xx;
1
0D, ug
1 Σ 1
DD, 0 w
110
(8)
00
2
00
00
By computing Equation 7 as described above as opposed to explicitly computing a stochastic Lanezos quadrature approximation (Ubaru et. al.] 2017), the gradient can also be efficiently computed using back propagation. Note that the gradient of Hutchinson's pseudoloss approximates the gradient the GP MLL, where the trace term is estimated using stochastic trace estimator. The accuracy of this approximation depends on the number of probe vectors t. Consoquently, gradient- optimization of the GP MIL can be approximated by optimizing Hutchinson's pseudoloss instead.

---

Published in Transactions on Machine Learning Research (08/2025)
Algorithm 1 SoftKI Regression. The procedure kmeans performs k-means clustering, batch splits the dataset into batches, and sof toftmax____ swwftaa_tttttteeeaeeeeaaeeeeeeee iataaaatttteoootiiiiooo produces a softmax interpolation matrix (see Section 3.1
Require: SoftKT GP hyperparameterss 0 = (33, f, E (R(xdd) T), kernel function k(x,
Require: Dataset D = {x (x(xx {.
Require: Optimization hyperparameters: batch size and learning rate
Ensure: Learned KT coefficients a.
Model Training
z kmeans(x,, (x,m) m)
for i= to epochs do
for (xb, Vb) in batch(D,b) b) do
4:
Σx,x f sof tmax interpolationn(x..22
5: K f- [k(2,; 2z) lis
SoftKI t- Exz ΣT
6:
Xn
xaz
7: + +1 11 Vo log 70 yb
X6:
SoftKK
XtXb
+ A)
b. Stabilized MLL (Section 3.2
end for
end for
Stabilized Inference (Section
10: Exz softmax interpollation(x z))
11: U U U22 t cholesky (Kzz)
A -1/2 Exa Knz
12: Q R. + QR
U
A A-1
13: a a R R-
14: return
Stabilized MLL. We combine SoftKI's MLL with Hutchinson's pseudoloss to arrive at the objective
log 0) = llog(y x; 0) when stable
log p(y x: otherwise
for optimizing SoftKl's hyperparameters. In practice, this means that we default to using SoftKI's MLL and fallback to Hutchinson's ' pseudoloss when numeric instability is encountered. This enables SoftKI to accurately recover approximations the MLL's gradient during gradient- based optimization, even when the current positioning of interpolation points would otherwise make the direct computation of the MLL's gradient unstable, while relying on the exact MLL as much as possible.
compute V Vlog ply Xx; 0) on minibatch of data of size b to perform stochastic optimization (sce line 7 of Algo- rithmn The stochastic estimate provides an unbiased estimator of the gradient. In this way, SoftKI can leverage powerful stochastic optimization techniques used to train neural networks such as Adam (Kingma & Ba 2014), and hardware acceleration such as graphics processing units (GPUs). The time complexity of evaluating V log p(y x):0) is cheap, costing mb) operations with space complexity (bm).. The lower GPUs since they have memory constraints compared to CPUs. performing stochastic gradient descent and not stochastic variational inference in SVGP (Hensman et al. 2013). In particular, we do not define distribution the interpolation points nor make variational approximation.
Stochastic optimization. Instead of computing V log p(y xx;0) the entire dataset
space complexity makes it easier
We emphasize that
3.3 Posterior Inference
Because K SoitKI is low-rank, can use the matrix inversion lemma to rewrite the posterior predictive of SoftKl as
- 'R -1,, SottKI K SoftKl (A -1 A A- -'Rey - - R -1 -1) Softkl)
8

---

Published in Transactions on Machine Learning Research (08/2025)
where Č = Kez KxxA AxxA-KKa "Kxx (see Appendix Al. To perform posterior mean inference, solve
Ca A
(9)
for weights Note that this is the posterior mean of SGPR with C replaced with Č and Kx replaced with Kxz- Moreover, note that C is simply C with Kxa replaced with Kx Since C is m m matrix, the solution of the system of linear equations has time complexity O(m²). The formation of C requires the multiplication of a m X n matrix with a n x m matrix which has time complexity O(nm²). Thus the complexity of SoftKI posterior mean inference is (nm²² since it is dominated by the formation of C. We refer the reader to the supplementary material for discussion of the posterior covariance (Appendix A
Solving with QR. Unfortumately, solving Equation for a can be numerically unstable. Foster et al. 20099
introduce a stable QR solver approach for a Subset of Regressors (SoR) GP (Smola & Bartlett, 2000) which we adapt to a SoftKI. Since the form of the SoftKI and SGPR posterior are similar, this QR solver approach also be adapted to improve the performance of SGPR (see Appendix B
Define the block matrix
A A Uzz
A
where U U.U Un Kun the upper triangular Cholesky decomposition of Kez so that
A A A RAx - + K =.
Let QR A be the QR decomposition of so that Q is a (n+m) x m orthonormal matrix and is m x m right triangular matrix. Then
Ca A
A (A-1/3
(QR)T(QR)a= = (QR)
(+)
A
11113)
Ra = Q
(10)
0
We thus solve Equation 10 for a via a triangular solve in O(m²) time. For additional experiments demon- strating the inability of other linear solvers to ofler comparable accuracy to the QR-stabilized solve detailed in this section, see Appendix
4 Experiments
We compare the performance of SIft against popular scalable GPs on selected data sets from the UCI data set repository (Kelly et al.] [0017], common GP benchmark (Section 4.1). Next, we test Soft on high-dimensional molecule data sets from the domain of computational chemistry (Section 4.2) Finally, we
explore the numerical stability of SoftKI (Section 4.3
4.1 Benchmark on UCI Regression
We evaluate the efficacy of SoftKI against other scalable GP methods on data sets of varying size n and data dimensionality d from the UCI repository (Kelly et al.] 2017). We choose SGPR and SVGP as two scalable GP methods since these methods applied to many data sets. It is not possible to apply SKI to most datasets that we test on due to scalability issues with data dimensionality, and so we omit it. For additional comparisons to other alternative SKI architectures, see Appendix
be applied in a relatively blackbox fashion, and thus, be
Experiment details. For this experiment, we split the data set into 0.9 for training and 0.1 for testing. We standardize the data have mean and standard deviation using the training data set. We use a Matérn 3/2 kernel and learnable output scale for each dimension (ARD). We choose 512 inducing

---

Published in Transactions on Machine Learning Research (08/2025)
Dataset
T
Exact GP
Soft m = 512
SGPR m = 512
SVGP m = 1024
3droad
391386
0.583 + 0.010
0.389 + 0.001
Kin40k
36000
8
0.039 + 0.001
0.169 + 0.008
0.177 0.004
0.165 + 0.005
Protein
41157
0.044 + 0.000
0.596 + 0.016
0.602 0.015
0.607 + 0.012
Houseelectric
1844352
11
0.047 0.001
0.047 = 0.000
Bike
15641
17
0.040
0.003
0.062 = 0.001
0.108 0.004
0.084 0.006
Elevators
14939
18
0.108 0.014
0.360 + 0.006
0.395 + 0.005
0.384 = 0.008
Keggdirected
43944
20
0.070 0.003
0.080 0.005
1.0000
0.082 = 0.004
Pol
13500
26
0.035 + 0.001
0.075 + 0.002
0.127 1 0.002
0.122 0.002
Keggundirected
57247
27
0.110 0.004
0.115 = 0.004
0.128 = 0.009
0.121 = 0.007
Buzz
524925
77
0.240 + 0.001
0.250 = 0.002
Song
270000
90
0.777 + 0.004
0.794 = 0.006
Slice
48150
385
0.022 = 0.006
0.520 = 0.001
0.082 = 0.001
Table 2: Test RMSE on UCI datasets. Best results obtained by approximate GP are bolded. Exact GP included for reference. Entries marked indicate datasets that triggered out of memory errors.
Datasct
d Exact GP
SoftKI m 512
SGPR = 512
SVGP m = 1024
3droad 391386
0.953 0.041
0.597 0.008
Kin40k
36000
1.636 + 0.006
0.055 0.043
-0.104 = 0.007
-0.082 + 0.007
Protein
41157
0.904 + 0.008
0.905 + 0.026
1.031 = 0.009
1.047 + 0.009
Houseelectric
1844352
11
1.274
-1.492 + 0.007
Bike
15641
17
1.662 + 0.065
- -0.379 0.016
-02996 = 0.011
-0.680 + 0.017
Elevators
14939
18
0.590
+ 0.075
0.406 0.021
0.587 = 0.006
0.586 0.009
Keggdirected
43944
20
nan
0.421 0.063
-08996 = 0.041
0..934 0.017
Pol
13500
26
1.750 + 0.020
-0.710 0.028
-0.428 + 0.002
-0.94 + 0.012
Keggundirected
57247
27
nan
0.302 + 0.136
-0.99 = 0.039
-0.590 + 0.021
Buzz
524925
77
0.276 + 0.348
0.130 + 0.008
Song
270000
90
1.179 + 0.008
1.288 = 0.002
Slice
48150
385
1.258 + 0.149
1.330 = 0.002
00.662 + 0.002
Table 3: Test NLL on UCI datasets. Best approximate results obtained by an approximate GP bolded.
Exact GP included for reference. Entries marked '_- indicate datasets that triggered out of memory errors. Entries marked nan indicate that the computations were unstable.
points for SoftKI. Following Wang al. (2019). we use m 512 for SGPR and m 1024 for SVGP. We learn model hyperparameters for SoftKI by maximizing Hutchinson's pseudoloss, and for SGPR and SVGP by maximizing the ELBO. For reference, we also include results for GP fit. with block CG, termed an Fxact GP. For Exact GP we use the KeOps (Charlier et al. 2021] library integration with GPyTorch to handle the additional scaling challenges of working with the full n X n kernel.
We perform 50 epochs of training using the Adam optimizer (Kingma & Ba 2014) for all methods with a learning rate of 7 0.01. The lcarning rate for SGPR is n 0.1 since we are not performing batching. We use default implementation of SGPR and SVGP from GPyTorch. For SofLKI and SVGP, use minibatch size of 1024. We report the average test RMSE (Table 2) and NLL (Table 3) for each method across three seeds. Additional timing information is given in Appendix
Results. We observe that SoftKT has competilive Lest RMSE performance compared to SGPR and SVGP, exceeding them when the dimension is modest (d æ 10). For the test NLL, we observe cases where SoftKI's test NLL smaller or larger relative to the test RMSE. As reminder, the NLL for a GP method consists of two components: the RMSE and the complexity of the model. Consequently, a small test NLL relative to the test RMSE indicates relatively small amount of noise in the data set. Conversely, large test NI.I.
10

---

Published in Transactions on Machine Learning Research (08/2025)
Dataset
71
d
Exact GP
SoftKT m=512
SGPR m=512
SVGP m=1024
Ac-ala3-nhme
76598
126
0.017 + 0.000
0.790 + 0.010
0.852 0.005
0.836 + 0.005
Dha
62777
168
0.016 0.000
0.886 0.012
0.897 0.012
0.882 + 0.010
Stachyose
24544
261
+
0.363 0.012
0.643 + 0.002
0.555 + 0.001
At-at
18000
354
0.020 1020
0.528 0.011
0.591 =
0.558 + 0.009
ag 9137
354
0.021
0.502
0.565 0.029
0.461 0.026
Buckyball-catcher
5491
444
0.018 0.001
0.153 0.006
0.393 + 0.008
0.307 0.015
na-otuue
4528
1110
0.031 + 0.001
1.001 = 0.048
0.045 + 0.000
Table 4: Test RMSE on MD22 datasets. Best results obtained by an approximate GP are bolded. Exact GP
results are included for reference. Exact GPs are roughly 20 875x slower than SoftKI (see Table 10 and thus, only run for 50 epochs.
Dataset
d
Exact GP
SoftKT m=512
SGPR m=512
SVGP m=1024
Ac-ala3-nhme
76598
126
0.000 + 0.000
1.188 = 0.014
1.367 + 0.003
1.356 = 0.002
Dha
62777
168
-0.975 = 0.844
1.299 = 0.014
1.418 = 0.006
1.409 = 0.007
Stachyose
24544
261
0.003 -1.330
0.522 = 0.059
1.123 + 0.002
1.012 + 0.003
At-at
18000
354
0.795 = 0.007
1.041 + 0.006
1.003 = 0.008
At-at-cg-og 9137
354
-1.426
0.000
0.742 = 0.010
0.990 = 0.024
0.800 = 0.016
Buckyball-catcher
5491
444
-1.567
0.003
-0.278 = 0.061
0.689 = 0.005
0.449 0.015
DW-nanotube
4528
1110
-1.992 + 0.000
-1.551 + 0.049
1.520 + 0.021
-0.954 + 0.009
Table 5: Test NLL on MD22 datasets. Best results obtained by an approximate GP are bolded. Exact GP
results are included for reference. Exact GPs are roughly 20 875x slower than SoftKI (sce Table 10 and thus, only run for 50 epochs.
relative the test RMSE indicates a large amount of noise in the datasct. We can see instances of this the bike dataset where SoftKI achieves comparable lest RMSE to SGPR and SVGP, but results in larger NLL. This indicates that SoftKI's kernel is more complex compared to the variational GP kernels.
4.2 High-Dimensional Molecule Dataset
Since SoftKT's are performant on high dimensional data sets from the UCI repository, we also test the performance on molecular potential energy surface data on the MD22 (Chmiela et al., (023) dataset. These are high-dimensional datasets that give the geometric coordinates of atomic nuclei in biomolecules and their respective energies.
Experiment details. For consistency, we keep the experimental setup the same for UCI regression but the number of epochs of training to 200 due to slower convergence on these datasets. We still use the Matérn kernel and not a molecule-specific kernel that incorporates additional information such as the atomic number of each atom (e.g., a Hydrogen atom) or invariances (e.g., rotational invariance). We standardize the data set to have mcan and standard deviation 1. We note that in an actual application of GP regression to this setting, we may only opt to center the largets to be mean 0. This is because energy is relative number that can be arbitrarily shifted, whereas scaling the distances between atoms will affect the physics. As before, we include results for Exact GPs as reference. Since they are roughly 20 to 875 times slower than SoftKI, so we only run them for 50 epochs (see Table 10
Results. Table 4 and Table 5 compares the test RMSE and test NLL of various GP models trained on MD22. We see that SoftKI competitive approximate method on these datasets, especially on the test NLL. GPs that fit forces have been successfully applied fit such data chemical accuracy. In our case, we do not fit. derivative information.
11

---

Published in Transactions on Machine Learning Research (08/2025)
Test R.MSEE
Test NLL
Dataset
MLL
Hutchinson Pseduoloss
MLL
Hutchinson Pseduoloss
3droad
nan
0.541 + 0.005
nan 0.055 + 0.043
1.004 + 0.175
Kin40k
0.169 = 0.008
0.183 + 0.006
0.07 = 0.029
Protein
0.596 = 0.016
0.602 + 0.016
0.904 + 0.03
0.914 + 0.028
Houseelectric
nan nan
0.049 + 0.002
nan nan
1.275 + 1.103
Bike
0.073 + 0.005
-0.101 + 0.179
Elevators
0.36 + 0.007
0.36 = 0.007
0.404 + 0.024
0.404 + 0.025
Keggdirected
0.08 + 0.005
0.082 0.004
0.293 0.169
0.628 = 0.289
Pol
0.075 + 0.002
0.084 + 0.001
-0.694 + 0.06
-0.783 + 0.033
Keggundirected
0.114 = 0.005
0.118 + 0.004
0.308 + 0.025
0.291 = 0.089
Song nan
Buzz
0.24 = 0.0
nan 1.328 0.045
0.599 + 0.555
0.777 = 0.004
0.777 + 0.004
1.179 + 0.008
1.178 + 0.007
Slice
0.022 + 0.006
0.044 0.005
0.607 + 0.058
Table 6: Comparison of SoftKI test RMSE and NLL when training with the SoftKI MLL logp(y x; O) versus the Hutchinson pseudoloss logí(y x; 0) Best results are bolded. Entries labeled nan, represent instances where using the exact MLL resulted in numerical instability in SoftKI (see Section 4.3)
4.3 Numerical Stability
In Section 3.2 we advocated for the use of Hutchinson's pscudoloss to overcome numerical stability issues that arise when calculating the SoftKI MLL. This adjustment is specifically address situations where the matrix S + A is not positive semi-definite in float32 precision. In other approximate GP models implemented in GPyTorch, this challenge is typically addressed by performing Cholesky decomposition followed by efficient low-rank computation of the log determinant through the matrix determinant lemma.
In our experience, even when using versions of the Cholesky decomposition that add additional jitter along the diagonal there are still situations where Kzz can be poorly conditioned, particularly when n large. We conjecture that this behavior originates from situations where the learned interpolation points of a SoftKI coincide at similar positions, driving the effective rank of Kzz down. In these situations, Hutchinson's 's pseudoloss offers a more stable alternative because it does not directly rely on the matrix being invertible.
In Table 6 we replicate the UCI experiment of Section 4.1 using the SoftKI MLL log x; and Hutchin- son's pseudoloss p(y 0) For most datasets, optimizing with the SoftKI MLL produces better test RMSE. However, more numericaly instability is also encountered, with SoftKI failing on four datascts that we lested. In these cases, we find that Cholesky decomposition fails.
5 Conclusion
In this paper, we introduce SoftKT, an approximate GP designed for regression on large and high-dimensional datasots. SoftKI combines aspects of SKI and inducing points methods to retain the benefits of kernel interpolation while also scaling to higher dimensional datasets. We have lested SoftKI on variety of datasets and shown that it is possible to perform kernel interpolation in high dimensional spaces in way that is competitive with other approximate GP abstractions that leverage inducing points. Exploring methods that enforce stricter sparsity in the interpolation matrices, such as through thresholding, could enable the use of sparse matrix algebra further improving the cost. of inference.
Acknowledgements Chris Camaño acknowledges support by the National Science Foundation Graduate Research Fellowship under Grant No. 2139433, and the Kortschak scholars program. Additionally, we would like to extend appreciation to Ethan N. Epperly for his thoughtful comments and advice on this work.
12

---

M.F. Hutchinson. A stochastic estimator of the trace of the inffuence matrix for laplacian smooth- ing splines. Communiccttin in Slatistics- Simulation and Compulation, 18: 18:1059 1076, 01 1989. doi: 11108000036109910088186666 Published in Transactions on Machine Learning Research (08/2025)
References
Andrew Adams, Jongmin Back, and Myers Davis. Fast high-dimensional filtering using the permutohedral lattice. Compul. Gruph. Forum, 29:753 762, 05 2010. doi: 10.111/1:1147--8552220000 200...00.000..1.....
Hans-Joachim Bungartz and Michael Griebel. Sparse grids. Acta numerica, 13:147-269, 2004. doi: 10.1017/ 1569995.000999
Benjamin Charlier, Jean Feydy, Joan Alexis Glaunès, François- David Collin, and Ghislain Durif. Kernel operations on the gpu, with autodiff, without memory overflows. Journal of Machine Learning Research, 22(74), 2021. URL http::/jmmt http://ir...rrrppppe//////////20
html
Stefan Chmiela, Valentin Vassilev-Galindo Oliver T Unke, Adil Kabylda, Huziel E Sauceda, Alexandre Tkatchenko, and Klaus-Robert Müller. Accurate global machine learning force fields for molecules with hundreds of atoms. Science Advances, 9(2), 2023. doi: 10. 10.1126/sciadv. ad00873.
John P Cunningham, Krishna V Shenoy, and Mancesh Sahani. Fast gaussian process mothods for point process intensity estimation. In Proceedings of the 25th international conference on Machine learning, pp. 192-199, 2008. doi: 10.1 10.1145/1390155 133000.88
Kurt Cutajar, Michael Osborne, John Cunningham, and Maurizio Filippone. Preconditioning kernel ma- trices. In International conference on machine learning. pp. //proceedings.co mmr.iingg.m.......m wtt.lr.rree////////ue
PMLR, 2016. URL https:
html
Leslie Foster, Alex Waagen, Nabeela Aijaz, Michael Hurley, Apolonio Luis, Joel Rinsky, Chandrika Satyavolu, Michael J. Way, Paul Gazis, and Ashok Srivastava. Stable and efficient gaussian process calculations. Journal of Machine Learning Research, 10(31):857 882, 2009. URL httt:://ialr.org/ppperrs//// org/papers/v10/
FosterO9a. html
Jacob Gardner, Geoff Pleiss, Kilian Weinberger, David Bindel, and Andrew G Wilson. Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration. Advances in neural informa-
tion processing systems, 31, 2018a. URL https://papers www..iiiss.....p.ee wttpcccccppcacaaeeeee fil.es//ppperr/21//aaall 27631779440783306447777722222222 -
Jacob Gardner, Geoff Pleiss, Ruihan Wu, Kilian Weinberger, and Andrew Wilson. Product kernel interpola- tion for scalable gaussian processes. In International Conference on Artificial Intelligence and Statistics, pp. 1407-1416. PMLR, 2018b. URL https::/pproceeding. http:/aww........eeeeeeeeeemeeee pttps///////araeeeerrr.eeee html
A. Girard. A fast 'monte-carlo cross-validation procedure for large least squares problems with noisy data. Numer. Math., 56(1):1-23 January 1989. ISSN 0029-599X. doi:
James Hensman, Nicolo Fusi, and Neil D. Lawrence. Gaussian processes for big data, 2013. URL https: //arxiv. www.:org/aabo http//9bb/330......
Sanyam Kapoor, Marc Finzi, Ke Alexander Wang, and Andrew Gordon Gordon Wilson. Skiing on sim- plices: Kernel interpolation on the permutohedral lattice for scalable gaussian processes. In International Conference Machine Learning, pp.
PMLR, 2021. URL httpp:://proceeedngggn mlr press/
v139/kapoor21a. html
Markelle Kelly, Rachel Longjohn, and Kolby Nottingham. The uci machine learning repository. https:
//rrchive. ics. wwc.suci.eee. edu. 2017.
R. Keys. Cubic convolution interpolation for digital image processing. IEEE Transactions on Acoustics, Speech, and Signal Processing, 1153)111, 1981. doi: 10.1109/TASSP.1985 1163711...11111111
13

---

Published in Transactions on Machine Learning Research (08/2025)
Diederik P. Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arXiv:1412. 6980, 2014. httts:://aaxiiiii http:/oorr/bb//11222/// 111.1/112......
ar preprint
David JC MacKay et al. Bayesian nonlinear modeling for the prediction competition. ASHRAE transactions. 100(2):1053 1062, 1994. doi:
18.
Wesley J. Maddox, Sanyam Kapoor, and Andrew Gordon Wilson. When are iterative gaussian processes reliably accurate?, 2021. URL.
wwwgg/as//22.... 15246
Wesley J Maddox, Andres Potapcynski, and Andrew Gordon Wilson. Low-precision arithmetic for fast
gaussian processes. In Uncertainty in Artificial Intelligence, pp. 1306 1316. PMLR, 2022. URL https:
//proceedings.nog mlr wtw../...rr.//////////022222 html
Joaquin Quinonero- Caunellro-Candela and Carl Edward Rasmussen, A unifying view of sparse approximate gaussian process regression. The Journal of Machine Learning Research, 6:1939 1959, 2005. URL https://jmlr:
wat.aionnn.ooocaaaaaaaa..... prggppaperss////uineoee qwwiiuuroonooaaaeeeee candela05a. html
Carl Edward Rasmussen and Christopher K. I. Williams. Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning). The MIT Press, 2005. ISBN 026218253X. doi: 10.7551 /mmtpress/ 3206.001.0001.
Norman Ricker. The form and laws of propagation of scismic wavelets. Geophysics, 18(1):10-40, 1953. doi: 10.1199/1.1437843.
Alex Smola and Peter Bartlett. Sparse greedy gaussian process regression. In T. Leen, T. Dietterich, and V. Tresp (eds.), Advances in Newral Information Processing Systems, vol- ume 13. MIT Press, 2000. URL httpp::/ www.eddiinnnnnuuiss.... www..ps.cc.pppprrrr.......... httppppesllleeppppr////////////e B31466000D6637777002255553333223
pdf
Edward Snelson and Zoubin Ghahramani. Sparse gaussian processes using pseudo-inputs. Advances in neural information processing systems, 18, 2005. URL httpss//papers: https//ccppppeerre/eppppppeee
html
Michalis Titsias. Variational learning of inducing variables in sparse gaussian processes. In Artificial intelli-
gence and statistics, pp. 567-574. PMLR, 2009. UtL/ppp::/prooeedng. mlr, press//55/titsias000 http:://oooooooooooo..ey/ss////s
html
Shashanka Ubaru, Jie Chen, and Yousef Saad. Fast Estimation of (r([(a)) via Stochastic Lanezos Quadrature. SIAM Journal on Matriz Analysis and Applications, 38(4):1075-100 1099, January 2017. ISSN 0895-4798, 1095-7162. doi: 10.1137/16M11044
Ke Wang, Geoff Pleiss, Jacob Gardner, Stephen Tyrce, Kilian Q Weinberger, and Andrew Gor- don Wilson. Exact gaussian processes on a million data points. Advonces in neural infor- mation processing systems, 32, 2019. URL http:://prroeedng....rr.s....... http://.......cccooppp/////////a html
Jonathan Wenger, Geoff Pleiss, Marvin Pförtner, Philipp Hennig, and John P. Cunningham. Posterior and computational uncertainty in gaussian processes, 2023. URL https://arxiv.c.
15449
Christopher K. I. Williams and Matthias Seeger. Using the nyström method to speed up kernel machines.
In Proceedings of the 13th International Conference on Newral Information Processing Systems, NIPS'00, pp. 661-667, Cambridge, MA, USA, 2000. MIT Press. URL https://papee: http:////naa//tttaaaaaaaaeeeeee.
catsscccpppppe files/
Andrew Wilson and Hannes Nickisch. Kernel interpolation for scalable structured gaussian processos (kiss- gp). In International conference on machine learning, pp. 1775-1784. PMLR, 2015. URL https:// httpa//iw..mm...eeee//////////// preas//37/silsonttt html
proceedings mlr.
14

---

Published in Transactions on Machine Learning Research (08/2025)
Andrew Gordon Wilson. Covariance kernels for fast automatic pattern discovery and ectrapolation with Goussion processes. PhD thesis, University of Cambridge Cambridge, UK, 2014. URL https://www.c.
cmu. cat.echhh/aanndeerr////dtreeeess
Luhuan Wu, Geoff Pleiss, and John P Cunningham. Variational nearest neighbor gaussian process. In Inter- national Conference on Machine Learning, pp. 24114-24130. PMLR, 2022. URL httpss://pocceedings.
mlr. wttppp/eevv/vv2222. http:/m222hhh.mm1
Mohit Yadav, Daniel Sheldon, and Cameron Musco. Kernel interpolation with sparse grids. Advonces in
Newral Information Processing Systems, 35:22883 22894, 2022. URL http://////iiaaappppppaaaaaaaeea
html
A
SoftKI Posterior Predictive Derivation
As reminder, the SoftKI postcrior predictive is
p([()) = N(K.2 Č- A
- -
SoftK SoftKT
- - S
(11)
where C Kez + K2xA Kxz. To derive this, we first recall how to apply the matrix inversion lemma to the posterior mcan of SGPR, which as reminder uses the covariance approximation Kxx = KSGPR Ky Kzx
Lemma A.1 (Matrix inversion with GPs).
Kuz(Ka K2xA -Kxx) "KxA KSSPP KK
SGPR
A)-1
(12)
Proof.
SCPR Xx
K SGPR
(identity)
XX
+ I = (K SGPR + A) - K SGPR
+ (K SGPR + A)-1A
(mult by (KSCPR A) A)-1)
xx
I SKSG A)
K SGPR
SGPR + A)
(rearrange)
SCPR K
SGPR
SGPR A)
SGPR XX
SGPR
SGPR
AA)-AA
x
x
xx
xx
(mult both sides by SGPR on left) xx
- K.rk-KK. - K.rK-JKa (aaxxKKae KYr + A) -K 22
K
SGPR (K SGPR
+ +A)
xx
(defn xx
SSPRR)
K. (K - K-IKex K
K- Kzx
A)
x K-1 K = .x
(KCPRIKSGPP SGPR AA-
(factor)
xx
K.z (kzz + Kax Kxz)
Knx =k
SSPR xx
+ +A) -A
(matrix inversion lemma)
Kuz + K (xxx)
KExA-1
SCPR
(K
SGPR
+ -A)
(mult A on right)
xx
Now, we extend the matrix inversion lemma to SIftKI.
Lemma A.2 (Matrix inversion with interpolation).
KKe + A (Kree) Kxx--1
Softkl SoftKI XX
A)-1
(13)
Proof. Recall SoftKI = ExeKre -xecKe Ka K Exx
Krek Kex
The result follows by applying
Lemma A.1 with Kxx replaced with K
15

---

Published in Transactions on Machine Learning Research (08/2025)
RMSE per Epoch
RMSE per Epoch
0.14
CG Solver 1e-4
CG Solver 1e-4
0.80
CG Solver 1e-3
0.13
CG Solver 1e-3
Solver le-2
Solver 1e-2
CG Solver le-1
CG Solver 1e-1
0.75
Cholesky Solver
0.12
Cholesky Solver
Direct Solver
Direct Solver
RMSE
QR Solve
RMSE
0.11
QR Solve
0.70
0.10
0.65
0.09
0.60
0.08
10
20
30
40
50
0
10
20
30
40
50
Epoch
Epoch
(a) Protein.
(b) Keggdirected.
Figure 3: Test RMSE of SoftKI models trained using different linear solvers described in Section on the keggundirected dataset (left) and the protein dataset (right).
The posterior covariance is computed by a simple application of the matrix inversion lemma. To compute it, observe that the same procodure for solving the posterior mean can also used for solving the posterior covariance by replacing y with KSftKI More concretely, x+
K SoftKI
KiSftKT A
K Ò-
R -1 K
SoftKT x
(14)
SoftKI
KSfftKl, A
K
SoftKI
SoftKI A Kxz - A -K
(15)
SofKKI Sofftkl
SoftKl
Softkl xx
Kxza,
(16)
x
where
Ca. - SoftKI
(17)
be solved for a, in the same way solved for the posterior mean (with y instead of Kotkk]) Since this depends the inference point cannot precompute the result ahead of time as could with the posterior mean. Nevertheless, the intermediate results of the QR decomposition can be computed once during posterior mean inference and reused for the covariance prediction. Thus, the time complexity of posterior covariance inference O(nm²) per test point.
B Alternative Methods for Posterior Inference
Section 3.3 details the adaptation of the QR stabilized linear solve for approximate kernels (Foster et al. 2009) to the SoftKI setting. In this section, we provide empirical evidence that illustrates how other linear solvers fail when confronted with datasets that generate noisy kernels. We focus on two datasets, protein and keggundirected.
As reminder the QR stabilized linear solve motivated by the challenge of solving the linear system:
Ča Kax A y.
(18)
where Ò is the estimated (potentially noisy) kernel matrix Ò Ke + - and R Exc is the interpolated kernel between interpolations points and training points
16

---

Published in Transactions on Machine Learning Research (08/2025)
Dataset
T
Soft m=512
SGPR m=512
SGPR m=512 QR
Pol
13500
26
0.075 + 0.002
0.127 + 0.002
0.121 + 0.002
Elevators
14939
18
0.36 + 0.006
0.395 + 0.005
0.392 + 0.005
Bike
15641
17
0.062 + 0.001
0.108 + 0.004
0.105 + 0.004
Kin40k
36000
8
0.169 + 0.008
0.177 + 0.004
0.169 + 0.004
Protein
41157
9
0.596 = 0.016
0.602 = 0.015
0.601 - 0.015
Keggdirected
43944
20
0.08 0.005
0.099 0.004
0.082 0.006
Slice
48150
385
0.022 = 0.006
0.52 = 0.001
0.513
0.001
Keggundirected
57247
27
0.115 + 0.004
0.128 0.009
0.124 0.007
3droad
391386
3
0.583 + 0.01
Song
270000
90
0.777 + 0.004
Buzz
524925
77
0.24 + 0.001
Houseelectric 1844352
11
0.047 = 0.001
Table Test RMSE on the UCI dataset using vanilla. SGPR with SGPR using the QR stabilized linear solve Reference results for SoftKI given in gray.
described in section 3.3
Method
Learning rate (n)
Kernel
Starting noise (3)
Inducing points (m)
Test RMSE
SGPR
0.1
Matern v = 1.5
0
128 K-means centroids
6 10-3
SKI
0.1
Matern v 1.5
0.1
20 x 20 grid
6 x 10 10-3
SoftKT
0.5
Matern == =1.
0.5
128 K-means centroids
2 x 10-3
Table 8: Experimental settings for the Ricker Wavlet experiment.
In this example we evaluate direct solver, the Cholesky decomposition method, and CG with convergence tolerances set to 10-1, 1 x 10-2, 1 10-3, and 1 x 10-4 Despite adjusting the tolerances for the CG method, experiments revealed that all solvers-except for the QR-stabilized training instability. Figure 3 depicts the test RMSE performance across the different solvers. The direct and Cholesky solvers, while more stable than the CG solvers, failed to produce reliable solutions. Similarly, the CG method did not converge to acceptable solutions within the tested tolerance levels. In contrast, the QR-stabilized solver consistently produced stable and accurate solutions justifying its choice as a correctional measure that stabilizes a SoftKI on difficult problems.
coatine-resulted in
Since the shape of the SoftKT posterior is similar to the SGPR posterior (see Lemma A.2 we also adapt the QR procedure used for SoftKI to SGPR. Table 7 illustrates the results of using QR solving for SGPR. We see that QR solving improves the test RMSE for SGPR. However, unlike SoftKI where QR solving is more important for training stability, we did not observe instability with an inversion SGPR + A using the matrix inversion lemma.
xx
C Additional Experiments
We provide supplemental data for experiments in the main text (Section C.I and Section C.2 We also report comparisons with SKI-based methods (Section C.3) and additional experimental results (Section C.4)
C.1 Ricker Wavelet Experiment
In Section 3.2 we present an illustrative comparison of SGPR, SKI and SoftKI on the task of learning the Ricker wavelet. Figure visualizes the absolute error given as surfaces. Each model was trained for 100 epochs using the Adam optimizer together with StepLR lcarning-rate schoduler. The datasct comprisod of 3. 000 training points and 200 test points. We obtain the initial inducing and interpolation locations for both SGPR and KT by clustering the training inputs with KMeans.
17

---

Published in Transactions on Machine Learning Research (08/2025)
Absolute
10
Figure 4: (Absolute Error Surfaces): Comparison of the absolute the training set domain from
Figure visualized as surfaces.
Dataset
m
Exact GP
SoftKI m = 512
SGPR m = 512
SVGP m 1024
3droad
391386
5.586 = 0.120
8.529 = 0.216
Kin40k
36000
8
0.720 + 0.057
0.453 + 0.012
0.019 + 0.000
0.823 + 0.010
Protein
41157
9
0.612 + 0.015
0.520 + 0.006
0.020 = 0.000
0.904 0.020
Houseelectric
1844352
11
26.066 + 1.335
40.330 +: 1.182
Bike
15641
17
0.112 = 0.001
0.227 0.011
0.013 = 0.000
0.367
Elevators
14939
0.261 0.012
0.178 + 0.004
0.013 = 0.000
Keggdirected
43944
20
3.904 0.208
0.554 0.004
0.021 = 0.000
0.950 0.011
Pol
13500
26
0.340 + 0.023
0.171 0.010
0.012 + 0.001
0.325 + 0.017
Keggundirected
57247
27
4.337 + 0.402
+.72240.004.
0.026 = 0.000
1.217 0.047
Buzz
524925
77
7.107 + 0.053
10.974 + 0.131
Song
270000
90
3.047 + 0.055
5.851 + 0.185
Slice
48150
385
1.363 + 0.004
0.029 - 0.000
1.069 + 0.014
Table 9: Timing per epoch of hyperparameter optimization in seconds.
C.2 Supplemental to Experiments
Hardware details, We run all experiments on single Nvidia RIX 3090 GPU which has 24Gb of VRAM. A GPU with more VRAM support larger datasets. Our machine an Intel 19-10900X CPU 3.70GHz with 10 cores. This primarily affects the timing of SoftKI and SVGP which use batched hyperparameter optimization, and thus, move data on and off the GPU more frequently than SGPR.
Timing. Table g compares the average training time per epoch in seconds to for hyperparameter opti- mization of SoftKI vs SVGP. Note that SGPR is much faster but does not support batched stochastic optimization, and thus, limited smaller datasets. We observe that SVGP and SoftKT have similar performance characteristies due the mini-batch gradient descent approach that both SVGP and SoftKI employ. SVGP requires slightly more compute since it additionally learns the parameters of a variational Gaussian distribution.
C.3 Comparison with SKI-based Methods
As we mentioned in the main text, is not possible to apply SKI to a majority of the datasets our experiments due to dimensionality scaling issucs. Nevertheless, there are some datascts with smaller dimen- sionality that SKI-variants such as SKIP (Gardner et al. 2018b) and Simplex-SKI (Kapoor et al, 2021 can be applied to.
18

---

Published in Transactions on Machine Learning Research (08/2025)
Dataset
Y1
d
Exact GP
SoftKT m=512
SGPR m=512
1=104
Ac- ac-ala3-nhme
76598
126
21.300 + 0.502
1.113 + 0.095
0.035 + 0.000
3.257 = 0.406
Dha
62777
168
36.841 + 0.375
1.195 + 0.017
0.031 + 0.000
2.364 = 0.078
Stachyose
24544
261
15.747 + 0.172
0.534 + 0.011
0.017
0.924 + 0.029
At-at
18000
354
2.669 1.69::0000 0.036
0.365
0.004
0.699 = 0.008
ag
9137
354
97.700 4.388
0.027
0.000
0.368 = 0.023
Buckyball-catcher
5491
444
40.875 3.186
0.175 0.005
0.011 +
0.221 0.007
Da-nanouube
4528
1110
298.096 + 1.270
0.344 0.026
0.012 0.002
0.198 0.002
Table 10: Timing per epoch of hyperparameter oplimization in seconds on MD22.
Dataset
d
Exact GP
SoftKI m=512
SKIP
Simplex-SKI
Protein
41157
9
0.511 0.009
0,632 + 0.012
0,817 0.012
0.571 - 0.003
RMSE
Houseelectric
1844352
11
0.054 0.000
0.052 + 0.0
0.079 0.002
Elevators
14939
18
0.399
0.011
0.423 + 0.011
0.447
0.037
0.510 0.018
Keggdirected
43944
20
0.083
0.001
0.089 + 0.001
0.487 + 0.005
0.095 + 0.002
Protein
41157
9
0.960
0.003
0.992 = 0.017
1.213 + 0.020
1.406 = 0.048
Houseelectric
1844352
11
0.207 + 0.001
0.251 + 0.005
0.756 + 0.075
Elevators
14939
18
0.626 = 0.043
0.372 + 0.004
0.869 + 0.074
1.600 + 0.020
Keggdirected
43944
20
0.838 - 0.031
0.254 + 0.156
0.996 + 0.013
0.797 + 0.031
Table 11: Comparison with SKT-based methods. Numbers for Exact, SKTP and Simplex-SKI are taken from Kapoor al. 2021). Best approximate GP numbers are bolded.
To compare against these methods, we attempt to replicate the experimental settings reported by (Kapoor et al. 2021) so that we can directly compare against their reported numbers. The reason for doing so is because Simplex-SKi relies on custom Cuda kernels for an efficient inplementation. Consequently, it is difficult to replicate on more current hardware and software stack as both the undcrlying GPU architoctures and PyTorch have evolved. From this perspective, we view the higher-level PyTorch and GPyTorch imple- mentation of SoftKI details.
practical strength since it can rely these frameworks to abstract away these
In the original Simplex-SKI experiments, the methods are tested on 4/9, 2/9, and 4/9 splits for training, validation, and testing respectively. For SoftKI, we simply train on 4/9 of the dataset, discard the validation set, and test on the rest. We report the test RMSE and test NLL achieved Table II We find that SoftKI is competitive with Simplex-SKI.
C.4 Lengthscale and Temperature
Figure gives the histograms the lengthscales learned by various methods in the experimention Sec- tion Figure gives the histogram of the temperatures learned by SoftKT. On datasets such as song and slice, we see that the majority of SoftKI's lengthscales for each dimension reach the maximum lengthscale of 5. The dimensiinality-speefff variation is instead captured in the temperature.
19

---

Published in Transactions on Machine Learning Research (08/2025)
20
AVEE
1e5p bafid
P32
25
162
250
(a) 3droad
(b) kin40k
(c) protein
(d) houseelectric
wa3
selt
sofld
33
(e)
bike
(f) elevators
g keggdirected
(h) pol
256
100
13
(1) keggundirected
(i)
buzz
(k) song
(l) slice
Figure 5: Histogram of lengthscale distributions learned by automatic relevance detection (ARD) by different approximate gaussian process methods on the uci dataset.
20
20
2.5
(a) 3droad
(b) kin40k
(c) protein
(d) houseelectric
2,0
035
(e) bike
(f) elevators
(g) keggdirected
(h) pol
20
keggundirected
(i) buzz
(k) song
(I) slice
Figure Histogram of Temperatures T' learned by SoftKI on the uci dataset.
20


Reviews:
[{"review_id": "hH3ogPFDtG", "reviewer": "Reviewer_HaUc", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strengths\n- The paper has a clear goal of allowing kernel interpolation to scale to datasets with higher-dimensional feature vectors. \n- To the best of my knowledge, KISS-like GPs are still the SoTA when it comes to using kernel interpolation in GPs.\n- The method is extensively evaluated and compared.\n\n### Weaknesses\n- The main idea of the paper (the use of the softmax for quadrature weights in (10)) lacks a derivation or even a strong explanation.\n- The paper makes some statements based on intuition; to strengthen the statements and for completeness, these should be supported by a formal argument (proof) or a reference. E.g.:\n  - p.6, \"they have a push-and-pull effect\"\n  -  p.8, \"The stochastic estimate provides an unbiased estimator of the gradient.\"\n- While the method removes the explicit dependence on the feature dimensionality (by removing the grid structure), it is unclear whether SoftKI can utilize the underlying manifold structure postulated for many practical problems. The efficiency could then be explained by, ie, [1].\n- The paper is too empirical. Bayesian inference through GPs is IMO a very useful tool for the ability to derive and explain algorithms. \n\n[1]: Rosa, Paul, et al. \"Posterior contraction rates for Matérn Gaussian processes on Riemannian manifolds.\" Advances in Neural Information Processing Systems 36 (2023): 34087-34121.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The authors develop a variation of the kernel interpolation approach to approximate a GP called SoftKI. SoftKI chooses quadrature points via optimization and assigns the quadrature weights based on a softmax over distances. The method is integrated with standard techniques from variational inference to allow for batched optimization and is tested on well-known datasets.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "### Major Issues\n- The method section largely reads like a laundry list of what needs to be done to make the method work. While those are all important points, more detailed derivations, especially on SoftKI and the stabilized MLL, are IMO in order. In its current form, the paper resembles a successful empirical study. With supporting theory, there might be subsequent work building on it. In particular\n  - \"4.2 Learning Interpolation Points\"\n  - (10)\n  - (17)\n- p.7 \"when numeric instability is encountered\", state the algorithm explicitly\n- The Nystroem method approximates the posterior. We approximate the Kernel with Methods like RFF. Before truncating, those representations are equivalent, see [2], after truncating, they are different and induce different approximation artifacts (variance starvation, loss of resolution).\n\n### Minor Issues\n- Figure 1 needs more interpretation. The distribution of inducing points is explained, the distribution of weights needs more dissemination.\n- The ELBO is never stated in its generic form; it might be good for the exposition to include it in (6).\n- \"2 Background\" and \"3 Related Work\" have no clear distinction.\n- For the Figures with multiple axes, please use subcaptions instead of the legend to label the axes.\n- Table 6, in the last row \"Slice\", is the wrong number bolded? \n- SGPR is still in the table in spite of \"We do not include SGPR\"\n- I cannot follow the line of thinking here: p.18, \"Second, the patterns in the interpolation points learned by SoftKI confirm our intuition that high-dimensional spaces can benefit from fewer interpolation points.\" Are there concrete examples? I cannot spot a qualitative difference.\n- The appendix is not well formatted, eg.: p.19ff trailing text under the figure, missing captions, sometimes the inducing points for SGPR seem to be missing. The captions are, in general, not very descriptive.\n\n### Questions \n- The authors perform an ablation wrt. the number of inducing points. Surprisingly, more interpolation points do not lead to better performance in general. They state that this might be due to overfitting. Does the \n- The paper puts a lot of emphasis on the low computational complexity of single gradient updates. Yet, the algorithm uses dense matrices. I was wondering why not discuss sampling from the SoftKI? It is explicitly interpreted as \na probability distribution over interpolation points for each datapoint. This could sparsify the gradient and robustify the model in a similar fashion to dropout.\n- Plotting inducing points along PC is often used to showcase how SoftKI is \"efficient\". In isolation, I don't know what is expected to be observed in those plots. Is there some intuition I am missing?\n\n\n[2]: Bach, Francis. \"On the equivalence between kernel quadrature rules and random feature expansions.\" Journal of machine learning research 18.21 (2017): 1-38.", "broader_impact_concerns": "I have no concerns about a negative broader impact.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Vq0II7CYn1", "reviewer": "Reviewer_GuQj", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "# Strengths: \nThe paper pursues a new branch in the SKI family (low-rank rather than sparse structured interpolation points), and gets it to work using a combination of contemporary tricks -- it looks like a solid contribution and a new model that should be in the practitioner's toolbox. The paper is well written, notation is clear, and benchmarks are extensive and show the method in favorable light. The new idea (learned softmax interpolation points) is potentially applicable alongside other structured SKI tricks to increase scale further. \n\n# Weaknesses: \n\nI think the paper is good overall -- I have two bigger comments on baselines and potentially confusing exposition, plus some minor points and suggestions. \n\n## The $m \\ll n$ (high $d$) vs $m \\gg n$ (low $d$) setting\n\nI was confused during my initial read of the paper because I kept looking for the part where SoftKI admits fast matmuls (as other SKI methods do) -- but that's not the point here, the \"fast\" matmul is simply matrix inversion lemma on a low-rank interpolation points matrix. My read on the various SKI methods is that they operate in the large-$m$ setting by trading off very good scaling in $m$ against unfavorable scaling in $d$, and attempt to come up with clever fixed interpolation schemes that make the scaling in $d$ less unfavorable. In contrast, the present approach accepts worse scaling in $m$ (by eschewing structure in the interpolation points) and removes the dependence on $d$ (by eschewing fixed interpolation schemes). Is that right? If so, this is all great, but it took me a few reads and going back to the SKI papers to realize this is what was going on, and making it explicit would help a reader (e.g. by reminding us that inducing point approaches often take $m \\ll n$ but SKI approaches can take $m \\approx n$ or even $m \\gg n$, or that it takes a different branch from other SKI approaches). It's all in there in the paper, and once I figured out what's going on I saw it, but better signposts would have helped. One possible idea is just a table with methods and their scaling in $m, n$, and $d$ which would make it clearer that the present contribution is an SKI-like method with SGPR-like scaling. \n\n## Comparing against other methods\n\nSGPR and SVGP are good \"reference\" baselines, and it's also nice to see the favorable behavior relative to Skip and Simplex-SKI in the supplement, but there are other GP methods that do scale to the sizes used in this contribution. For example, VNNGP (which the present paper cites) has linear scaling in both $n$ and $m$ and no explicit dependence on $d$ (though the nearest-neighbor algorithm incurs a one-time cost dependent on $d$), and as such can operate on a n=372K d=54 dataset which is not far off from what is studied here. GPnn (Allison et al. 2023) has similar scaling and tests on datasets with $d$ up to 378. CaGP (Wenger et al. 2024) has favorable scaling in both $m$ and $d$ and runs experiments on datasets with n up to 1.8M. I don't expect the present contribution to compare against all of them (and under the TMLR rubric, I don't think comparison against *any* of them is needed to make the point that the present contribution is correct and interesting), but it would strengthen the paper, particularly because my sense is that SoftKI may implicitly find a set of neighbors for each point during optimization. It would also make the comparison to where the interpolation points are (section 5.3) potentially more interesting. \n\n## Minor points\n* SKI vs KISS-GP: Wilson and Nickisch make a distinction between SKI (the approximation given in expr [7] in the present paper) and KISS-GP (SKI with cubic interpolation and with structure-exploiting linear algebra). The review in section 2.3 does not make that distinction, and only reports the complexity of the Toeplitz approach, whereas the KISS-GP paper also discusses Kronecker. \n* It's true that naive GP inference is $O(n^3)$ complexity (as claimed in the introduction) but I think by now the fast MVM iterative solver approach has been around long enough to where we can claim $O(n^2)$ is the baseline for modestly large $n$. I don't think there's a need to background this under \"efforts have been made to accelerate exact GP inference\" especially considering that those tricks underpin the standard SKI approach. \n* The notation for a Gaussian distribution ($\\mathcal{N}(\\cdot \\mid \\mu, \\mathbf{\\Sigma})$) is defined after expression 2 even though it is used in the model likelihood and expr 1 above. So it should possibly be moved up, though I'm also not sure that it needs to be defined for an ML audience. \n* Fig 1 caption should explicitly describe what the magenta segments are (they are connecting the points used for interpolation?). \n* I'm not sure what to make of Fig 3 and section 5.3. It looks like the distribution of SoftKI points doesn't overlap as much with the data over the first two components as SGPR / SVGP, the conclusion isn't much beyond \"they're different\" and the discussion is vague. I think it could easily be moved to the appendix in favor of bringing up the Skip and Simplex-SKI results.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "The contribution presents a new method in the SKI family where rather than exploiting sparsity and structure in the interpolation matrices of points selected a priori, it learns the interpolation points using an approximate sparsity constraint. While this brings back quadratic dependence on the number of interpolation points, it removes the dependence on data dimensionality, which allows the method to scale to higher-dimensional datasets. Further, the fact that interpolation points are learned lets the approach perform well on these datasets even without requiring very many interpolation points. A number of additional techniques are used to make inference stable and tractable.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "See discussion above -- I think my main requested change is addressing the settings of high $m$ low $d$ vs low $m$ high $d$ and trying to situate the new contribution in a way that is less confusing for readers. I also strongly encourage addressing all the minor points. I think an additional baseline or two against contemporary methods with favorable scaling would also be beneficial but I don't think it's critical, especially given the TMLR rubric.", "broader_impact_concerns": "No concerns.", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Xswt1j8a3h", "reviewer": "Reviewer_GRfH", "summary": "", "questions": "", "limitations": "", "rating": null, "confidence": null, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strengths:\n\n-\tI appreciate Sections C.3 and C.4 as I understand that the complexity of the algorithms has made some amount of exact reproduction impossible over the course of 3+ years.\n\n-\tI like the molecule experiments and the heavy interpretation of them, which are a nice change from just the traditional UCI experiments.\n\n-      Overall, the technique is pretty straightforward and well explained. Good job.\n\nWeaknesses:\n\n-\tWhile I appreciate all of the numerical advances in this paper, I strongly suspect that the reason why the numerics are such a pain under the hood is the fundamental degeneracy of the model itself. The subset of regressor approximations can be viewed as finite rank GPs or Bayesian linear models with |z| features (see Yadav et al, ’21 [https://proceedings.mlr.press/v130/yadav21a/yadav21a.pdf] for a modern description).\n\n-\tThere’s not a lot of comparisons to exact GPs on the datasets where they can scale reasonably well (see below).", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "Yes", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "paper_topic_and_main_contributions": "", "excitement": null, "reproducibility": null, "ethical_concerns": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "summary_of_contributions": "This paper extends the structured kernel interpolation for GP kernels by replacing the sparse matrix W with a softmax matrix and then learning the location of the inducing points. As the resulting GP is low rank, the numerics of the GP matter so the authors use the Hutchinson’s pseudo-loss (from Maddox et al, ’22, Wenger et al, ’22) when the gradients are unstable. For predictive means and covariances, due to the degeneracy, the authors also use a QR decomposition for the solve instead of conjugate gradients.", "claims_and_evidence_yn": "Yes", "claims_explanation": "", "audience_yn": "Yes", "audience_explanation": "", "requested_changes": "Questions:\n\n-\tHow is the change in the value of the loss dealt with on the optimizer side when the loss used flips to the Hutchinson’s pseudo loss? \n\n-\tHave the authors tried directly optimizing the pseudo loss for the entire trajectory (like in Maddox et al, ’22)? As the gradients are the same then the optimization results should be identical, up to numerical differences?\n\n-\tSKIP and then Pleiss et al, ’18 (https://arxiv.org/abs/1803.06058) give efficient methods for predictive variances (more so than the algorithm described in Appendix A). Could these methods be applied to speeding up predictive variance computations (e.g. for NLLs) here?\n\n-\tTables 1, 2, 3 and 4: why not include comparisons with exact GPs here (on at least the <50k data sets)? For Tables 1 and 2, I believe that these are basically the same experiments as in Wang et al, ’20, Gardner et al, ’19, and Maddox et al, ’22 (for the floating point experiments).\n\n-\tIn some cases, the pivoted Cholesky pivots (essentially) (Burt et al, ’20, JMLR[ https://www.jmlr.org/papers/v21/19-1015.html]) are better at initializing GPs in the inducing point setting, did the authors try that instead of k means?\n\n-\tFollowing up from the connection with Yadav et al, ’21, is there a more numerically stable way to compute the predictive mean using the Bayesian linear regression formulation?\n\n-\tAs this approximation ends up being really subset of regressors (Foster et al, ’09; Smola & Bartlett, ’00), why not compare to just that technique with both fixed and learned inducing points?", "broader_impact_concerns": "n/a", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)