Derivative Gaussian Processes on a Two-Direction Budget

Hyunseok Seung; Matthias Katzfuss

stat.MLOct 7, 2026

Derivative Gaussian Processes on a Two-Direction Budget

Authors: Hyunseok Seung, Matthias Katzfuss

Organizations: University of Wisconsin–Madison

Abstract

Gradient observations promise more accurate Gaussian process (GP) surrogates, but the cost of incorporating them has long stood in the way of realizing that promise. We propose a derivative GP with a budget of just two directions per observed gradient. One direction focuses on each gradient's direct contribution to target prediction, while the other aggregates its indirect contributions through correlations with the conditioning function values. Within a Vecchia approximation, where each prediction conditions on $m$ nearby inputs in $d$ dimensions, this construction represents their $md$ gradient coordinates using at most $2m$ directional derivatives, giving $\mathcal{O}(m^3)$ dense factorization cost per prediction target. For general conditioning sets, we bound the posterior approximation error relative to using full gradients and characterize when the error is small or the approximation is exact. In simulations, our method matches the accuracy of a leading exact gradient-reduction method at equal conditioning set size. Because its cost grows much more slowly with that size, it can use conditioning sets well beyond the memory limit of the exact method, reaching lower prediction error with a small fraction of the time and memory. Notably, our method can exploit gradient observations while requiring less computation time or memory than function-only GP baselines.

Figures & tables

Figure 1: Posterior fidelity under different directional-derivative choices as $m$ increases for a Matérn- $5/2$ GP in $d=50$ . We compare LITE with methods using $\mathbf{u}_{ij}$ or $\mathbf{u}_{ij}$ paired with a random direction in $\operatorname{col}(\mathbf{D}_{i})$ . KL divergence between the TERA and approximate posteriors is averaged over prediction inputs, with all methods using the same conditioning inputs.

Figure 2: Relative compound symmetry (CS) error (solid) and expected KL divergence from the full-gradient posterior to LITE ’s (dashed) as $d$ increases. We use an isotropic Matérn- $5/2$ kernel, $n=1024$ uniform inputs in $[0,1]^{d}$ , $200$ targets, maximin ordering, and $m=30$ nearest preceding inputs.

	Factor size	Time	Memory
Full grad	$m(d+1)$	$\mathcal{O}(m^{3}d^{3})$	$\mathcal{O}(m^{2}d^{2})$
TERA	$m^{2}+m$	$\mathcal{O}(dm^{2}+m^{6})$	$\mathcal{O}(dm^{2}+m^{4})$
LITE	$3m$	$\mathcal{O}(dm^{2}+m^{3})$	$\mathcal{O}(dm+m^{2})$

Table 1: Per-target cost of one Vecchia conditional factor with

m

conditioning inputs in

d

dimensions. TERA’s entries assume

m\leq d

; when

m>d

, the full-gradient factor is smaller.

Panel of Figure 3 — Figure 3: RMSE under the assumed GP against total prediction time for $256$ targets and peak GPU memory (top). Wall-clock time and peak GPU memory of one training step against $m$ on the Buckyball-catcher dataset (bottom). Point labels in the top panels indicate $m$ . Prediction uses $B=32$ , while the bottom panels use $B=32$ and $B=1{,}024$ . TERA exceeds the $80$ GB GPU memory for $m\geq 160$ in prediction; missing TERA points in the bottom panels exceed the GPU memory.

Figure 4: Test RMSE per atom, wall-clock time, and peak GPU memory on six MD22 molecular datasets. LITE achieves the lowest runtime and memory usage on every benchmark, running $2.8$ – $3.5\times$ faster than TERA while using only $13.3$ – $30.2\%$ of its memory. Its RMSE is higher than TERA’s but lower than that of the other methods. Standard GP and DSoftKI run out of memory on an 80 GB GPU when applied to DHA, whereas LITE requires only $1.06$ GB. Results are averaged over five independent runs, with error bars indicating standard errors.

Figure 5: Simple regret on Ackley ( $500D$ , left), Levy ( $800D$ , center), and best observed value on LassoDNA ( $180D$ , right) versus iteration. With only two directions per gradient, LITE substantially outperforms the function-only baselines, VBO and TuRBO-1, on all three tasks. It achieves the lowest final regret on Levy and a final value comparable to TERA on LassoDNA, while TERA attains lower regret on Ackley. Curves average five independent runs, with shaded bands indicating standard errors.

Appendix figures & tables7 assets

Supplementary material from the paper’s appendix.

Appendix

GP prediction
Setting	Value
Training/test inputs	$100,000$ / $256$
Input dimension	$d=500$
Input design	Sobol in $[0,1]^{d}$
Kernel	ARD Matérn- $5/2$
Signal variance	$1$

Appendix

Table 2: Settings for conditioning set and batch size scaling.

Setting	Value
Train/test split	90/10
Seeds	1,27,42,86,99
Kernel	isotropic SE
Initial lengthscale	1.0
Initial signal variance	1.0
Initial value-noise variance	$10^{-3}$

Appendix

Table 3: Settings for MD22 regression.

Benchmarks	Ackley ( $d=500$ ), Levy ( $d=800$ ), LassoDNA ( $d=180$ )
Domain	$[0,1]^{d}$
Iterations / initial design	$100$ iterations / $30$ Sobol points
Seeds	$1,\ 27,\ 42,\ 86,\ 99$
Noise parameters	$10^{-3}$ , fixed on synthetic tasks; initialized at $10^{-3}$ and learned on LassoDNA

