Vanilla-SPDE Exchange in Spatio-temporal GPs

• Vanilla-SPDE Exchange is a hybrid inference method that leverages the equivalence between kernel-based Gaussian process regression and SPDE state-space formulations for separable spatio-temporal models.
• The method employs static GP regression at a forecast boundary and then transitions to dynamic SPDE propagation to mitigate the cubic costs of dense prediction grids.
• By integrating these approaches, Vanilla-SPDE Exchange achieves significant computational and memory savings while ensuring exact inference under Matérn-3/2 kernel conditions.

Searching arXiv for the exact term and related papers to ground the article in current sources. Vanilla-SPDE Exchange most directly denotes the hybrid exact inference scheme introduced for spatio-temporal Gaussian processes in "Dynamic Gaussian Processes and the Vanilla-SPDE Exchange" (Zhang et al., 30 Jun 2026). In that usage, it exploits an exact equivalence between the standard kernel formulation of Gaussian process regression and an SPDE/state-space formulation, performs static GP regression only at a forecast boundary, and then switches to the SPDE representation for forward propagation and sampling. The construction is designed for exact inference with separable spatio-temporal kernels, especially when observation locations and prediction locations are largely disjoint and when dense prediction grids make standard GP sampling prohibitively expensive (Zhang et al., 30 Jun 2026).

1. Formal setting and motivation

The basic inferential setting is a spatio-temporal Gaussian process on a domain $X=S\times T$, with noisy observations

$$y_i = f(\boldsymbol x_i) + \varepsilon_i,\qquad \boldsymbol x_i=(s_i,t_i),\qquad \varepsilon_i\sim N(0,\sigma_{\text{obs}}^2),$$

and prior $f\sim \mathrm{GP}(\mu,k)$. For test points $\boldsymbol x^*\in X^M$, standard GP regression yields

$$f(\boldsymbol x^*)\mid \mathcal D \sim N\!\big(\mu_{y^*\mid\mathcal D},\Sigma_{y^*\mid\mathcal D}\big),$$

with

$$\mu_{y^*\mid \mathcal D} = \mu(\boldsymbol x^*) + K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}(y-\mu(y)),$$

$$\Sigma_{y^*\mid \mathcal D} = K_{**} - K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}K_*.$$

Using standard dense linear algebra, factorising or inverting $K$ is $\mathcal O(N^3)$; regression scales as $\mathcal O(N^3 + N^2M + NM^2)$; and drawing $$y_i = f(\boldsymbol x_i) + \varepsilon_i,\qquad \boldsymbol x_i=(s_i,t_i),\qquad \varepsilon_i\sim N(0,\sigma_{\text{obs}}^2),$$0 posterior samples adds $$y_i = f(\boldsymbol x_i) + \varepsilon_i,\qquad \boldsymbol x_i=(s_i,t_i),\qquad \varepsilon_i\sim N(0,\sigma_{\text{obs}}^2),$$1 (Zhang et al., 30 Jun 2026).

These costs are particularly restrictive in spatio-temporal problems with dense prediction grids, where $$y_i = f(\boldsymbol x_i) + \varepsilon_i,\qquad \boldsymbol x_i=(s_i,t_i),\qquad \varepsilon_i\sim N(0,\sigma_{\text{obs}}^2),$$2. In that regime, the full cost for regression plus $$y_i = f(\boldsymbol x_i) + \varepsilon_i,\qquad \boldsymbol x_i=(s_i,t_i),\qquad \varepsilon_i\sim N(0,\sigma_{\text{obs}}^2),$$3 samples becomes

$$y_i = f(\boldsymbol x_i) + \varepsilon_i,\qquad \boldsymbol x_i=(s_i,t_i),\qquad \varepsilon_i\sim N(0,\sigma_{\text{obs}}^2),$$4

which is prohibitive when $$y_i = f(\boldsymbol x_i) + \varepsilon_i,\qquad \boldsymbol x_i=(s_i,t_i),\qquad \varepsilon_i\sim N(0,\sigma_{\text{obs}}^2),$$5, $$y_i = f(\boldsymbol x_i) + \varepsilon_i,\qquad \boldsymbol x_i=(s_i,t_i),\qquad \varepsilon_i\sim N(0,\sigma_{\text{obs}}^2),$$6, or $$y_i = f(\boldsymbol x_i) + \varepsilon_i,\qquad \boldsymbol x_i=(s_i,t_i),\qquad \varepsilon_i\sim N(0,\sigma_{\text{obs}}^2),$$7 are large (Zhang et al., 30 Jun 2026). State-space SPDE formulations remove the cubic dependence on the number of time points, but exact inference remains cubic in the spatial state dimension, and that cost worsens when observation and prediction locations are disjoint because the SPDE state must include all unique spatial locations involved in either set (Zhang et al., 30 Jun 2026).

