Scaling Gaussian Process Regression with Full Derivative Observations
Authors: Daniel Huang
Venue: Accepted by TMLR ()

Abstract: We present a scalable Gaussian Process (GP) method called DSoftKI that can fit and predict full derivative observations. It extends SoftKI, a method that approximates a kernel via softmax interpolation, to the setting with derivatives. DSoftKI enhances SoftKI's interpolation scheme by replacing its global temperature vector with local temperature vectors associated with each interpolation point. This modification allows the model to encode local directional sensitivity, enabling the construction of a scalable approximate kernel, including its first and second-order derivatives, through interpolation. Moreover, the interpolation scheme eliminates the need for kernel derivatives, facilitating extensions such as Deep Kernel Learning (DKL). We evaluate DSoftKI on synthetic benchmarks, a toy n-body physics simulation, standard regression datasets with synthetic gradients, and high-dimensional molecular force field prediction (100-1000 dimensions). Our results demonstrate that DSoftKI is accurate and scales to larger datasets with full derivative observations than previously possible.

Published in Transactions on Machine Learning Research (01/2026)
Scaling Gaussian Process Regression with Full Derivative Observations
Daniel Huang Base26, CA, USA
dan@base26labs..om com
Reviewed on OpenReview: https: // openreview. net/ fttl/form? id- fbonXp38r9
Abstract
We present a scalable Gaussian Process (GP) method called DSoftKI that can fit and predict full derivative observations. Tt extends SoftKT, a method that approximates kernel via soft- max interpolation, to the setting with derivatives. DSoftKI enhances SoftKI's interpolation scheme by replacing its global temperature vector with local temperature vectors associated with each interpolation point. This modification allows the model to encode local directional sensitivity, enabling the construction of scalable approximate kernel, including its first and second-order derivatives, through interpolation. Moreover, the interpolation scheme elimi- nates the need for kernel derivatives, facilitating extensions such as Deep Kernel Learning (DKL).. We evaluate DSoftKT on synthetic benchmarks, toy n-body physics simulation, standard regression datasets with synthetic gradients, and high-dimensional molecular force field prediction (100-1000 dimensions). Our results demonstrate that DSoftKT is accurate and scales to larger datasets with full derivative observations than previously possible.
1 Introduction
A convenient feature of using Gaussian Process (GP) to approximate functions is its ability to incorporate derivative observations since the derivative of GP is also a GP. This enables more informative modeling when derivative data is available (e.g., in the physical sciences). However, utilizing GP with derivative observations, abbreviated GPwD, for regression faces significant scalability challenges. In particular, vanilla GPwD inference has time complexity 3d)3d3) where n the mumber of data points in a dataset and d is the dimensionality of the data. Thus, GPwDs have scaling challenges in both n and
To alleviate these challenges, scalable GP regression methods such as Stochastic Variational GP (SVGP) (Hensman et al., 2013) have been extended to the setting with derivatives (DSVGP) (Padidar et al. mani
2021). It uses < inducing points C'ninonero-Canded & Rasmussen 2005 Snelson & Ghahra-
2005 and achieves a time complexity of (mmdd²) for posterior inference. Since the cubic scaling in d can be prohibitive, DSVGP with directional derivatives (DDSVGP) (Padidar et al. 2021) introduces < d inducing directions, an analogue of inducing points for dimensions, to achieve a time complexity of O for posterior inference However, this comes the cost of directly predicting derivatives and introducing further approximations.
In this paper, present a GPwD method called DSoft that can fit and predict all derivative information (Section 4. It is extension of SoftKI (Camaño & Huang 2025). scalable GP method that approxinates kernel via softmax interpolation from m & n interpolation points whose locations are Iearned. Thus, it blends aspects of kernel interpolation (Wilson & Nickisch, 2015) which introduces kernel interpolation from a structured lattice with variational inducing point methods (Titsias 2009 where inducing points are adaptod to the datasct. To handle gradient information, DS modifies SoftKT's interpolation schome by replacing its global temperature vector with local temperature vectors associated with each interpolation
This assumes mp² d which no longer holds around 2000 for 500 and p 2. See Appendix for details.
2Cod available at httpp://gitthbl wow.baase2llass/teofff.ooo

---

Published in Transactions on Machine Learning Research (01/2020)
Branin
Vani la
DSVGP
DDSVGP
DSoftKi (Ours)
Valuo
Figure Comparison GPwD regression with vanilla GPwD, DSVGP, DDSVGP, and DSoftKI on the Branin surface (2D). We plot the contours the surface, the gradient with respect the first argument (V1), and the gradient with respect to the second argument (V2). Vanilla GPwD is accurate but intractable for sizable and/or d. DSVGP forms nice approximation of the original surface but is intractable for large d. DDSVGP uses p 2 directional derivatives to enhance scalability but loses fidelity in modeling the surface. DS provides an accurate and scalable approximation of the surface.
point. This allows the model to encode local directional sensitivity (Section 4.1 enabling the construction of a scalable approximate kernel, including its first and second-order derivatives, through interpolation. Notably, this design eliminates the need for kernel derivatives, facilitating extensions such as Deep Kernel Learning (DKL) (Wilson al. (0166) which use learned kernels (see Appendix D This contrasts with methods such DSVGP DSVGP/DDSVP. which introduce separate inducing points for values and gradients well require computing kernel derivatives, leading to scaling challenges in both n and (see Appendix D for more details). As a result, DT achieves time complexity of ((m2nd) for posterior inference. Because each datapoint consists of d values to fit, we can view KI as achieving similar complexity to an approximate GP (e.g., Titsias 2009)) that is fit to r2d+ + 1) datapoints.
We evaluate DSoftKI on synthetic functions with known derivatives as a bascline compare against existing GP regression and GPwD regression methods. Since KI enables derivative modeling with GPs at larger and than previously possible, also evaluate its efficacy high-dimensional molecular force field modeling task (d around 100-10000 which requires predicting gradients. To test broader range of datasets, we also validate performance on a toy n-body physics simulation and standard UCI regression benchmarks with synthetic derivatives (Appendix C Our experiments show that DSoftKT is a promising method for GPwD regression with full derivative prediction, both in terms of accuracy (e.g., see Figure I and scaling to larger n and d than previously possible (Section
2 Related Work
One strategy for scaling GPwD regression to extend an existing scalable GP method to the setting with derivatives. In this vein, both Stochastic Variational Gaussian Process (SVGP) (Hensman al., 2013) and Structured Kernel Interpolation (SKT) WWilson & Nickisch 2015) have been extended to the setting with
2

---

Published in Transactions on Machine Learning Research (01/2026)
Method
Predict V? Kernel V?
Posterior Inference
Vanilla
yes
yes
dn)
DSVGP
yes yes (m))
DDSVGP
no yes
(m)
DSoftKI (ours)
yes
no
O(m2nd)
Table A comparison of GPwD methods across several dimensions. As a reminder, n the number of points, d is the dimensionality of the data, m is the number of inducing/ internnggtterroolation points, and p the number of inducing directions. The proposed method, DSoftKI, scalable, predicts full derivative observations, and does not require computing first or second-order derivatives of the kernel to construct the GPwD kernel. The kernel approximations made by each method are discussed in more detail in Appendix A.1
derivatives. Since we also follow this strategy, we briefly summarize these methods here and refer the reader to Appendix E for more background on these related works.
DSVGP Padidar et al. 2021 extends SVGP to the setting with derivatives by modifying the SVGP varia-
tional approximation to the GPwD kernel, resulting in a kernel matrix of size m(d + 1) x m(d + 1) where m. is the number of inducing points (Snelson & Ghahramani 2005; Quinonero Quinonero-Candela & Rasmussen 2005
Like SVGP, the inducing points in DSVGP are learned by optimizing an evidence lower bound (ELBO) that can be computed with stochastic variational inference. The time complexity of computing the ELBO for each miníbatch of optimization has time complexity of Om)?d3) and the time complexity of posterior inference is dd) Since this can be prohibitive for large d, DDSVGP (Padidar et al., 2021) utilizes < d directional derivatives to further improve the time complexity of computing the ELBO per minibalch and posterior inference to O(m²p") since the d-dimensional gradients have been reduced to p-dimensional directional dcrivatives.
DSKT (Eriksson et al.j 2018 extends SKT (Wilson & Nickisch, 2015) the seiting with derivatives by approximating the gradients of the SKI interpolation kernel. Like SKI, the resulting DSKI kernel has structure that enables fast matrix-vector multiplications (MVMs), and consoquently, conjugate gradient (CG) methods to perform GP inference. Unlike SKI which uses a cubic interpolation scheme, DSKI uses quintic interpolation scheme to better handle gradient observations, resulting in time complexity of a single MVM of (nd66' + m l). Thus, the scaling in the dimension d is even worse than SKI, limiting the application of DSKI to even smaller dimensions. There are other kernel interpolation methods that improve the dimensionality scaling (e.g., see Kapoor et al. (2021 0: Yadav et al. (2023) that can be used besides SKT. However, to the best of our knowledge, these works have not been extended to the setting with derivatives.
Beyond the two works above, there has been relatively little additional work done in GPwD regression compared to GP regression. De Roos et al. (2021) tackles GPwD regression in high dimensions but is limited to low n settings. GPyTorch (Gardner et al. 2021 opens the possibility for scalable GP inference on GPUs and provides support for many standard kernels by hard-coding their first and second-order derivatives but does not contain a standard implementation of a scalable GPwD. GP regressions have been scaled with CG solvers [Wang et al., 2019) on multi-GPU hardware and stochastic gradient descent (Lin et al.) 2023). However, to the best of our knowledge, this has not been studied in the setting with derivatives. Our work makes significant improvements in the size of n and d that can be handled in GPwD regression with full derivative observations.
3 Background
We also review SoftKI since our method extends it to the setting with derivatives (Section 3.3 We begin by introducing notation that will be used throughout this paper.
In this section, we review background on GPs (Section 3.1 and GPwDs (Section 3.2
Notation. The notation A [)(i defines a matrix whose i-th and j-th entry Auj g(i,)) for some function g defined indices and j. Given list of vectors (ar Rdy indexed by i, define the vector
3

---

Published in Transactions on Machine Learning Research (01/2026)
x == (zdla (i.e., g(i. xa) to be the flattened length nd column vector. Similarly, x [my]1 is the equivalent row vector. The notation f(x) [f(=)]aa maps a function f over flattened vector x. Given lists of vectors (x, E [== and (z' E )==) j=1. define the matrix Kxx' = [6(21,)))6 so that it maps function k over the pair of x [2,/s1 and x' [all.
3.1 Gaussian Processes
A Gaussian process (GP) is a random (continuous) function f Rd - R". It is defined by a mean function Je: Rdd + IR" and a positive semi-definite function TRd x IR + R(xx)) called a kernel function. A GP has Gaussian fiiiee-dimensional distributions so that f(x) ~ N(xx, Kxx) for any x where M(xx, Kxx) indicates a Gaussian distribution with mean jxx = =((x) and covariance mairiz Kxx = [2)1s, Without loss of
generality, we will assume Px 0 since we can shift the mean of Gaussian.
To perform GP regression on the labeled dataset {( Ut) 24 CR, yi ER].., use the generative process
f(x) Kxx) y fl ~ ~(((x)),A)
(GP)
(likelihood)
where f GP and y is [xx perturbed by Gaussian noise with covariance A = B8 The noise 3 is example of GP hyperparameter. Others include the kernel lengthscale E and scale y. The hyperparameters can be set by maximizing the marginal log likelihood (MLL) of a GP
logply xx:0 -N 0, Kxx (1,7) A(B))
where we have explicitly indicated the dependence of Kxx on l and 7, and A on 8 for hyperparameters 0 (,). The time complexity of computing the MLL is O
Once the hyperparameters have been set, we can perform posterior inference. The posterior predictive distribution has the closed-form solution (Rasmussen & Williams 2005
= (Kx + A) K.. - K.x + A) "xx))
where ( I), 1) is notation for the probability density function (pdf) of a Gaussian distribution with mean a and covariance P. The time complexity of posterior inference O(n²) which the complexity of solving the system of linear equations in n variables (Kxx A)a = for a, i.e., a = (Kxx + A)-
3.2 Gaussian Processes with Derivative Information
If f Rd + R is a GP with kernel k Rd x gd + R, then Vf Rd + IRd is also GP with kernel k/ Rd x Rdd - (d(xx) defined as
Ok y) dz,0y; 1.j
(x
vzk(z,))
(3)
Consequently, we can construct GP with derivative observations (GPwD), random vector-valued function i R. - [Rd+1 defined as
f(x) = vf(x)) (1(ee)
f(x)
(4)
that simultaneously models a function ff and its gradient Vf. It has a jointly Gaussian distribution f(x) N(O,, Kxx) where Kyx (e(x1,,;))e and
((10))-((aa k(2xx') x²) (om(ena") la
a lis Fk(e, (on
k(2,2')
= -V,,) k(o, x') (1 V.)
(5)
where I is the identity operator that Kxx is nld + 1) ned + matrix.

---

Published in Transactions on Machine Learning Research (01/2026)
To perform GPwD regression the labeled dataset {(r,, ye, dy:) as E R Ms R, R where each dy: is gradient label, we use the generative process
f(x) ~ AO0 Kxx)
(GPwD)
y flf(x) ~ (
(likelihood)
where y lys dy. n(d+ 1) 1 vector of values and gradient labels and
A 1000 0 )
(6)
is a n(d 1) x n(d+1] diagonal matrix of noises, 83 for the function value and &g for the gradients. GPwD's hyperparameters can also be set by maximizing the MLL
A
log x; = Ny 0, Kxx + A(B,,
(7)
where 8 8= (e, Be. (g). The time complexity of computing the MLL d())d²)
The posterior predictive distribution
p(((+)
+ +A) A K.x A) Kxx.)
(8)
which takes the same form GP regression where replace the corresponding variable with its version. The complexity of posterior inference is O(1ad which is the complexity of solving system of linear equations in n(d+ 1) variables. Notably, GPwD regression has an asymptotic dependence on the data dimensionality d since each partial derivative contributes an additional linear equation.
3.3 Soft Kernel Interpolation
Soft kernel interpolation (SoftKT) Camañño & Huang. 2025 is an approximate GP method that uses an approximate kernel
SoftKI xx
Exx K zn Enx
(9)
where Exz= = [ob(x:))is
(x(x) = exp (- 1 2 T 2;
(10)
Σm k== exp exp (- 2 2x
performs softmax interpolation between m & n interpolation points z whose locations are learned, Ø is a Hadamard division (f.e., element-wiee division), and T Rd a learnable temperature vector akin to using automatic relevance detection (ARD) (MacKay et al., 1994) to set lengthscales for different dimensions. Thus, it is method that combines aspects of kernel interpolation from SKT and adaptability of inducing point locations to data in variational inducing points method to obtain scalable GP method.
Although the interpolation points bear resemblance to inducing points, they are not associated with cor- responding inducing variables that have normal priors which are introduced for the purposes of variational optimization. Instead, interpolation points are selected based on distances to data, and so encode the geom- etry of the underlying data that is most useful for interpolation (Figure 2 They are learned by optimizing combination of the SoftKI MLL, and an approximate MLL when numeric instability arises due the current placement of interpolation points. More concretely, the objective is
ply
x; 0) when numerically stable
logIly x; 0)
(11)
logply
x)
otherwise
where log p(y x; the SoftKI MLL and log p(y pseudoloss [Maddox al. 2022). is defined as
x; 0) is an approximate MLL termed Hutchinson's
log x)
Ug Doun
Σ 1 (Dow,) ==1
(12)
5

---

Published in Transactions on Machine Learning Research (01/2026)
Branin
Inducing Points (SVGP)
Interpolation Points (SoftKI)
Value
Figure 2: Comparison of learned interpolation points locations versus learned inducing point locations (only those unit square shown) on the Branin surface. Interpolation point locations encode geometric structure in the data that is most useful for interpolation whereas inducing point locations reflect the normal priors placed their associated inducing variables.
where Dg K (L, 7, 2) A3),, un ur are solutions to the equations Du un ... ur... ui)
(yw1 ... wi), and wy for 1 j < I are Gaussian random vectors normalized to have unit length. Ob- serve that
ÖD, a0
aDg
Ve logply x;
Ug
ug
+ tr Dø
(13)
00
7D,
1 Σ
7Dg
Un
Uo +
00 W
(14)
00
j=l
Ve logāly x)
(15)
that the gradient of the Hutchinson's 's pseudoloss approximately equal to the gradient of the MLL when
approximated with Hutchinson's stochastic trace estimator [Girard 1989 Hutchinson 1989.. The time complexity of computing the Hutchinson's 's pseudoloss per minibatch of size b is O + m²)
The posterior predictive distribution
-- SoftkI
RRxAA
- C- K Sofki)
(16)
where Rxn = Exx= and C Kzn + R A xx- It is similar the posterior predictive of Sparse Gaussian Process Regression (SGPR) (Tissias 2009), the difference being that we replace kernels in the SGPR posterior with the interpolated versions. The time complexity of posterior inference is O(m²n).
4 Method
In this section, extend SofKKI to work with derivatives. First, we introduce the DSoftKT kernel (Sec- Next, we discuss posterior inference (Section 4.2
tion 4.1
Finally, we discuss the role of value and gradient
noise hyperparameters in DSoftKI (Section 4.3
4.1
Soft Kernel Interpolation with Derivatives
The DSoftKl kernel takes the same form DSKI, i.e.,
DSoftKI xx
Exe Kzg Enx so Kxx
(17)
where
Σxe = [ (Go((eo)) ij
= ((TT,) of(zz)
(18)

---

Published in Transactions on Machine Learning Research (01/2026)
Styblinski Tang
Shared T
Individual T
Value
Figure 3: Reconstructed Styblinski Tang surface using shared temperature (original SoftKI scheme) versus individual temperature (proposed scheme) during interpolation for DSoftKI. We also overlay the learned interpolation points (scaled and translated to fit the unit square) in green to illustrate the differences in learned interpolation points.
Unpacking the definition, we obtain
Oig (20011)0015) = ((000) 12 (0a,,) (Kz((00s)))
Jij JijKardi
DDSofKK xx
(19)
where we have abbreviated d3(x;) oij and VoL(x;) doy reduce clutter. Notably, the computation of the DSoftKI kernel does not require computing the first-order second-order derivatives of the kernel. Instead, it is approximated via interpolation. Consequently, it is tractable to perform hyperparameter optimization with first-order gradient methods. Additionally, it can support a wide range of kernels, including learned kernels as DKL. This contrasts with methods such as DDSVGP where the kernel and its gradients are hard-coded in practice alleviate the computational cost and to enable tractable hyperparameter optimization with gradient-basse methods. We refer the reader to Appendix D for demonstration of DKL with DS.ftKI.
While the DSoftKT kernel affords computational savings, it also increases the burden of kernel interpolation to approximate kernel's first and second-order derivatives. To solve this challenge, we adapt the softmax interpolation schome to additionally account for directional orientation of interpolation points relative to the data. More concretely, we replace the global temperature vector with local temperature vectors associated with each interpolation point. this way, interpolation points that are close each other can nevertheless model variations surface curvature as influenced by relevant gradient observations. The original softmax interpolation scheme would be unable to distinguish these variations since the global temperature vector enforces a shared orientation for every interpolation point.
To implement the above idea, we associate each interpolation point z with a corresponding learnable temperature vector Tş E Rd for 1 1 k < m as in
cxp cxp | T - 21
x<(x))
(20)
Th k= exp (- || 3 Th 2x
This introduces x m extra parameters (dxm) (Ti ... T..) the DSoftKI model that can be learned via hyperparameter optimization. The gradient of the softmax interpolation has closed-form solution
V Vz,) = Σ ol(:)))05f of (zx))
It Ø Tí 2
0 T
(21)
||, 2 Tx
2
k=1
which can be obtained in the standard way via the chain-rule. The j-th component of the term
24 TT =4
Tzj Tkj
2 Th
(22)
Th 2z +
T2 2 Th 2 +

---

Published in Transactions on Machine Learning Research (01/2026)
Algorithm 1 DSoftKI regression adapts the SoftKI algorithm to handle the additional computational challenges of the DSoftKI kernel and learning local temperature vectors T(mxd).
Require: DS hyperparameters (6, (1,,,2//xx T(mxd), 3u, 9g)
Require: Dataset (x.
Require: Optimization hyperparameters:: batch size b, number of epochs E, and learning rate 1).
Ensure: Learned DSoftKI coefficients a.
1: for i 1= 1 to E do
for Xs, yt in batch((x, ý), b) do
P batch splits the dataset into chunks of size b
3:
L, L +- Cholesky(K,)
4:
F + ExxL
Low rank representation of K (4,7, T) FF"
D3 [FFT + (B. 3y)]l=
I delays computation of result until needed Stabilized DSoftKT MLL (Equation 25
6:
0 x-- 0 + 1/Vo A.ologil, P(ys Xb; Do)
7: end for
8: end for
9: a f solve Ča = R 'y)
> Solve system of linear equations (Appendix A.2
10: return
where have added small factor e > 0 to the denominator to avoid division by zero, is proportional to the direction from the scaled interpolation point T& @ z to the data point 2; along dimension j, with scaling factor 1/T As a result, points that are nearby can nevertheless have different influences on the interpolation strength of the gradients depending on their respective dircctions. Figure 3 illustrates the differences in surface reconstruction using a shared versus individual temperature vector across interpolation points. We refer the reader to Appendix B.5 for further discussion on the role of temperatures and their connection with lengthscales.
4.2 Posterior Inference
Algorithm summarizes the adaptation of the SoftKI algorithm to work with the DT kernel. Since DSoftKT introduces different interpolation scheme and more learnable parameters, this introduces different optimization dynamics and computational concerns.
Hyperparameter optimization. We learn the locations of the interpolation points by using stochastic gradient descent on the stabilized DSoftKI MLL. The MLL of DSoftKI is
logp(ý/>; 0) (ý 0.)
(23)
where D3 K DSoftKl (f, 7,,7,2, T) + A(B,, B,) for 0 = (e, 7, T, Be, Bg). We refer the reader to Appendix B.1 xx
for details on hyperparameter initialization. For a minibatch of size b, the resulting matrix is of size b + 1) x b(d + 1) as opposed to b x in KI. For sizable and d, this would be intractable to compute with since it requires solving a system of linear equations.
As result, we decompose DSofkK FFT where Kzz LL I is a Cholesky decomposition and F ExxL. This changes the space requirement from O(b2d?) for direct representation to O(bdd). In practice, we expect < bd, since m and b are around the same order of magnitude. Since the Cholesky decomposition is necessary to retain a tractable representation of the DSoftKT kernel, we perform it in double-precision when we encounter numerical instability. In our experience, we rarely encounter numerical instability due to the current placement of interpolation points in DSoftKI. We believe that this is the case since interpolation points can additionally take into account directional information DSoftKT so that points that are close by distanco-wise can nevertheless be oriented to prioritize different directions basod on gradient observations.
The factored representation is used to compute the DSoftKI MLL using lazy representation
D3 [FFT" + (3.,,,)1
(24)
that directly stores F and A(Se, and I- delays computation of any intermediate result until needed. To take advantage of this representation, use low-rank multivariate Gaussian distribution which
8

---

Published in Transactions on Machine Learning Research (01/2026)
d
SoftKT
SVGP
DSVGP
DDSVGP
DSoftKT
Branin
2
0.004 + 0.0
0.018 + 0.004
0.088 + 0.062
0.176 + 0.031
0.003 + 0.001
Six-hump-camel
2
0.026 + 0.003
0.05 + 0.003
0.101 + 0.031
0.669 + 0.04
0.015 + 0.006
Stybliinski-tang
2
0.025 + 0.001
0.05 + 0.003
0.101 + 0.01
0.125 + 0.06
0.012 + 0.002
Hartmann
6
0.05 + 0.002
0.164 + 0.006
0.335 + 0.009
0.346 = 0.019
0.011 + 0.004
Welch
20 0.01 + nan
0.065 = 0.001
0.578 = 0.013
0.003 + 0.003
Table 2: Test RMSE (best bolded) on selected synthetic datasets. One of the runs for SoftKI encountered mumerical instability (hence +nan). indicates a timcout for DSVGP.
avoid working with the full bid + 1) x b(d + 1) matrix Do via the Woodbury matrix identity and matrix determinant lemma. This reduces the challenge of computing the MLL to computing the determinant and inverses of the m x m capacitance matrix F+ A (,,,,,,) to form. When numerical instability is encountered during the computation of DSoftKI's MLL, we compute Hutchinson's 's pscudoloss. The time complexity of single MVM is O(mbd). The stabilizod DSoftKI MLL is thus
instead which takes time complexity O(m²bd)
log p( x; when numerically stable
log log 0)
(25)
log p(y x; otherwise
B.2 contains an ablation of DSoftKI using the DSoftKI MLL compared to the Hutchinson's 's pseudoloss.
Appendix
Posterior inference. Once have learned the interpolation points, we construct the posterior pre- dictive distribution. It is SoftKI's posterior, with the corresponding variables replaced with the (.) versions
((()))8)=((Re)) -RRKKKKx RzxA
-'A--1 -ly, DDDoftKt DSSfft - Redd - ")
(26)
where Rxz= = Exe and Ć == Kez + K -Kx Since each data point in GPWD regression introduces d + 1 values to fit, we have a system of n(d + 1) equations in m variables. Consequently, solving the system of linear equations takes ((ndm²) time and space. For large d, this can exceed GPU memory limits. Consoquently, we sometimes solve these equations on a CPU. It would be an interesting direction of future work to see how this can be improved such as with alternative linear solvers or utilizing multi-GPU hardware.
4.3 The Role of Value and Gradient Noises
Since DSoftKT forms its approximate kernel and its derivatives via interpolation, the choice of value and gradient noises become intertwined. To examine this further, we can unpack the posterior mean equation and see that
Ča
Σ
k(zar Zz)Waz
(27)
where
n
1
wj
32 Jist +
(Dor,) )" )ddy dy:
(28)
32
for 1 a < m. This indicates that cach posterior cocfficient jointly influences the reconstruction of function values and their gradients as opposed to introducing separate posterior coefficients for values and gradients. The influence of fitting values versus gradients on the weights is determined by the ratio $2/32. A ratio of /3 = d suggests that values and gradients are equally weighted because each gradient has d components. A ratio larger than d thus prioritizes gradients while ratio smaller than prioritizes values. We investigate the impact of the ratio the value and gradient test RMSE more in Appendix B.3
9

---

Published in Transactions on Machine Learning Research (01/2026)
SoftKT
SVGP
DSVGP
DDSVGP
DSoftKI
Branin
-0.725 + 0.017
0.116 0.0
-1.044 + 0.885
-0.532 + 0.068
-4.432 + 0.162
Six-hump-camel
-0.743 + 0.057
0.123 + 0.001
-1.539 + 0.159
-0.64 = 0.058
-2.175 + 0.237
Styblinski-tang
-0.54 + 0.082
0.123 + 0.001
-1.182 + 0.166
-1.064 + 0.188
-2.296 + 0.045
Hartmann
-0,494 + 0.111
0.221 0.001
-1.362 + 0.029
-1.324 + 0.027
-3.16 + 0.325
Welch
-0.699 + nan 0.141 = 0.0
0.706 = 0.055
-1.14 + 4.795
Table 3: Test NLL (best bolded) selected synthetic datasets. One of the runs for SoftKI encountered mumerical instability (hence +nan). indicates a timcout for DSVGP.
d
SoftKI
SVGP
DSVGP
DDSVGP
DSoftKI
Branin
0,194 + 0.005
0.339 + 0.149
1.318 0.02
1.945 + 0.022
2.555 0.199
Six-hump-camel
0.257 + 0.092
0.373 + 0.104
1.31 0.008
1.955 + 0.062
2.434 0.087
Styblinski-tauy
0.21 = 0.006
0.335 = 0.14
1.316 0.018
1.91 = 0.016
2.539 0.137
Hartmann
0.199 + 0.003
0.337 + 0.126
9.305 + 0.113
1.934 + 0.045
2.544 + 0.091
Welch
20
0.213 + 0.002
0.418 + 0.123
1.944 + 0.063
4.452 + 0.192
Table 4: Wall-clock training time seconds per epoch (best GPwD bolded) on selected synthetic dataset. indicates timeout for DSVGP.
5 Experiments
In this section, evaluate DSoftKI synthetic functions (Section 5.1 to obtain baseline of comparison and high-dimensional molecular force field (Section 5.2) to test scale. We refer the reader to Appendix C for more experiments with toy n-body simulation and on the UCI datasct (Dua & Graff 2017) with synthetic
gradients.
Baseline GP and GPwD methods. We use the default GPyTorch implementations GGardner et al. 2018 of SVGP, DSVGP, and DDSVGP with 2 inducing directions as baseline variational GPs. We use the PLL
Jankowiak et al.) 2020) modification the ELBO for the variational GPs as recommended by (Padidar 2021b to allow separate noise parameters for function values (B?) and gradients (82). We also use SoftKT baseline GP where also adopt the DSoftKI interpolation scheme. We do not compare against DSKI since the dimensionalitt of the datasets are too high. We use the RBF kernel (with scale) using automatic relevance determination (ARD) lengthscales with the exception of DDSVGP which does not support it. Unless otherwise stated, we use m = 512 inducing points for all mothods. We optimize all hyperparamotors using the Adam (Kingma & Ba 2014) optimizer. All GPs and GPwDs use single-precision floating point numbers, except for DSoftKI which uses Couble-precision occasionally enhance the stability of Cholesky decomposition described previously.
et al.
5.1 Regression with and without Derivative Information
Our first experiment tests DSoftKT on selected synthetic test functions ranging from d= d = 20 with known derivatives following (Padidar et al. 2021) so that we can compare against existing GP regression and GPwD regression. For each function, we generate 20000 datapoints, using 10000 points for training and reserving 10000 points for testing.
Effective learning rate. Since the training data that each GP observes varies across methods due to how each method utilizes derivative information, setting up a fair comparison between each method is not straightforward. For instance, even comparing the variational methods SVGP, DSVGP, and DDSVGP is nuanced since the amount of training data cach method encounters is different at n, n(d + 1), and n(p + 1) where p the number of inducing directions respectively. To control for this, we choose to use the notion of effective learning rale Aef defined Aeff DArase where D is a method-specific number of derivative
10

---

Published in Transactions on Machine Learning Research (01/2026)
DSoftKT RMSE
DSoftKI NLL
A RMSE
A NLL
Branin
0.002 + 0.0
-4.266 + 0.009
-0.000
0.166
Six-hump-camel
0.163 = 0.02
-0.392 + 0.112
0.148
1.783
Styblinski-tang
0.026 + 0.002
-1.829 + 0.171
0.014
0.467
Hartmann
0.336 + 0.034
0.336 = 0.092
0.326
3.495
Welch
0.027 = 0.008
-1.311 + 0.405
0.024
-0.171
Table 5: DSoftKT uses the original interpolation scheme proposed SoftKI. A RMSE and NLL give the increase in RMSE and NLL respectively compared to the DSoftKI interpolation scheme.
Ac-ala3-nhme
Dha
At-at
Stachyose
0.0910
0.0783
0.0758
0.0627
0.0867
0.0761
0.0729
0.0606
0.0823
0.0
0.079
0.0701
0.0584
la PMMSS A
0.0059
0.0064
0.0068
0.0054
0.0057
0.0060
0.0040
0.0045
0.0051
0.0039
0.0043
0.0041
Buckyball- catcher
Double-walled- nanotube
At-at-cg-cg
0.0535
0.0482
0.0304
softki
svgp
0.0521
0.0444
0.0293
ddsvgp
dsoftki
0.0507
0.0405
..0283
0.0025
1000
0.0035
0.0017
(200
0.0036
0.0006
0.0040
0.0073
RMSE atom
Figure 4: Test RMSE per atom test gradient RMSE per component obtained by various methods
MD22 datasct. Bottom left is best.
dimensions and Absse is base learning rate. The learning rate that we use for hyperparameter optimization is thus Aet which we obtain for a constant Abase and minibatch size across methods. In this way, we keep all the training data, epochs, and minibatch size constant across methods while requiring hyperparameter optimization to take step in proportion to the amount of data that it encounters. We use D 1 for GP regression, D d for DSVGP, D = d+1 for DDSVGP, and D /3 (d5()))) + for DSoftKI. We set 83 d3? so that D = for KI.
Results. Table 2 gives the test root mean-squared error (RMSE) and Table g test negative log-likelihood (NLL) averaged across three runs. We use Abase 0.01 and minibatch size of 1024, scttings that are known work well for SVGP on this benchmark, to set Aed for all other methods. We were able to run DSVGP on the Welch dataset due to time constraints since d 20 and each minibatch of hyperparameter optimization has time complexity ((.)d))
Once we control for the effective learning rate, we see that fitting derivatives helps improve DSoftKI's performance relative to SoftKI, its non-derivative base, as measured by test RMSE performance. We believe that this increase performance is due to DSoftKT's modified interpolation scheme which jointly takes into
2afaf
11

---

Published in Transactions on Machine Learning Research (01/2026)
n
d
SoftKT
SVGP
DDSVGP
DSoftKI
Ac-ala--nhme
76598
126
0.894 + 0.006
1.051 + 0.018
7.907 + 0.004
82.665 + 0.06
Dha
62777
168
0.887 + 0.021
0.951 + 0.027
6.544 + 0.012
89.155 + 0.223
At-at
18000
180
0.336 + 0.006
0.427 + 0.001
2.626 + 0.01
27.236 + 0.037
Stachyose
24544
261
0.475 + 0.005
0.555 0.007
3.344 + 0.021
54.181 + 0.142
At-at-cg-og
9137
354
0.355 = 0.008
0.423 0.007
1.988 = 0.006
27.102 = 0.036
Buckyball-catcher
5491
444
0.244
0.007
0.295
0.004
1.249
+ 0.007
20.764
+ 0.111
Double-walled
4528
1110
0.409 + 0.007
0.434 0.004
1.906
+ 0.003
45.379 + 0.007
Table 6: Wall-clock training time in seconds per epoch (best. GPwD bolded) on MD22 dataset. We abbre- viate bonoll-wwaled-nanotte (Figure as Double-walled due to space.
account distances and directional information via learned temperature hyperparameters. In Table S we run DSoftKT with the SofTKT interpolation scheme, and that it performs worse on several datasets. For the variational methods, we observe that SVGP performs the best, followed by DDSVGP which fits two directions, and followed last by DSVGP which fits all derivative information once effective learning rate is controlled for. As a reminder, we use Aed to control for the amount of data each method sees during a step of hyperparameter optimization and not as a recommendation to obtain the best possible results. We conjecture that the difference in optimization dynamics between DSoftKI and DSVGP DSVGP/DDSVGL /DDSVGP is because the former jointly models values and gradients with each interpolation point whereas the latter approaches introduce separate inducing points for values and caddients/directioon
Another trend that we observe we control for Aet is that GPwDs tend to have lower NLLs compared to their non-derivative counterparts. For methods based on SVGP, this is somewhat surprising since the variants that fit derivatives also have higher test RMSE, indicating that the versions that fit derivatives result in simpler models. On the other hand, for DSoftKI, we observe that the NLL is positively correlated with its test RMSE relative to the performance of SoftKI. This suggests that the additional information is used to increase the complexity of the model when also improves the fit. of the model.
Table 4 reports the (wall-clock) training time (in seconds) per epoch of training on selected synthetic datasets. One epoch of training is one iteration through the dataset. We use a single RTX 6000 Blackwell Pro Max- Q edition GPU. We observe that the non-derivative methods such as ST and SVGP are faster than the derivative-fitting counterparts. This is expected since we fit d times more data in GPwD regression. For methods that predict full derivative observations, we see that DSoftKI has higher per-epoch training time compared to DSVGP, another method that can predict full derivative observations, on datasets with = 2. However, as scales, DSVGP's cubic scaling in d leads to prohibitive cost and it becomes infcasible to run DSVGP on the Welch (d 20) dataset. DSoftKI also has higher per-epoch training time than DDSVGP, particularly for large d. This is expected given DSoftKI's Omm'nd) complexity versus DDSVGP's O(m²p") with p 2 so that its per-epoch time does not scale with increasing d. However, DDSVGP only predicts = 2 directional derivatives.
5.2 High-Dimensional Molecular Force Fields
In this section, we evaluate the ability of DSoftKT to fit gradient information. Towards this end, we use the MD22 datasot (Chmicla et al., 2023), datasct of molocular energy surfaces where inputs are molccular configurations (d (d= 168 to 1110), outputs are energies, and gradients of the energy surface correspond to the physically-meaneee
quantity of (negative) forces. Consequently, fitting this information enable to model the physical dynamics.
Obtaining gradient predictions. To obtain the gradient prediction for SVGP, Soft and DDSVGP, we take the gradient of the posterior predictive mean to obtain the predicted gradient value since none of these methods prodict the gradient directly. For DSoftKI, we use the prediction of the gradients and not the gradient of the prediction. In an exact GPwD, these predictions would be equivalent. We are not able scale DSVGP to this setting.
12

---

Published in Transactions on Machine Learning Research (01/2026)
Results. Figure 4 reports the test value RMSE (per atom) vs. the gradient RMSE (per component) averaged across 3 runs for each method on each molecule. We refer the reader to Appendix B.5 for more experimental details. The bottom left corner is the best, since that is where there is minimum error on both values and gradients. We observe that DSoftKI performs well, particularly when the dimension of the dataset is large. One surprising observalion we make that lower value error is not correlated with lower gradient error. Upon examining the test. RMSE curves (sce Appendix B.5), we see that many GPs that do not fit full derivative information such as SVGP increase the test RMSE error for derivative prediction training progresses. Conversely, DSoftKI improves the test RMSE error both values and gradients training progresses for most datasets. Further investigation of the tradeoffs of fitting gradients is warranted for GPwD regression, since the additional computational complexity of fitting derivatives should be weighed against if obtaining accurate gradients predictions is required or not.
Table reports the (wall-clock) training time (in scconds) per epoch of training on the MD22 dataset. We include SoftKI and SVGP as non-derivative methods for completeness. As reminder, DVSGP does not scale to this setting and use DDSVGP with p = 2 inducing directions. As expected, DSoftKT has higher per-epoch training time than DDSVGP, particularly for large d. The additional time complexity is justified in cases such as molecular dynamics where full forces are required since the shape of the surface also needs to be fit. Notably, DSofKIT' training time scales roughly linearly with d (e.g., 20.8 seconds at d = 444 versus 45.4 seconds at = 1110), consistent with its ((m2nd) complexity. In contrast, DDSVGP's time remains nearly constant across dimensions since its (m""p") complexity is independent of d. However, this assumes mp² > d, which we observe begins to break when comparing its timing on the Double-walled dataset with the Buckyball-cat dataset.
6 Conclusion
In this paper, we introduce GP that can fit and predict full derivative information called DSoftK1. It has posterior inference time complexity of O(ndm²) and hyperparameter optimization time complexity of O(m² + bbnd,, and thus, scales to larger n and d than supported by previous GPwD methods while retaining the ability to fit and predict full derivative information. We have evaluated KT on GPwT) regression tasks and shown that it accurate, both in terms of test RMSE and test NLL, as well as modeling gradients. While promising, there certain limitations of our work that should be investigated further.
First, while DSoftKT is more scalable than existing GPwD methods that fit derivatives, further work can be done to improve its computational complexity. In particular, the method is currently bottlenecked by the computation of the DI MLL. Improved approximations of the MLL or a different objective could further scale the efficiency of fitting derivative observations. Second, a deeper investigation of the tradeoffi between fitting derivative information or not when they are available is warranted. For some applications, the added difficulty of fitting gradient observations may not warrant. improved surface modeling. Third, while we demonstrate that DSoftKI supports DKL out of the box (Appendix D), achieving computational speedups by projecting into a lower-dimensional feature space, further exploration of learned kernels for GPwD regression is promising direction.
References
Chris Camaño and Daniel Huang. High-dimennional gaussian process regression with soft kernel interpo- lation. Transactions on Machine Learning Research, 2025. ISSN 2835-8856. URL https::/openrerrie: nett/forun?hdd+Ub22111
Stefan Chmiela, Valentin Vassilev-Galindd, 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, :(2))ead)) 2023.
Filip De Roos, Alexandra Gessner, and Philipp Hennig. High-dimension gaussian process inference with derivatives. In International Conference Machine Learning, pp. PMLR, 2021.
13

---

Published in Transactions on Machine Learning Research (01/2026)
Dheeru Dua and Casey Graff. UCI machine learning repository, 2017. URL http:://rrchiieeeeeee uci. edu/
m1
David Eriksson, Kun Dong, Eric Lee, David Bindel, and Andrew G Wilson. Scaling gaussian process regression with derivatives. Advances in neural information processing systems, 31, 2018.
Jacob R. Gardner, Geoff Pleiss, Ruihan Wu, Kilian Q. Weinberger, and Andrew Gordon Wilson. Product 08903
kernel interpolation for scalable gaussian processes, 2018. CRL
Jacob R. Gardner, Geoff Plciss, David Bindel, Kilian Q. Weinberger, and Andrew Gordon Wilson. Gpytorch: Blackbox matrix-matrix gaussian process inference with gpu acceleration, 2021.
A. Girard. A fast 'monte-carlo cross- validatioi procedure for large least squares problems with noisy data. Numer. Math., 56(1):1-23) 23, January 1989. ISSN 0029-599X. doi: 10.1007/BFO1395777.
James Hensman, Nicolo Fusi, and Neil D Lawrence. Gaussian processes for big data. arXiv preprint arXiv1309.6600 2013.
M.F. Hutchinson. A stochastic estimator of the trace of the influence matrix for laplacian smooth- ing splines. Communication in Statistics- Simulation and Computation, 18:1059-1076, 01 1989. doi: 10. (0.1080/036100 100001000888000
Martin Jankowiak, Geoff Pleiss, and Jacob Gardner. Parametric Gaussian process regressors. Tn Hal Daumé III and Aarti Singh (eds.), Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pp. 4702-4712. PMLR, 13-18 Jul 2020. URL https:
//proceedings.co mlr. www.......s/..ee. press/vv19/ lankooiak20a..hte..oo. html
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 on Machine Learning, pp. 5279-5289. PMLR, 2021.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. arxiv preprint arXiv:1412.6966 2014.
Jihao Andreas Lin, Shreyas Padhy, Javier Antorán, Austin Tripp, Alexander Terenin, Csaba Szepesvári, Jos6 Miguel Hernández-Lobato, and David Janz. Stochastic gradient descent for gaussian processes done right. arXiv preprint arXiv:2810.2058 2023.
David JC MacKay al. Bayesian nonlinear modeling for the prediction competition. ASHRAE transactions, (00(2):1033-106) 1994.
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.
Misha Padidar, Xinran Zhu, Leo Huang, Jacob Gardner, and David Bindel. Scaling gaussian processes with derivative information using variational inference. Advances in Neural Information Processing Systems, 34:6442-6453 2021.
Henry C Plummer. On the problem of distribution in globular star clusters. Monthly Notices of the Royal Astronomical Society, 71(5):460-4M) 470, 1911.
Joaquin Hainonero-Candell and Carl Edward Rasmussen, A unifying view of sparse approximate gaussian process regression. The Journal of Machine Learning Rescarch, 6:1939-1959, 2005.
Carl Edward Rasmussen and Christopher K. 1. Williams. Gaussian Processes for Machine Learning (Adaptive Computation and Machine Learning). The MIT Press, 2005. ISBN 026218253X.
Edward Snelson and Zoubin Ghahramani. Sparse gaussian processes using pseudo-inputs. Advances in neural information processing systems, 18, 2005.
14

---

Published in Transactions on Machine Learning Research (01/2026)
Michalis Titsias. Variational learning of inducing variables in sparse gaussian processes. In Artificial intel- ligence and stalislics, pp. 567-574. PMTR, 2009.
Pauli Virtanen, Ralf Gommers, Travis E. Oliphant, Matt Haberland, Tyler Reddy, David Cournapeau, Ev- geni Burovski, Pearu Peterson, Warren Weckesser, Jonathan Bright, Stéfan van der Walt, Matthew Brett, Joshua. Wilson, K. Jarrod Millman, Nikolay Mayorov, Andrew R. J. Nelson, Fric Jones, Robert Kern, Eric Larson, C Carey, îlhan Polat, Yu Feng, Eric W. Moore, Jake VanderPlas, Denis Laxalde, Josef Perktold, Robert Cimrmann Ian Henriksen, E. A. Quintero, Charles R. Harris, Anne M. Archibald, António H. Ribeiro, Fabian Pedregosa, Paul van Mulbregt, and SciPy 1.0 Contributors. SciPy 1.0: Fundamental Algorithms for Scientific Computing in Python. Nature Methods, 17:261-272, 2020. doi: 10.1038/s41592-0
Ke Wang, Geoff Pleiss, Jacob Gardner, Stephen Tyree, Kilian Q Weinberger, and Andrew Gordon Wilson. Exact gaussian processes million data points. Advances in neural information processing systems, 32, 2019.
Andrew Wilson and Hannes Nickisch. Kernel interpolation for scalable structured gaussian processes (kiss- gp). In International conference on machine learning, pp. 1775-1784. PMLR, 2015.
Andrew Gordon Wilson, Zhiting Hu, Ruslan Salakhutdinov, and Eric P. Xing. Deep kernel learning. In Arthur Gretton and Christian C. Robert (eds.), Proceedings of the 19th International Conference on Artificial Intelligence and Statistics, volume Cadiz, Spain, 09 May 2016. PMLR. URL https: httppp::///ooceiing..oo. mlr, presss/55lwiilson666 htt.e/iiag........ereee/////aaaa html
Proceedings of Machine Learning Research, pp. 370-378,
Mohit Yadav, Daniel Sheldon, and Cameron Musco. Kernel interpolation with sparse grids, 2023. URL https://arxiv arg/abs/ ab///3//.../.... 2305. 14451
A
Method
We describe the KT method in more detail, including comparing its kernel structure with related GPwD methods (Section A.1 detailing the solving procedure (Section A.2), and its implementation (Section A.3)
A.1 Comparing Gradient Kernels
For a closer comparison of exact GPwDs, DDSVGP, DSVGP, and DSoftKI, we write out the kernel matrix between a test point * and the training data x or inducing points z in more detail. For an exact GPwD, the kernel directly uses the training data and the derivatives of the kernel used exactly below
K = [(t.) ((,,,,)) (I VHA
(29)
= ( 2x), K(R,,/VT
(.k((t V.K(
(30)
For DDSVGP (Le., with directional derivatives), the kernel uses the inducing points and additional vari- alional parameters to approxinate the training data. Additionally, il uses directional derivatives of the underlying kernel function as below
((+.) Ae k(+,,z)) (1 σf)
K DDSVGP
(31)
[
) 2,) (
k(H,2;)
(32)
15

---

Published in Transactions on Machine Learning Research (01/2026)
For a DSVGP, the kernel uses the inducing points, additional variational parameters, and ordinary derivatives of the underlying kernel as below
K
DSVGP
(
1 A(#, 24) (T Vz]
(33)
[( (v)1)) v.k(+,z))
k(M;2))VV 2,k(+,,z)V
(34)
For DT the kernel uses the inducing points, an interpolation function, and gradients of the interpolation function as below
DSoftKI -z
( o{({)) nj
Kza (z)
(1 v
(35)
o Kea (17
(v oz(x)) (v.ol(+)) (od(x)))
(36)
A.2 Solving
The procedure
a solve Ca = -zx
(37)
is implemented by solving
-1/2y
Ra D
(38)
for where
A 1-1/2K
QR
(39)
L
is the QR. decomposition and Kzz LL" is the Cholesky decomposition since Č (QR) This has time complexity (nm²² to compute.
A.3 Implementation
We implement DSoftKT in PyTorch (version 2.4.1) and GPyTorch (Gardner et al. 2018) (version 1.12) so that it can leverage GPU acceleration. DSoftKI currently only leverages single GPU. Arithmetic for KI is implemented in single-precision floats. There are a few places where the Cholesky decompasition of K is utilized in posterior inference. To stabilize these computations, we add a small jitter to the diagonal, when nceded. In extreme cases when Kzz values are extremely small, we may switch to doublo-precision floats to perform a Cholesky decomposition, before converting the results back into single precision.
Recall that we use CG implement the Hutchinson pseudoloss. We use PyTorrh implementation of CG descent. method following (Maddox et al. 2022). We CG tolerance of le - 5. We rank-10 pivoted Cholseky decomposition to obtain preconditioner for CG descent.
B Additional Experimental Results
We detail the hyperparamoters used (Section B.1
We also provide supplemental ablations (Section B.2 and Section B.3 and supplemental experimental results for the experiments described in the main text
(Section B.4 and Section B.5
B.1 Hyperparameters
For the experiments in this paper, we initialize hyperparameters as follows:
16

---

Published in Transactions on Machine Learning Research (01/2026)
Test RMSE
Test NLI.
d
Exact MLL
Hutchinson
Exact MLL
Hutchinson
Branin
2 0.002 + nan
0.029 + 0.0
-4.527
0.208 + 0.003
Six-hump-cames
2
0.014 + 0.002
0.198 = 0.018
-2.185 + 0.245
0.627 = 0.014
Styblinski-tang
2
0.013 + 0.004
0.035 + 0.001
-2.283 + 0.065
0.726 + 0.08
Hartmann
6
0.011 + 0.003
0.408 + 0.016
-3.145 + 0.318
1.159 + 0.007
Welch
20
0.618 + 0.006
1.373 + 0.002
Table 7: Comparison of using Exact MLL vs. Hutchinson's pscudoloss. The best test RMSE and test NLL are bolded. The notation +-nn indicates that all runs failed.
unstable whereas the nolation indicates that
Derivative Noise vs. Noise
Branin
1.1
Sixhump came!
Skyblirski-taang
Hartmann
0.9
0.6
()
Figure Effect of noise.
lengthscales fj 1,
output scale 1,
interpolation points z via k-means clustcring on training data,
lemperatures Tı₃ = for all entries, value noise 83 0.1, and
gradient noise 8 0.1 d where d is the dimensionality of the dataset.
As discussed in the main text, all hyperparameters
(e, 7, (mxd)) T(mxd),, ,0,
(40)
are learned jointly via pradient-basod optimization of the stabilized DSoftKI MLL (Equation 25 using
Adam (Kingma & Ba 2014).
B.2 Exact MLL vs. Hutchinson's Pseudoloss.
Table reproduces the results of Section 5.1 using the Fxact MLL compared to Hutchinson's pseudoloss. We see that when using DSoftKI's MLL, there instability in training on Branin and Welch. However, using Hutchionson's 's pseudoloss exclusively leads to worse results. As result, we opt for an objective that the DSoftKT MLL, when possible, and Hutchinson's pseudoloss otherwise described in Section 4.2
B.3 Effect of Noise Levels
The DSoftKI method balances the fitting value information with gradient information by controlling the relative ratio of 32 B3. Figure shows the results of varying this ratio on the synthetic benchmark gain more insight into how to set this ratio. We see that the relatative ratio 83 82 dependent on the dataset. The experiments reported in the paper had relative factor of 1.
17

---

Published in Transactions on Machine Learning Research (01/2020)
Surface
Interpolation Points
Veriational Inducing Pointss
1.0
4
sondd dsoftei
pigp
dsvgp
3-
ddsvgp
2-
0.6
1-
0-
0-
-1
-1-
-2
-2"
0.0
-3
-3
0.0
0.2
0.8
0.0
1.0
(a) Branin.
Surface
Interpolation Points
Variational Inducing Points
1.0
softil dsoftii
svgp
dsvgp
2
ddsvgp
1-
0.5
0-
-2
-2
0.0
0.0
0.2
0.6
0.8
1.0
-2
(b) Six-hump-can camel.
Surface
Interpolation Points
Variational Inducing Points
1.0
softkl dsoftk
dsygp ddsvgp
2
2
0.6
1
0.4
0-
0.2
-2
-2-
0.0
-3
-3
0.0
0.2
0.4
0.6
1.0
(c) Styblinski-tang
Figure 6: Learned interpolation and inducing points across various GP and GPwD regression modh
B.4 More Synthetic Benchmark Results
Data normalization. We scale data {(a, Vi, dy:) as ((hypercube(x:). (y: - p)/a. dy: where hypercube: : _R [0, 1jd maps into the unit hypercube, á is the empirical mean of the values, and & is the empirical standard deviation of the values and the derivatives. This ensures that the values and derivatives are measured in the same units, and follows standard practice.
Gradient fitting results. Table g reports the gradient fitting errors for the experiment described in Section 5.1 synthetic functions using scaled effective learning rate. We observe that DT has the
18

---

Published in Transactions on Machine Learning Research (01/2026)
d
SoftKT
SVGP
DSVGP
DDSVGP
DSoftKT
Branin
2
0.17 + 0.024
0.237 + 0.029
0.357 + 0.207
1.526 + 0.058
0.07 = 0.01
Six-hump-camel
2
1.069 + 0.156
1.206 + 0.057
1.545 + 0.071
7.032 + 0.187
0.345 + 0.084
Styblinskii-ang
2
1.008 + 0.11
0.998 + 0.02
0.86 + 0.133
0.895 + 0.165
0.319 + 0.057
Hartmann
6
0,429 + 0.232
0.377 + 0.007
0.506 + 0.01
0.527 = 0.01
0.028 + 0.008
Welch
20 0.011 + nan
0.052 + 0.0
0.269 = 0.002
0.001 + 0.001
Table 8: Gradient fitting errors synthetic functions using scaled effective learning rate. One of the runs for SoftKI encountered numerical instability (hence +nan). indicates a timcout for DSVGP.
DSVGP
DSoftKI
DSoftKI (Sg/(d3)) 10d)
Branin
0.323 + 0.821
281.505 + 72.881
0.397 + 0.061
Six-hump-camel
0.33 + 0.248
53.827 + 49.856
8.921 7.412
Styblinski-tang
0.665 + 0.121
41.907 + 14.305
4.929 2.126
Hartmann
0.047 = 0.014
10.755 + 0.916
0.591 0.265
Welch
-1.329 1.906
0.138 + 0.161
Table 9: Uncertainty quantification over gradient predictions.
lowest gradient errors of the methods. The relative size of B3 and 83 controls the degree to which value fitting is prioritized compared to gradient fitting.
We also compare the gradient test RMSE of the base GP and its extension with derivatives. For gradient test RMSE, we sum up the errors across the dimensions. Thus, we should normalize by the number of dimensions to compare across datasets. We observe that it is dataset dependent whether or not DSVGP outperforms SVGP. As reminder, DDSVGP does not fit all the derivative information, and thus, we derive its gradient prediction from its value prediction. We observe that derivative information helpful for DSoftKI, the extension of SoftKI to the setting with derivatives.
Uncertainty quantification of derivatives. We only report test NLLs for DSVGP and KT since these are the only methods that provide uncertainty estimates for full derivative information. In general, we find that DSoftKI has larger test NLLs compared to DSVGP. This is consistent with our earlier findings that DSVGP gives lower test NLL compared to SVGP, even though its test RMSE is worse
Visualization. Figure visualizes the learned interpolationiindues /induuing points for each method. We make a couple of observations on the learned points. First, of the interpolation points learned by SoftKI and DSoftKI compared to the inducing points learned by the variational approaches. In particular, the variational approaches have additional learned parameters that can help represent the structure in the dataset. In contrast, Soft and DSoftKI have to rely solely on the inter- polation points to represent the structure the dataset. Second, we observe that the interpolation pearoolation/inducinn /inducing points learned those GPs that fit derivatives and those that do not are different. In particular, DSoftKI and SIftKI capture different geometric structure in the underlying data. The inducing points learned by DSVGP are more spread out compared to those learned by DDSVGP and SVGP.
that there are differences between the structure
B.5 More Force Field Results
Dataset details. More concretely, we have a labeled data set ((x;, Mi, dy:) where v; encodes the Curlesian coordinates of the atomic muclei in the molecule Angström's, y/4 encodes the energy the molecule in kcal/mol, and dy; encodes the negative force exerted on each atom in the molecule. The dimensionality of 2; is 3A where A is the number of atoms in the molecule since each atom takes 3 coordinates to describe in 3Dsspace. For example, 42 atom system will have 126. We consider all available molecules.
19

---

Published in Transactions on Machine Learning Research (01/2026)
lest RMSS Ac-Ao---i144
HNSE on DHA
ROSAL on AI-AT
lest RMSS stachyose
0.38
0.50
0.20
0.54
0.45
0.38
0.40
softi
02c/g2
0.53
0.04
0.25
40
15
30
Opoch
Opoch
Eonch
(a) Ac-ala3-nhan
(b) Dha.
(c) At-at.
(d) Stachyose.
Text RMSE on buesyball-catcher
Tist RNSE or dbw.le--nilldd-aaaeee
RAP
0.35
0.50
0.5
0.45
2.0
0.40
0.4
5p3p
1.5
0.35
1.0
0.30
0.
(hofth
0.5
15
40
20
20
20
30
40
50
10
70
35
(e) At-at-cg-og.
Buckyball-catcherr
(g) cuuble-vallednnanttiee
Figure 7: Test RMSE curves.
ANSE
ANSE stachyose
11.4
11.66
16.2
13.0
13.4
12.9
18.0
10
soft:
12.0
seftel
setrd
15.6
s0tod
12.5
36
60
99
49
(a) Ac-ala3 1lln-nhme.
(b) Dha.
(c) At-at.
(d) Stachyose.
Test 9 RASE on AT-AT-CC-CC
Test V RMSE on buskyball-etthher
Test AMSE on double inluuennertueee
18.8
33.5
33.0
2000
ddsyp
ddripp
DAL
(e) At-at-cg-cg.
Buckyball -catcher.
(8) Doublleswallee-nanotubee
Figure 8: Test gradient. RMSE curves.
Data normalization. synthetic experiments. Instead of using hypercube, we scale the Cartesian coordinates as 2,/3 which has the effect of changing units.
We use the same data normalization scheme for energies and forces as we did in the
Training curves. Figure 7 and Figure g give the test RMSE, test gradient RMSE, and their standard deviations at each epoch of training the MD22 dataset. As reminder, did not tune the learning rates batch sizes to obtain the best possible performance. Rather, we control the effective learning rate so that we can compare GP and GPwD regression. SVGP is in blue, SoftKI is in green, DDSVGP is in orange, and DSoftKI is in Red. On Ac-Ala3-NHMe, we observe that methods that fit derivative information overfit both values and gradients. On several datasets such as АT-AT and Stachyose, we observe that DSoftKI overfits on test RMSE but does not on test gradient RMSE. On other dataset such as AT-AT-CG- Seekyybal-cattchen and waalle-wlled--nnntte and gradients. In contrast, DDSVGP tends to overfit, as measured by test gradient RMSF. This is somewhat
we observe that DSoftKI is able to fit both the values
20

---

Published in Transactions on Machine Learning Research (01/2026)
12500
25000
30000
50000
softki
20030
softki
softiii
40000
softti
7500
Count
15030
Count
30000
5000
10030
20000
10000
2500
10000
0
10
(a) Ac-ala3-nhie
(b) Dha.
(c) At-at.
(d) Stachyose.
60000
diofttii
softii
50000
dsaftki softici
125000
denttii
40000
100000
Count
40000
30000
Count
75000
20000
50000
10000
25000
0.5
1.0
0.5
1.0
3.5
(e) At-at-cg-cg.
Buckyball -catcher
(g)
Doublee-allle--annttue
Figure 9: Histogram of learned lemperatures for DSoftKI and SoftKI on the MD22 dataset. The learned
temperatures have similar distributions.
surprising, since it does fit directional derivative information. We also observe that SoftKI and SVGP cannot fit gradients even though they do fit the values. This suggests that GPs that do not fit derivatives do not capture the shape of the surface well.
Temperatures and lengthscales. Figure and Figure 10 give the histograms of lemperatures and length- scales that are learned by DSoftKI and SoftKI. As a reminder, we modified SoftKI also have different temperatures per interpolation point as did for DSoftKI to have more fair comparison.
The lemperature Tę E Rd for interpolation point = controls how strongly that point responds to variations along each input dimension during interpolation. A smaller temperature component Tej makes the interpo- lation weight of(x)) more sensitive to changes in dimension j, effectively sharpening the directional influence. Conversely, larger temperatures smooth out the response along that dimension.
Figure ssoows that the learned temperature distributions are similar between DT and SoftKI, indicating that the interpolation scheme (which is what is affected by the temperature) is primarily adapted to the data geometry. Moreover, observe that DSoftKI learns larger kernel lengthscales compared to SoftKI (Figure TO). This suggests that when gradient information is available, the model can rely on the temperature vectors to capture local directional variation, allowing the kernel lengthscales to increase and model smoother global structure. In other words, the temperatures and lssathscales play complementary roles: temperatures handle local directional sensitivity while lengthscales control global smoothness.
C Additional Datasets
We describe experiments additional datasets inchuding toy n-body simulation (Section C.1 and the UCI dataset with synthetic gradients (Section C.2 to test DSoftKI a broader range of datasets.
C.1 Toy N-Body Dataset
Computational physics and chemistry provide settings where high-dimensional gradient observations arise naturally. Consider an n-body system with particles (e.g., atoms or point masses). Let q = (gI...., gn) denote the positions of each particle where each di R3, and let p (P1,..., PP) denote the momenta where
21

---

Published in Transactions on Machine Learning Research (01/2026)
125
dhoftki
asoftii
20
dsoftki
dsoftki
20
Count
Count
40
10
50
25
20
5-
0
2
Lengthscale
Lengthscale
Lengchscale
Lecgee
(a) Ac-ala3-nhme.
(b) Dha.
(c) At-at.
(d) Stachyose.
50
asoftki softki
csaftki softki
200
csoftki
softki
40
150
Count
Count
30
Count
50
100
20
20-
50-
10
0
Lengchsccle
Lengchscate
Lengthscalle
(e)
At-at-cg-cg.
Buckyball-catcher
Double--allled--aatttee
Figure 10: Histogram of lengthscales MD22 dataset for DS and SoftKI. The lengthscales learned by the DSoftKI tend to be larger than the lengthscales learned by SoftKI.
d
DDSVGP
DSoftKI
V-DDSVGP
V-DSoftKI
nbody-4
0.058 = 0.005
0.059
0.002
1.453 + 0.083
0.708
0.034
n
36
0.027 + 0.001
0.035 + 0.002
1.5 = 0.059
1.032
0.045
n
48
0.03 = 0.0
0.028 0.002
1.368 + 0.117
1.016
0.092
nbody-10
60
0.026 0.003
0.024 0.0
1.343 + 0.037
1.104 + 0.04
Table 10: Test RMSE (best bolded) on n-body datasets.
V-DDSVGP and V-DSoftKI denote gradient
RMSE (best bolded).
each Pt E R3. The dynamics are governed by Hamilton's 's equations
dq 0H dp 0H
(41)
dt 3p
dt ôq
where (q. p) is the Homiltonian. Thus, given observations of the Hamiltonian H and its gradients VHI, we can learn the system's dynamics, i.e., dq 2 and de
dp More concretely, it is GPwD regression problem where
x = (q, y = H (x), and
Vy= OH 0H
(42)
0q 0p
The input dimensionality is 6n since each particle contributes 3 position and 3 momentum coordinates. A learned surrogate for the Hamiltonian enables efficient prediction of the system's dynamics without expensive numerical integration, which can be intractable for large
We construct an n-body gravitational dataset nbody-n as follows
The Hamiltonian is
|
Gm,m;
H(q.
(43)
2m;
1 2
1=1
37hhi is classical system, MD22 dataset. (Section which is quantum system.
22

---

Published in Transactions on Machine Learning Research (01/2026)
d
DDSVGP
DSoftKT
n
24
0.98 + 0.016
1.68 + 0.019
nbody-6
36
0.995 + 0.029
2.802 + 0.048
nbody-8
48
0.97 + 0.006
3.646 + 0.059
nbody-10
60
0.976 = 0.001
4.591 + 0.162
Table 11: Training time per epoch (in seconds, best bolded) on n-body Hamiltonian datasets.
n
d
DDSVGP
DSoftKI
V-DDSVGP
V-DSoftK1
Kin40k
18000
8
0.867 + 0.002
0.864 + 0.006
0.399 + 0.001
0.393 + 0.001
Protein
20578
9
0.81 + 0.006
0.787 + 0.003
2.778 + 0.087
2.655 + 0.283
Bike
7820
17
0.629 + 0.026
0.421 + 0.006
0.763 = 0.004
0.624 = 0.012
Elevators
7470
18
0.821 = 0.013
0.797 + 0.009
0.449 = 0.002
0.448 = 0.001
Pol
6750
26
0.898
0.007
0.702 + 0.01
0.383
0.013
0.324 - 0.012
Slice
19260
385
0.357
0.005
0.03 = 0.008
0.165
0.007
0.094
0.004
Table 12: Test RMSE (best bolded) on selected UCI datasets. V-DDSVGP and V-DSoftKI denote gradient RMSE (best bolded).
where my is the mass of particle G is the gravitational constant, and is a softening parameter that prevents singularities when particles approach each other. The first term is the kinetic energy and the second is the gravitational potential energy with Plummer softening (Plummer 1911). In this toy example, the Hamiltonian has analytic gradients which may not exist for more complex systems.
We generate datasets for {4, 6, 8, 10) particles, corresponding to input dimensions d E {24, 36, 48, 60). To generate physically realistic configurations, we simulate 100 trajectories by integrating Hamilton's 's equations using a high-order Runge-Kutta method (DOP853) implemented in SciPy (Virtanen et al.) 2020). Initial positions are sampled from (0,02²), initial momenta from (O, 0.5²) scaled by particle mass, and masses uniformly from [0.5, 2.0]. We apply center-of-mass corrections for numerical stability. From each trajectory, we sample configurations to obtain 10000 total samples per setting, computing the Hamiltonian and its gradients analytically. We filter out the 5% of samples with the largest gradient norms for numerical stability, leaving 9500 samples. For this toy dataset, we use G and c c= 0.1.
Table 10 reports the results of using a 90/10 training-lesting split on each dataset. DSoftKI and DDSVGP achieve similar value RMSE, but DSoftKI substantially outperforms DDSVGP gradient prediction (e.g.. 0.708 RMSE versus 1.453 V-RMSE for n 4). This demonstrates DSoftKT's advantage in applications requiring accurate full gradient predictions, such as learning physical dynamics from energy observations, and justifies the higher computational cost (Table II
C.2 UCI Dataset
To evaluate scaling on diverse input distributions boyond physics datascts, we construct GPwD regression benchmarks from standard UCI regression datasets (Dua & Grafl 2017) by generating synthetic gradients.
We select six UCI datasets spanning a range of dimensions: Kin40k (d = 8), Protein (d = 9), Bike (d= (7), Elevators (d 18), Pol (d = 26), and Slice (d = 385). For each dataset, we estimate gradients using a k-nearest neighbor finite difference approximation: for each point with label yr, we compute directional derivatives along the directions to its k =3 3 nearest. neighbors and average the resulting gradient estimates, i.e.,
1 Σ Vy
yj
y1 Fz
(xz Xz)
(44)
JEN.(0
where (i) denotes the nearest neighbors of point This introduces noise into the gradient observations, providing a test of robust to gradient estimation error.
23

---

Published in Transactions on Machine Learning Research (01/2026)
"
d
DDSVGP
DSoftKT
Kin40k
18000
8
2.257 + 0.022
1.715 + 0.1
Protein
20578
9
2.593 + 0.05
2.054 + 0.067
Bike
7820
17
0.977 + 0.003
1.317 + 0.067
Elevators
7470
18
0.994 + 0.024
1.297 + 0.017
Pol 6750
26
0.876 + 0.014
1.663 0.077
Slice
19260
385
2.378 + 0.047
71.279 + 0.076
Table 13: Wall-clock training time (in seconds) per epoch (best bolded) on selected UCI datasets.
DSoftKI
DSoftKI-DKL
Ac-ala3-nin
5.98e-03 + 2.90e-05
6.93e-03 + 9.40e-05
Dha
5.85c-03 + 8.80c-05
6.19c-03 = 7.40c-05
At-at
4.01e-03 1.04e-04
5.63e-03 + 5.40e-05
Stachyose
3.99e-03 + 2.90e-05
5.29e-03 = 1.30e-05
At-at-cg-og
2.58e-03 + 5.60e-05
3.98e-03 + 2.14e-04
Buckyball-catcherr
1.81e-03 + 5.40e-05
4.09e-03 = 1.17e-04
9.33e-04 + 2.70e-05
2.90e-03 + 4.59e-04
Table 14: Normalized test RMSE per atom (best bolded) on MD22 dataset.
Table 12 reports the results. DS outperforms DDSVGP on both value and gradient prediction across datasets. The Slice dataset (d = 385) is particularly notable: DSoftKI achieves 0.030 RMSE versus DDSVGP's 0.357, ratong strong performance on high-dimensional structured data. The higher train- ing time for DSoftKI on Slice (Table 13) reflects its O(m2nd) complexity, but is justified by substantially better accuracy.
Deep Kernel Learning
One advantage of DSoftKI over DSVGP and DDSVGP is that it does not require computing first or second- order derivatives of the kernel (Section 4.1). This enables the use of learned kernels, such as Deep Kernel Learning (Wilson et al. 2016) (DKL), without hard-coding kernel derivatives so that standard automatic differentiation techniques can be applied.
To demonstrate this capability, we evaluate DSoftKI with a deep kernel learning (DSoftKI-KLL) setup on the MD22 molecular datasets. We use 2-layer MLP feature extractor with hidden dimension 64, output dimension 24, and tanh activations. The RBF kernel is applied in the learned feature space. We train all hyperparamnters jointly using Adam, with a learning rate of 0.002 for the GP hyperparameters and 0.0005 for the neural network weights.
Table 14 Table 15 and Table T6 report the results. DSoftKI-DKL, performs comparably DS, with slightly higher RMSE on both values and gradients. We emphasize that this experiment intended as a demonstraaion that deep kernel learning works out of the box with DT, rather than a claim that it im- proves performance these particular datasets. Notably, DSoftKI-DKL offers computational speedups over DSoftKI since the interpolation points are projected into a lower-dimensional feature space. Furthermore. the neural network gradients can be computed efficiently via Jacobian-vector products, whereas DSVGP and DDSVGP would require hard-coding the feature extractor Jacobians or face intractable computation.
E
Additional Background on DDSVGP and DSKI
In this section, we describe in more detail how variational inducing point methods and SKI have been extended to the setting with derivatives.
24

---

Published in Transactions on Machine Learning Research (01/2026)
DSoftKT
DSoftKI-DKI.
Ac-ala3-nhie
8.27c-02 + 1.73c-04
8.62c-02 + 1.00c-05
Dha
7.41e-02 + 2.59e-04
7.62e-02 +: 2.27e-04
At-at
02 + 8.50e-05
7.49e-02 + 1.25e-04
Stachyose
5.86e-02 + 6.00e-05
6.19e-02 + 1.04e-04
At-at-cg-og
5.08c-02 = 2.02e-04
5.32c-02 = 1.10c-04
Buckyball-catcher
4.09e-02
1.78e-04
4.75e-02
+ 1.66e-04
nouble-wallle-nanotuue
2.84e-02
5.10e-05
3.28e-02 = 1.04e-03
Table 15: Normalized gradient test RMSE per dimension d (best bolded) on MD22 dataset.
d
DSoftKT
DSoftKT-DKT.
Ac-ala3-nhme
76598
126
82.665
0.06
58.98 + 0.516
Dha
62777
168
89.155
0.223
48.736 + 0.56
At-at
18000
180
27.236
0.037
13.958 + 0.025
Stachyose
24544
261
54.181 = 0.142
19.324 = 0.369
At-at-cg-og
9137
354
27.102 = 0.036
12.574 = 0.339
Buckyball-catcher
5491
444
20.764 + 0.111
7.655 + 0.146
Double-walled-nanotuie
4528
1110
45.379 = 0.007
7.143 + 0.088
Table 16: Comparing (wall-clock) training time (in seconds) per epoch (best bolded) on MD22 dataset for DSoftKI versus DSoftKI-DKL.
Variational inducing points. An inducing point method (Snelson & Ghahramani 2005 Quinonero-
Candela & Rasmussen 20055) introduces set of m < inducing points (24 R and associated inducing variables (z)
(w: R a proxy for the given inputs x and outpuls respectively.
In an inducing point method such as Sparse Variational Gaussian Processes (SGPR) (Titsias) 2009) and Stochastic Variational Gaussian Process (SVGPs)
(Hensman et al. 2013), the inducing points and variables
are related the dataset as
(1(x) ~ ~N -N(0, Kez KZx
(inducing)
Kxx
Kxx
y v1 ((()))))()))
(likelihood)
Note that the marginalization of the inducing variables reduces the model to the standard GP model. To make posterior inference tractable, SVGP makes variational approximation ((((x,, u) p(f(x) [))((u) where q(u) (uu m, S), treating m (mx1) and (mm)m) as additional learnable variational parameters. More concretely, = (e, 7, B, z, m, S). They can be learned by maximizing a lower bound on the ML. called the evidence lower bound (ELBO), defined as
ELBO( Σ Eat KL(g(u) (lp( z))
(45)
which cquivalently minimizes the KL-divergence (L(q(u) [| |z)). The ELBO can be optimized in mini- batches, resulting a time complexity of 0(m²). Given learned hyperparameters, the optimal posterior distribution
df(+)) K.zK m, K.. K.,KWl K Kxx)KKK, + A)
(46)
The complexity of inference is O(m²).
A SVGP can be extended to the setting with derivatives, i.e., a DSVGP (Padidar et al. 2021 More concretely, define inducing variable (u: Vuf)² that models the value and its gradient at.
25

---

Published in Transactions on Machine Learning Research (01/2026)
inducing point 2. A DSVGP introduces learnable parameters m(xx), dm (dml1), S(mxm), dS dS(1xm), and (dg(dmxdm) to define variational posterior ((ū) = -N m, S) where
m; and $
[
Siz ds.,
dSi
(47)
dm.
d2
:1
The parameters can be learned by maximizing the DSVGP ELBO which obtained by replacing the SVGP ELBO with the respective (-) versions.
Again, it is lower-bound the MLL. The time complexity of computing the ELBO Omm2d²). The posterior distribution is that of SVGP with all of the matrices replaced with their respective () versions. The time complexity of posterior inference is O(m²d since Keu is a (m(d + 1)) x (m(d + 1)) matrix.
minibatch is
DDSVGP (Padidar et al. 2021) utilizes diroctional derivatives oys = VTVf(x) where V projects onto p
dimensional subspace. This results in the modified matrices
K = (ov.
k(za: 26) (I %%,)
(48)
@v..
ab
K,z
k(*, kx#,20) (I 2%r)
(49)
16
which has additional learnable variational parameters ( (d(xp)) j for projecting. The m(p + 1) x m(p + 1) matrix Kzz = [k(2a, ab)las can be formulated efficiently with Hessian-vector products in time complexity dpm) The modified posterior distribution and ELBO are obtained by changing the respective (-) vari- ables to the (+) versions. Posterior inference be computed in time O(m²p"), assuming mp² > d to account for the cost of forming Kzz, since the d-dimensional gradients have been reduced to p-dimensional directional derivatives.
Structured kernel interpolation. One problem with inducing point methods the requirement that m < n. SKT (Wilson & Nickisch 2015) overcomes this limitation by using the approximate kernel
K W Kx Wex R Kxx
(50)
where z are cleverly chosen interpolation points on a pre-defined lattice, Wxx [w,((;))s matrix of interpolation weights with interpolation function (x,), and Wxz = WT xx. The SKI posterior is
lx, (N([()] ( + +A) K., K, + - KKx)
(51)
which is the GP posterior with Kxx replaced with xx
The posterior is tractable to compute since the lattice structure given by the interpolation points enables fast matrix-vector multiplies (MVMs) so that conjugate gradient (CG) methods can be used to solve large systems of linear equations. The complexity of a single MVM (m4 + mo m) with a cubic interpolation function (Wilson & Nickisch 2015) where is the number of interpolation points. However, the asymptotic dependence on d limits the application of SKI to low dimensions.
DSKI EEriksson et al. 2018) extends the SKI Wilson & Nickisch 2015 method to create a scalable GP regression method that handles derivative information. To accomplist iis this, DSKT uses the following approx- imate kernel
K xx
DSKI
WxaK Wax Kxx
(52)
where
) wet (2:)
(odee)) (wud(r))
(
W XM
(53)
ij
That is, DSKT approximates the kernel by differentiating the interpolation matrix. The posterior distribution is the SKI posterior with the SKI kernel replaced with the DSKI kernel. Like SKI, DSKI leverages fast MVMs enabled by the lattice structure and CG methods perform GP inference. The complexity of a single MVM is Onddd + m log m). The scaling in the dimension is even worse than SKI, limiting the application of DSKI to even smaller dimensions.
26


Reviews:
[{"review_id": "9xMyJSWG5z", "reviewer": "Reviewer_Dm1z", "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": "", "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": "Key contribution: the paper presents a new scalable GP method for fitting functions when noisy measurements of both function evaluations and the corresponding gradients are available. The method is a generalization of SoftKI (Camaño & Huang, 2025). As its predecessor, the method does not require computing kernel first and second order gradients, which improves computational complexity; unlike SoftKI, the interpolation-point specific temperature vectors allow more flexibility in modeling the gradients, which requires changes in the implementational considerations that keep the computations scalable. Small dimensional  and high dimensional data is used to show the benefits brought about by the introduction of the temperature hyperparameters.", "claims_and_evidence_yn": "Yes", "claims_explanation": "The methodology and its relations to previous methods is rigorously presented in the main text\n\nThe text is also thorough in justifying the exact setup for the numerical comparisons, the complexities entailed in the comparison and the limitations of the approach.", "audience_yn": "Yes", "audience_explanation": "The results are of broad relevance for GP applications with high dimensional data and gradient information (e.g. in the natural sciences); of particular utility might be scenarios when kernels need to be flexible (Wilson et al., 2016).", "requested_changes": "I found the text generally very readable and the descriptions precise and rigurous. \n\nA few suggestions of improvement:\n\n1) in section 3.3 softmax: i found the /T notation is weird for a d-dimensional vector T\n2) would be useful to remind people of the dimensionality of variables in section 4  \n\nPractical notes: \n- any specific considerations about the initial values of the hyperparameters?\n- in empirical data i expect the gradient measurements to be generally noisier  that the primary observations, how does that affect the estimation\n- since the text emphasizes that the method is one of the few that can be used with learned kernels it would be useful to include a concrete demonstration of that \n- was not sure where the heteroskedastic noise was for, please explain", "broader_impact_concerns": "No concerns", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "cffNLnvcBv", "reviewer": "Reviewer_csfB", "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": "", "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 introduces a GP capable of fitting and predicting full derivative information, called DSoftKI, by extending SoftKI. This very recent method approximates a kernel via softmax interpolation from learned interpolation point locations. The method achieves a posterior inference time complexity of $O(ndm^2)$, allowing it to scale to larger values of $n$ and $d$ than previous GPwD approaches while still retaining the ability to model full derivative information. Here $m$ is the number of inducing points (Quinonero-Candela & Rasmussen, 2005; Snelson & Ghahramani, 2005). The authors evaluate DSoftKI on GPwD regression tasks and demonstrate that it achieves strong accuracy, as measured by both test RMSE and test NLL, while effectively modeling gradients.\n\nIn general, the idea of the paper is not novel. The time complexity of $O(ndm^2)$ directly comes from the original approach, SoftKI.  However, the approach that DSoftKI handles derivative observations is non-trivial, and it entails several substantive methodological contributions.", "claims_and_evidence_yn": "Yes", "claims_explanation": "- To implement the softmax interpolation scheme, they associate each interpolation point with a corresponding learnable\ntemperature vector as in Equation 20. \n- DSoftKI introduces a different interpolation scheme and more learnable parameter at Algorithm 1. \n- Section 4.3 introduces the choice of value and gradient noises.", "audience_yn": "Yes", "audience_explanation": "The topic of Scaling Gaussian Process Regression using derivative observations is important. Although the proposed DSoftKI method builds upon SoftKI, extending it to handle derivative observations is non-trivial and entails several methodological contributions.", "requested_changes": "- How to estimate the variable $T$?\n- It is unclear why DSoftKI is more efficient than DDSVGP? As the author explained in the introduction, introducing $p <<d$ from DDSVGP comes at the cost of directly predicting derivatives and introducing further approximations. However, their method also adds many approximations.", "broader_impact_concerns": "No", "additional_comments": "", "extra_scores": {}, "extra_text": {}}, {"review_id": "Ue6ebLDcZM", "reviewer": "Reviewer_sBME", "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": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "No", "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 focuses on Gaussian Process Regression in a particular setting - when the derivatives of the target function are observed (GPwD), and proposes a novel technique called DSoftKI to achieve it more efficiently than existing techniques that is particularly relevant for high-dimensional data. This technique is summarized in Algo 1, and argued to be more efficient than SoftKI, an existing technique to achieve the same. The heart of Gaussian Process lies in the Kernel, and hence the core approach here is \"Soft Kernel Interpolation with Derivatives\" with the help of interpolation points, that also introduces \"temperature\" variables. In a regression setting, the variables are estimated in a Gradient-descent based algorithm, where the gradients of the objective function (p(y|x,M) is estimated using Hutchinson pseudoloss. In the proposed approach, this is calculated efficiently, using a Cholesky Decomposition and a low-rank reconstruction of the Kernel matrix over the training points. The posterior inference formula is calculated accordingly. Experimental analysis compare the RMSE and Negative Log-Likelihood achieved by the proposed method vis-a-vis benchmarks on benchmark synthetic datasets, and also on high-dimensional molecular force fields.\n\nIn general, the paper and the contributions are sound, though the main weaknesses in my opinion are as follows:\n\nThe technical contribution may be on the lighter side, as it seems to be a (non-trivial) increment over SoftKI\nThe experimental validation is mainly on simulated datasets, except one molecular force field data", "claims_and_evidence_yn": "No", "claims_explanation": "My main concern is that most of the experiments are on synthetic datasets except one molecular force field data.  Also, there are insufficient experimental results related to computational efficiency, which is the main motivation of the paper", "audience_yn": "Yes", "audience_explanation": "Gaussian Process regression is an important approach in many domain, and there are many open computational questions related to it that interests the ML community", "requested_changes": "My main suggestion would be to consider a couple of more datasets from different scientific domains where Gradient information of the target function is known, and show the results. Furthermore, the main claim of the proposed algorithm is computational efficiency, so there should be more detailed comparisons of running time for the proposed approach against the benchmarks.\n\nOne feature of the proposed approach is the additional temperature variables. The appendix does show some histograms of these variables, but we would like to understand their physical significance - which interpolation points are assigned higher/lower temperatures, and how these affect the final result.", "broader_impact_concerns": "NA", "additional_comments": "", "extra_scores": {}, "extra_text": {}}]

Author rebuttal:

Decision: Accept (poster)