Appendix

Table 4: Settings for BO synthetic objectives and LassoDNA benchmark.

Figure 6: Test RMSE versus conditioning set size on six MD22 datasets, with all methods using the covariance parameters learned by TERA at $m=20$ after one epoch. Parameters remain fixed throughout prediction. Increasing $m$ provides little improvement for LITE and the function-only Vecchia GP, and LITE ’s RMSE remains higher than TERA’s at matched $m$ .

Figure 7: Test RMSE versus conditioning set size on six MD22 datasets, using the same TERA estimates as in Figure 6 . Solid curves keep these parameters fixed; the dashed curve trains LITE for one additional epoch at each $m$ , starting independently from the same estimates. Further training improves LITE on four datasets, changes little on AT-AT-CG-CG, and worsens prediction on Buckyball-catcher.

Figure 8: Test RMSE versus training step on six MD22 datasets. Both methods use $m=20$ , training and prediction batch sizes of $32$ , and five training epochs. Additional updates do not uniformly improve predictive accuracy, with RMSE increasing after early minima on some datasets. Curves average five independent runs.

Figure 9: Wall-clock time and peak GPU memory on Ackley ( $500D$ , left), Levy ( $800D$ , center), and LassoDNA ( $180D$ , right). Time includes objective evaluations and gradient computation where used, including hypergradient computation on LassoDNA. Bars average five independent runs, with error bars indicating standard errors.

Explore similar work

Jun 1, 2026stat.ML

Scalable Derivative Gaussian Processes via Exact Gradient Reduction

Gradient observations can substantially improve Gaussian process (GP) surrogates, particularly in high-dimensional settings where function evaluations are expensive. However, exact inference with

n

function values and

n

full gradients in

d

dimensions scales cubically in the joint state size, imposing an intractable

\mathcal{O}(n^3 d^3)

computational bottleneck. We introduce TERA, a highly scalable derivative GP method based on target-specific exact gradient reduction. We prove that for stationary kernels, the gradient components orthogonal to the directions connecting the target and conditioning points are conditionally independent of the target function value; consequently, the exact conditional density is fully characterized by at most

m^2

directional derivatives once a conditioning set of size

m

is specified. By using these reduced, dimension-free conditionals as local factors in a Vecchia approximation, TERA effectively decouples

n

and

d

from the dense matrix inversion. This reduces the per-target evaluation cost to

\mathcal{O}(dm^2 + m^6)

time and

\mathcal{O}(dm^2 + m^4)

memory, leaving the underlying derivative GP model mathematically unchanged. Empirical evaluations demonstrate that TERA achieves state-of-the-art predictive accuracy while operating orders of magnitude faster than standard derivative GPs. Crucially, both computation time and peak GPU memory remain essentially flat with respect to

d

, enabling highly scalable inference in high-dimensional spaces.

Hyunseok Seung, Matthias Katzfuss

Department of Statistics · University of Wisconsin–Madison · Madison, WI, 53706

Gaussian Process Surrogate Gradient Gradient Conditional Distribution

Ask PDF

Sep 15, 2026cs.LG

Online Gradient Computation for Warping Gaussian Process Transformations

Warped Gaussian processes (GPs) handle non-Gaussian observations by mapping them into a latent standard GP via a parametric transformation called warping. Existing streaming variants, however, either optimize the warping parameters periodically or sacrifice analytical tractability for a higher model capacity. To bridge this gap, we show that the gradient of the instantaneous negative log-likelihood of a warped GP admits an exact recursive computation. Based on this result, we propose a novel online method for warped GPs that jointly updates the latent GP moments and optimizes the warping parameters.

Emilio Ruiz-Moreno, Konstantinos Slavakis, Baltasar Beferull-Lozano

SIGIPRO Department, Simula Metropolitan Center for Digital Engineering, Oslo, Norway · SURE-AI Center, Simula Research Laboratory, Oslo, Norway · Institute of Science Tokyo, Yokohama, Japan

Gaussian Process Wasserstein Gradient Flows Gradient Streaming

Ask PDF

May 20, 2026stat.ML

Conditioning Gaussian Processes on Almost Anything

Gaussian processes (GPs) offer a principled probabilistic model over functions, but exact inference is restricted to the linear-Gaussian regime. We establish an explicit equivalence between GPs and a class of linear diffusion models, recasting predictive sampling as an ODE with closed-form Gaussian dynamics and a likelihood-dependent guidance term that admits a simple Monte Carlo approximation. In the linear-Gaussian setting, we recover standard GP conditioning exactly; beyond conjugacy, the same machinery handles any conditioning statement admitting point-wise likelihood evaluation -- including non-linear physics, and, for the first time, natural language via large language models. Whitening isolates the irreducible non-Gaussian dynamics, minimising Wasserstein-2 transport cost and eliminating numerical stiffness. The result is a general-purpose GP inference scheme requiring no bespoke derivations. Together, these results provide a general mechanism for incorporating the full richness of real-world knowledge as conditioning information, opening a new frontier for the probabilistic modelling of real-world problems.

Henry Moss, Lachlan Astfalck, Thomas Cowperthwaite +5

Lancaster University · University of New South Wales · University of Cambridge +1

Gaussian Process Diffusion Models

Ask PDF