Vanilla-SPDE Exchange is the response to that asymmetry. It uses the static GP formulation where it is advantageous—regression against scattered observation locations—and the SPDE formulation where it is advantageous—forward forecasting and sampling over long temporal grids (Zhang et al., 30 Jun 2026).

2. Static GP and SPDE/state-space equivalence

The method relies on the established exact equivalence between kernel GP models and state-space/SPDE representations for temporal Matérn kernels with half-integer smoothness. For a one-dimensional stationary Matérn GP with smoothness $$y_i = f(\boldsymbol x_i) + \varepsilon_i,\qquad \boldsymbol x_i=(s_i,t_i),\qquad \varepsilon_i\sim N(0,\sigma_{\text{obs}}^2),$$8, the process can be represented as the stationary solution of an SDE

$$y_i = f(\boldsymbol x_i) + \varepsilon_i,\qquad \boldsymbol x_i=(s_i,t_i),\qquad \varepsilon_i\sim N(0,\sigma_{\text{obs}}^2),$$9

For Matérn-$f\sim \mathrm{GP}(\mu,k)$0, with $f\sim \mathrm{GP}(\mu,k)$1, the corresponding two-dimensional state $f\sim \mathrm{GP}(\mu,k)$2 satisfies

$f\sim \mathrm{GP}(\mu,k)$3

Its exact discrete-time transition is

$f\sim \mathrm{GP}(\mu,k)$4

with

$f\sim \mathrm{GP}(\mu,k)$5

(Zhang et al., 30 Jun 2026).

For separable spatio-temporal kernels

$f\sim \mathrm{GP}(\mu,k)$6

with $f\sim \mathrm{GP}(\mu,k)$7 Matérn-$f\sim \mathrm{GP}(\mu,k)$8, the state-space extension takes the form

$f\sim \mathrm{GP}(\mu,k)$9

and, after time discretization,

$\boldsymbol x^*\in X^M$0

The key property is exactness: the finite-dimensional distributions of $\boldsymbol x^*\in X^M$1 under the SPDE/state-space model coincide with those of the kernel GP with covariance $\boldsymbol x^*\in X^M$2 (Zhang et al., 30 Jun 2026).

This equivalence sits within a broader SPDE-as-GP tradition. General SPDE-based stationary random field constructions recover Matérn and Markovian models through the spectral relation $\boldsymbol x^*\in X^M$3 (Vergara et al., 2018). Graph representations then extend the Matérn SPDE construction to unstructured point clouds through sparse GMRF precisions built from graph Laplacians (Sanz-Alonso et al., 2020). Vanilla-SPDE Exchange specializes that general perspective to exact spatio-temporal GP inference with an explicit static-to-dynamic handoff.

3. Mechanics of the exchange

For Matérn-$\boldsymbol x^*\in X^M$4, the SPDE state at each spatial location is

$\boldsymbol x^*\in X^M$5

Because Gaussian processes are closed under linear operators, this induces an augmented vector-valued GP with block kernel

$\boldsymbol x^*\in X^M$6

Vanilla-SPDE Exchange uses this augmented kernel to infer the posterior state at a single forecast boundary time $\boldsymbol x^*\in X^M$7 on the prediction spatial set $\boldsymbol x^*\in X^M$8. Denoting

$\boldsymbol x^*\in X^M$9

static GP regression yields

$$f(\boldsymbol x^*)\mid \mathcal D \sim N\!\big(\mu_{y^*\mid\mathcal D},\Sigma_{y^*\mid\mathcal D}\big),$$0

This posterior is then reinterpreted as the filtering distribution of the SPDE state at time $$f(\boldsymbol x^*)\mid \mathcal D \sim N\!\big(\mu_{y^*\mid\mathcal D},\Sigma_{y^*\mid\mathcal D}\big),$$1 on the prediction locations (Zhang et al., 30 Jun 2026).

The exchange occurs at that point. Once $$f(\boldsymbol x^*)\mid \mathcal D \sim N\!\big(\mu_{y^*\mid\mathcal D},\Sigma_{y^*\mid\mathcal D}\big),$$2 and $$f(\boldsymbol x^*)\mid \mathcal D \sim N\!\big(\mu_{y^*\mid\mathcal D},\Sigma_{y^*\mid\mathcal D}\big),$$3 are obtained, the explicit kernel representation is discarded and the state is propagated using the SPDE/state-space dynamics. For one prediction step,

$$f(\boldsymbol x^*)\mid \mathcal D \sim N\!\big(\mu_{y^*\mid\mathcal D},\Sigma_{y^*\mid\mathcal D}\big),$$4

and for any future time $$f(\boldsymbol x^*)\mid \mathcal D \sim N\!\big(\mu_{y^*\mid\mathcal D},\Sigma_{y^*\mid\mathcal D}\big),$$5,

$$f(\boldsymbol x^*)\mid \mathcal D \sim N\!\big(\mu_{y^*\mid\mathcal D},\Sigma_{y^*\mid\mathcal D}\big),$$6

The physical process is recovered with

$$f(\boldsymbol x^*)\mid \mathcal D \sim N\!\big(\mu_{y^*\mid\mathcal D},\Sigma_{y^*\mid\mathcal D}\big),$$7

and posterior samples are drawn as

$$f(\boldsymbol x^*)\mid \mathcal D \sim N\!\big(\mu_{y^*\mid\mathcal D},\Sigma_{y^*\mid\mathcal D}\big),$$8

(Zhang et al., 30 Jun 2026).

In operational terms, the method consists of four stages: augment the GP to the SPDE state, regress statically only at the forecast boundary, exchange to the SPDE representation, and propagate forward by Kalman prediction on the prediction grid. The construction is exact, not approximate, within the model class where the SPDE representation is exact (Zhang et al., 30 Jun 2026).

4. Computational profile and empirical behavior

Let $$f(\boldsymbol x^*)\mid \mathcal D \sim N\!\big(\mu_{y^*\mid\mathcal D},\Sigma_{y^*\mid\mathcal D}\big),$$9 be the number of observations, $$\mu_{y^*\mid \mathcal D} = \mu(\boldsymbol x^*) + K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}(y-\mu(y)),$$0 and $$\mu_{y^*\mid \mathcal D} = \mu(\boldsymbol x^*) + K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}(y-\mu(y)),$$1 the numbers of unique spatial and temporal observation locations, $$\mu_{y^*\mid \mathcal D} = \mu(\boldsymbol x^*) + K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}(y-\mu(y)),$$2 and $$\mu_{y^*\mid \mathcal D} = \mu(\boldsymbol x^*) + K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}(y-\mu(y)),$$3 the numbers of unique spatial and temporal prediction locations, $$\mu_{y^*\mid \mathcal D} = \mu(\boldsymbol x^*) + K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}(y-\mu(y)),$$4 the number of unique spatial locations across data and prediction, and $$\mu_{y^*\mid \mathcal D} = \mu(\boldsymbol x^*) + K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}(y-\mu(y)),$$5 the number of posterior samples. The computational comparison is as follows (Zhang et al., 30 Jun 2026):

Method	Regression cost $$\mu_{y^*\mid \mathcal D} = \mu(\boldsymbol x^*) + K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}(y-\mu(y)),$$6	Total cost $$\mu_{y^*\mid \mathcal D} = \mu(\boldsymbol x^*) + K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}(y-\mu(y)),$$7
Vanilla	$$\mu_{y^*\mid \mathcal D} = \mu(\boldsymbol x^*) + K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}(y-\mu(y)),$$8	$$\mu_{y^*\mid \mathcal D} = \mu(\boldsymbol x^*) + K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}(y-\mu(y)),$$9
SPDE	$$\Sigma_{y^*\mid \mathcal D} = K_{**} - K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}K_*.$$0	$$\Sigma_{y^*\mid \mathcal D} = K_{**} - K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}K_*.$$1
VaSE	$$\Sigma_{y^*\mid \mathcal D} = K_{**} - K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}K_*.$$2	$$\Sigma_{y^*\mid \mathcal D} = K_{**} - K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}K_*.$$3

The decisive feature is that VaSE keeps the regression stage static but only over the observation set together with one prediction time slice. It therefore avoids the cubic dependence on $$\Sigma_{y^*\mid \mathcal D} = K_{**} - K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}K_*.$$4 that dominates vanilla GP sampling, while also avoiding the $$\Sigma_{y^*\mid \mathcal D} = K_{**} - K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}K_*.$$5 dependence that penalizes pure SPDE formulations when observation and prediction locations are disjoint. Its sampling cost matches the SPDE formulation, $$\Sigma_{y^*\mid \mathcal D} = K_{**} - K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}K_*.$$6, which is linear in the number of prediction time points (Zhang et al., 30 Jun 2026).

The reported experiments use a synthetic spatio-temporal GP on a two-dimensional spatial domain, with a 2D RBF spatial kernel, a Matérn-$$\Sigma_{y^*\mid \mathcal D} = K_{**} - K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}K_*.$$7 temporal kernel, a spatial prediction grid of $$\Sigma_{y^*\mid \mathcal D} = K_{**} - K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}K_*.$$8 points, and $$\Sigma_{y^*\mid \mathcal D} = K_{**} - K_*^\top (K+\sigma_{\text{obs}}^2 I_N)^{-1}K_*.$$9 time points on $K$0 at increment $K$1. Observations are collected on $K$2 at increment $K$3, with total $K$4 ranging from 400 to 2000 under both overlapping and disjoint observation-location scenarios. The task is to draw $K$5 posterior samples for $K$6 over the full spatial grid (Zhang et al., 30 Jun 2026).

The numerical findings align with the complexity analysis. VaSE is consistently faster than vanilla GP. In the disjoint setting, pure SPDE inference may initially be competitive, but its walltime overtakes vanilla GP as $K$7 grows because its regression cost scales with $K$8, whereas VaSE depends on $K$9 and $\mathcal O(N^3)$0 only. Peak memory usage is highest for vanilla GP, because it stores dense Gram matrices over the full prediction grid; SPDE and VaSE have similar, significantly lower memory usage (Zhang et al., 30 Jun 2026).

5. Broader interpretations in adjacent SPDE literature

The phrase itself is not standardized across arXiv. A plausible broader reading is that it names an exchange between a baseline or “vanilla” SPDE representation and a second representation—conditioning, fluctuation, pricing, or numerical structure—that is more convenient for a particular task.

In "SPDE bridges with observation noise and their spatial approximation" (Nunno et al., 2021), the vanilla object is the unconstrained linear parabolic SPDE

$\mathcal O(N^3)$1

and the exchange is from that original mild solution to an SPDE bridge conditioned on noisy final-time data $\mathcal O(N^3)$2. The bridge is represented explicitly through Gaussian conditioning operators, and the paper shows that, for sufficiently rough observation noise, the spatial discretization error of the bridge is essentially the same as that of the corresponding discretization of the original SPDE (Nunno et al., 2021).

In "SPDE limit of the global fluctuations in rank-based models" (Kolli et al., 2016), the phrase can be read in yet another sense: a large interacting particle system with rank-based coefficients has a deterministic hydrodynamic limit described by a porous medium type PDE, while its global fluctuations converge to a linear parabolic SPDE with additive noise. The limiting field is therefore a canonical fluctuation SPDE obtained by exchanging a microscopic interacting system for a Gaussian macroscopic fluctuation equation (Kolli et al., 2016).

In derivative pricing, "On Stochastic Partial Differential Equations and their applications to Derivative Pricing through a conditional Feynman-Kac formula" (Das et al., 2021) establishes an exchange between standard risk-neutral valuation and a backward SPDE representation. The inner conditional expectation in an iterated pricing formula solves a backward SPDE driven by a backward Brownian motion in an enlarged backward filtration, thereby replacing part of a vanilla pricing problem by an SPDE object (Das et al., 2021).

In score-based generative modeling, "Connecting SPDE to SGMs" (Seo, 15 Jan 2025) starts from the vanilla Fokker–Planck PDE associated with an SGM diffusion and introduces an additive-noise SPDE

$\mathcal O(N^3)$3

with existence and uniqueness proved in weighted Sobolev spaces inspired by Ornstein–Uhlenbeck operators. There the exchange is from a deterministic density evolution to a stochastic one, motivated by discretization issues (Seo, 15 Jan 2025).

These uses do not define a single universal doctrine. They do, however, suggest a common pattern: a “vanilla-SPDE exchange” often means replacing a baseline deterministic or kernel-based description by an SPDE description when the latter yields better conditional structure, better scaling, or a more transparent limiting object.

6. Limitations and prospective developments

As presented in its exact GP-inference form, Vanilla-SPDE Exchange is model-class specific. The temporal kernel must be Matérn with half-integer smoothness, so that an exact SDE/SPDE representation exists. The method is formulated for Gaussian observation noise and linear observation models. Its computational savings are strongest when prediction times are numerous and observation and prediction locations are largely disjoint; it still retains quadratic or cubic dependence on the number of spatial prediction locations $\mathcal O(N^3)$4, so extremely large spatial grids may require additional sparsity or low-rank structure (Zhang et al., 30 Jun 2026).

The method is also modular rather than exhaustive. Hyperparameter learning is treated as external: one may optimize via static GP marginal likelihood or use SPDE-based likelihoods where appropriate. The exchange itself is agnostic to that choice. The paper explicitly notes that VaSE is orthogonal to sparse GP approximations, so inducing-point, Vecchia, or structured-kernel-interpolation methods could be inserted into the static regression stage, while approximate state-space representations could replace exact Matérn SPDEs when exact mappings are unavailable (Zhang et al., 30 Jun 2026).

A broader implication is that VaSE is best viewed not as an isolated algorithmic trick but as a specific member of a larger SPDE equivalence program. In that program, the same Gaussian law may be realized as a kernel GP, a state-space model, a mild SPDE solution, or a sparse GMRF approximation on graphs. Vanilla-SPDE Exchange identifies a particularly useful point of transition between those realizations: static regression at a forecast boundary, followed by dynamic propagation in the SPDE state. Within the exact Gaussian setting, this produces exact posterior samples up to numerical precision and clarifies where the cubic costs of GP inference can be avoided without giving up exactness (Zhang et al., 30 Jun 2026).