Spatiotemporal Gaussian Process Regression

• Spatiotemporal Gaussian Process Regression is a probabilistic framework that models continuous latent functions over space and time using Gaussian process priors and tailored covariance kernels.
• It employs both separable and nonseparable kernel designs to capture spatial, temporal, and space–time dependencies, enhancing interpretability and adaptability across various applications.
• Efficient computational methods, including spectral reduction, state-space modeling, and sparse approximations, enable scalable inference and reliable uncertainty quantification for large datasets.

Searching arXiv for the cited papers and closely related spatiotemporal GPR work. arxiv_search(query="spatiotemporal Gaussian Process Regression functional Gaussian process Bayesian nonparametric analysis (Duan et al., 2015)", max_results=5) arxiv_search(query="spatiotemporal Gaussian Process Regression (Duan et al., 2015)", max_results=10) Spatiotemporal Gaussian Process Regression (GPR) is a probabilistic framework for modeling continuous latent functions over space and time, typically written as $f(s,t)$ or $f(x,t)$, from discrete and noisy observations. It places a Gaussian process prior on the latent field, encodes spatial, temporal, and space–time dependence through a covariance kernel, and yields posterior Gaussian predictions with uncertainty quantification. Within this general formulation, the literature spans separable and nonseparable kernels, spectral reduced-rank constructions, state-space and Kalman realizations, manifold-valued formulations, nonstationary covariance learning, robust and non-Gaussian likelihoods, and application-specific models for count data, dynamical systems, and large-scale environmental and biomedical inference (Süzen et al., 2016, Duan et al., 2015, Zhang et al., 2022, Ruiz-Medina et al., 22 Mar 2026, Laplante et al., 4 Feb 2025).

1. Core probabilistic formulation

A standard setup uses observations

$$y_i = f(x_i,t_i) + \varepsilon_i,\qquad \varepsilon_i \sim \mathcal{N}(0,\sigma_n^2),$$

with irregularly sampled spatiotemporal inputs $D=\{(x_i,t_i,y_i)\}_{i=1}^n$. Writing $z=(x,t)$, one places a prior

$f(\cdot)\sim \mathrm{GP}(m(\cdot),k(\cdot,\cdot)).$

For test inputs $z_*$, the posterior predictive distribution is Gaussian with mean

$$\mu_* = m(z_*) + K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}(y-m(X)),$$

and variance

$$\sigma^2_* = k(z_*,z_*) - K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}K_{X*}.$$

Hyperparameters are commonly learned by maximizing the log marginal likelihood

$$\log p(y\mid \theta) = -\frac{1}{2}(y-m)^\top K^{-1}(y-m)-\frac{1}{2}\log|K|-\frac{n}{2}\log 2\pi,$$

with $f(x,t)$0. A recurring advantage of this formulation is that irregular sampling in time does not require resampling; sparse, irregular observations are handled directly through the covariance matrix (Süzen et al., 2016).

When the kernel is sufficiently differentiable, posterior inference extends to derivatives by differentiating the kernel. For profile inference, one may write

$f(x,t)$1

with analogous formulas for temporal and mixed derivatives. This is central in settings where gradients, rather than levels alone, are the inferential target (Cao et al., 7 Jul 2025).

2. Kernel design and dependence structure

A common construction is the separable spatiotemporal kernel

$f(x,t)$2

which factorizes spatial and temporal dependence. In practice, this includes spatial squared exponential or Matérn terms multiplied by temporal kernels such as Ornstein–Uhlenbeck, Matérn, periodic, or exponentially decaying components. Separable kernels are attractive because they admit Kronecker structure, state-space realizations, and relatively interpretable length-scales in space and time (Zhang et al., 2022, Süzen et al., 2016).

Nonseparable structure is introduced when the covariance depends jointly on space and time. One example is a composite ambulance-demand kernel

$f(x,t)$3

which adds marginal spatial and temporal terms to a space–time interaction term. Another is the Gneiting class, used as a canonical nonseparable construction in several spatiotemporal settings. In the spectral FGP setting, nonseparability arises when the spatiotemporal spectral density $f(x,t)$4 depends jointly on spatial and temporal frequencies rather than factorizing as $f(x,t)$5 (Nabarro et al., 2018, Duan et al., 2015).

A central limitation of classical separability appears in the temporal evolution of spatial dependence (TESD), defined as the zero time-lag spatial covariance

$f(x,t)$6

Under a separable kernel, TESD reduces to $f(x,t)$7, so time only rescales a fixed spatial pattern. The generalized spatiotemporal GP (gSTGP) addresses this by making the spatial Mercer eigenvalues time-dependent: $f(x,t)$8 This yields a non-stationary, non-separable covariance targeted specifically at changing spatial dependence over time (Lan, 2019).

Anisotropy is treated in several equivalent ways. In Euclidean settings it can be encoded by axis-specific ranges, by inserting a spatial metric $f(x,t)$9 into the phase of a spectral basis, or by using an elliptical norm $$y_i = f(x_i,t_i) + \varepsilon_i,\qquad \varepsilon_i \sim \mathcal{N}(0,\sigma_n^2),$$0. In spatiotemporal models, different spatial axes may have distinct ranges $$y_i = f(x_i,t_i) + \varepsilon_i,\qquad \varepsilon_i \sim \mathcal{N}(0,\sigma_n^2),$$1, while the temporal axis has its own $$y_i = f(x_i,t_i) + \varepsilon_i,\qquad \varepsilon_i \sim \mathcal{N}(0,\sigma_n^2),$$2 (Duan et al., 2015).

3. Spectral, state-space, and large-scale computation

The Functional Gaussian Process (FGP) gives a spectral reduced-rank construction of a Gaussian process by projecting from a fixed, finite set of discrete spectral frequencies into continuous input space. With symmetric frequencies $$y_i = f(x_i,t_i) + \varepsilon_i,\qquad \varepsilon_i \sim \mathcal{N}(0,\sigma_n^2),$$3, basis functions $$y_i = f(x_i,t_i) + \varepsilon_i,\qquad \varepsilon_i \sim \mathcal{N}(0,\sigma_n^2),$$4, and a shared complex Gaussian spectral vector $$y_i = f(x_i,t_i) + \varepsilon_i,\qquad \varepsilon_i \sim \mathcal{N}(0,\sigma_n^2),$$5, the latent field is

$$y_i = f(x_i,t_i) + \varepsilon_i,\qquad \varepsilon_i \sim \mathcal{N}(0,\sigma_n^2),$$6

For inputs $$y_i = f(x_i,t_i) + \varepsilon_i,\qquad \varepsilon_i \sim \mathcal{N}(0,\sigma_n^2),$$7, the covariance takes the form

$$y_i = f(x_i,t_i) + \varepsilon_i,\qquad \varepsilon_i \sim \mathcal{N}(0,\sigma_n^2),$$8

This is a valid stationary Gaussian process by construction, is connected to Bochner’s theorem through a discrete spectral measure, and supports fast multiplication by $$y_i = f(x_i,t_i) + \varepsilon_i,\qquad \varepsilon_i \sim \mathcal{N}(0,\sigma_n^2),$$9 and $D=\{(x_i,t_i,y_i)\}_{i=1}^n$0 via FFTs. The resulting weight-space posterior avoids $D=\{(x_i,t_i,y_i)\}_{i=1}^n$1 inversion, with cost scaling as $D=\{(x_i,t_i,y_i)\}_{i=1}^n$2, while lattice and sparse-FFT variants give $D=\{(x_i,t_i,y_i)\}_{i=1}^n$3 or $D=\{(x_i,t_i,y_i)\}_{i=1}^n$4 transforms (Duan et al., 2015).

A different scalability route is the state-space realization of separable spatiotemporal kernels with rational temporal power spectral density. In that setting,

$D=\{(x_i,t_i,y_i)\}_{i=1}^n$5

and the temporal kernel is realized as a finite-dimensional linear state-space model. Kalman filtering then provides linear-in-time inference. One line of work proves a “Kalman representer theorem”: for fixed observed spatial locations, the filtered state at time $D=\{(x_i,t_i,y_i)\}_{i=1}^n$6 is a sufficient statistic to compute the minimum variance estimate of the process at any $D=\{(x_i,t_i,y_i)\}_{i=1}^n$7 over the entire spatial domain. Another exploits the SVD of $D=\{(x_i,t_i,y_i)\}_{i=1}^n$8 so that transformed innovation covariances become diagonal and transformed state covariances remain block-diagonal, reducing complexity from $D=\{(x_i,t_i,y_i)\}_{i=1}^n$9 to $z=(x,t)$0 under a fixed spatial grid and separable kernel (Todescato et al., 2017, Zhang et al., 2022).

Two further large-data approximations occupy complementary positions. The Gaussian Process Random Field (GPRF) partitions data into local blocks and couples exact local GP marginals through pairwise potentials, yielding a tractable approximation to the full GP marginal likelihood. With sparse local neighborhoods, complexity becomes $z=(x,t)$1, where $z=(x,t)$2 is the number of blocks and $z=(x,t)$3 a typical block size. Random projection methods instead project all data onto a lower-dimensional subspace through a sketch matrix $z=(x,t)$4, approximating the covariance by

$z=(x,t)$5

This yields $z=(x,t)$6 training and $z=(x,t)$7 prediction, and the paper derives a KL-divergence bound from the Frobenius approximation error (Moore et al., 2015, Banerjee et al., 2011).

4. Nonstationary, manifold, and geometry-aware extensions

FGP admits a nonstationary extension, NS-FGP, through a mixture with spectral dependency. Its distinctive feature is that all mixture components share a single spectral vector $z=(x,t)$8, while the local spectral modulation varies with location through stick-breaking weights. The process remains jointly Gaussian, and cross-component covariance stays nonzero: $z=(x,t)$9 This construction is explicitly contrasted with standard mixture models that imply conditional independence across components. In practice it clusters data based on local stationarity while preserving a global Gaussian dependence structure (Duan et al., 2015).

The gSTGP framework generalizes nonstationary space–time covariance in a different way. Rather than treating $f(\cdot)\sim \mathrm{GP}(m(\cdot),k(\cdot,\cdot)).$0 as a single undifferentiated input, it varies the eigenvalues of the spatial kernel over time, producing a quasi Kronecker sum structure. The paper establishes posterior contraction for the covariance operator itself and reports that Model II learns TESD accurately while remaining computationally efficient through sparse block structure and low-dimensional spectral trajectories (Lan, 2019).

On compact Riemannian manifolds, time-adaptive functional Gaussian Process regression uses invariance under isometries to diagonalize the spatial covariance in the Laplace–Beltrami eigenbasis. For a connected, compact two-point homogeneous space $f(\cdot)\sim \mathrm{GP}(m(\cdot),k(\cdot,\cdot)).$1, the covariance expansion is

$f(\cdot)\sim \mathrm{GP}(m(\cdot),k(\cdot,\cdot)).$2

The coefficient processes are Gaussian with diagonal cross-mode covariance, and truncation at $f(\cdot)\sim \mathrm{GP}(m(\cdot),k(\cdot,\cdot)).$3 gives exact dimension reduction in the spectral domain. The paper studies logarithmic truncation $f(\cdot)\sim \mathrm{GP}(m(\cdot),k(\cdot,\cdot)).$4 and power-law truncation $f(\cdot)\sim \mathrm{GP}(m(\cdot),k(\cdot,\cdot)).$5, with Empirical Bayes updates performed sequentially in time (Ruiz-Medina et al., 22 Mar 2026).

Geospatial trajectory modeling supplies a related geometry-aware specialization. Trajectories may be modeled as one-dimensional time-indexed position with separate single-output GPs for longitude and latitude, or as spatiotemporal fields with input $f(\cdot)\sim \mathrm{GP}(m(\cdot),k(\cdot,\cdot)).$6. The chapter emphasizes great-circle distance, local tangent-plane coordinates, derivative GPs for velocity and acceleration, and the use of white kernels to represent known GPS accuracy (Nguyen et al., 2021).

5. Likelihoods, robustness, and structured priors

Although the basic spatiotemporal GPR model is Gaussian, several important variants are driven by non-Gaussian likelihoods or more structured priors. For count data, ambulance demand is modeled as a log-Gaussian Cox process with

$f(\cdot)\sim \mathrm{GP}(m(\cdot),k(\cdot,\cdot)).$7

where $f(\cdot)\sim \mathrm{GP}(m(\cdot),k(\cdot,\cdot)).$8 is an exposure offset. With inducing points and stochastic variational inference, the ELBO contains the exact expectation

$f(\cdot)\sim \mathrm{GP}(m(\cdot),k(\cdot,\cdot)).$9

which makes Poisson inference computationally practical (Nabarro et al., 2018).

Robustness to outliers has been developed in state-space spatiotemporal GPs through the robust and conjugate GP framework. Instead of a heavy-tailed likelihood, it uses generalized Bayes with a weighted score-matching loss and inverse multiquadric weights

$z_*$0

The resulting filtering and smoothing updates remain Gaussian and conjugate. Adaptive choices

$z_*$1

remove manual tuning, and the paper proves bounded posterior influence under these weights (Laplante et al., 4 Feb 2025).

A different prior-based strategy appears in spatio-temporal structured sparse regression. There the latent signal is sparse, with spike-and-slab variables $z_*$2, a spatial GP prior on support logits $z_*$3, and a temporal GP prior on evolving spatial means $z_*$4: $z_*$5 Expectation propagation yields a Gaussian–Bernoulli approximation for $z_*$6, with $z_*$7 per iteration for the hierarchical model, compared with $z_*$8 for the one-level GP baseline (Kuzin et al., 2018).

Physics-informed spatiotemporal GPR uses linear operators to encode mechanistic knowledge directly in the kernel. In the systems-biology example,

$z_*$9

so if $$\mu_* = m(z_*) + K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}(y-m(X)),$$0 has GP prior covariance $$\mu_* = m(z_*) + K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}(y-m(X)),$$1, then

$$\mu_* = m(z_*) + K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}(y-m(X)),$$2

This yields a multi-output GP over protein concentration and mRNA expression without explicitly solving the PDE, and posterior inference over $$\mu_* = m(z_*) + K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}(y-m(X)),$$3 is carried out by Hybrid Monte Carlo (Niu et al., 2016).

6. Resolution diagnostics and evaluation methodology

Recent work has made the smoothing behavior of GPR explicit through a low-pass filter interpretation. For stationary kernels on a regularly sampled lattice and negligible aliasing, the expectation of the GP estimator satisfies

$$\mu_* = m(z_*) + K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}(y-m(X)),$$4

This defines an effective cutoff frequency $$\mu_* = m(z_*) + K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}(y-m(X)),$$5 through $$\mu_* = m(z_*) + K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}(y-m(X)),$$6, and for Matérn kernels yields an explicit expression in terms of the signal-to-noise rate $$\mu_* = m(z_*) + K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}(y-m(X)),$$7. The paper also introduces

$$\mu_* = m(z_*) + K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}(y-m(X)),$$8

an information-theoretic measure of the effective number of data points supporting an inferred value or derivative. In pedestal-profile inference, these quantities quantify the trade-off between over-fitting and over-regularization; for derivatives, $$\mu_* = m(z_*) + K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}(y-m(X)),$$9 near 1–2 indicates high sensitivity to a small number of measurements (Cao et al., 7 Jul 2025).

Evaluation of spatiotemporal GPR under sparse irregular sampling can itself be nontrivial when no dense reference trajectory exists. One practical method is to use the GP posterior to reconstruct a fine-scale regularly spaced series, fit a secondary autoregressive model to that reconstruction, and score one-step-ahead forecasts against the next observed value. In the reported experiments, the secondary model was seasonal ARIMA(1,1,1)(1,1,1), and the resulting MAPE-AR metric was used to compare GP reconstructions of Ornstein–Uhlenbeck and fractional processes across 3%, 5%, and 7% sparsity regimes (Süzen et al., 2016).

7. Applications and empirical comparisons

The empirical literature is dominated by large-scale environmental, biomedical, and operational applications. In stationary simulations with $$\sigma^2_* = k(z_*,z_*) - K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}K_{X*}.$$0 irregular locations, FGP achieved the most accurate parameter estimates and fastest time for Matérn $$\sigma^2_* = k(z_*,z_*) - K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}K_{X*}.$$1 and squared exponential kernels; for the squared exponential case, the reduced dimension $$\sigma^2_* = k(z_*,z_*) - K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}K_{X*}.$$2 cut runtime by approximately $$\sigma^2_* = k(z_*,z_*) - K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}K_{X*}.$$3 with negligible loss. In a nonstationary localized squared exponential experiment, NS-FGP converged to 3 clusters aligned with spatially varying ranges and achieved RMSE $$\sigma^2_* = k(z_*,z_*) - K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}K_{X*}.$$4 and MAD $$\sigma^2_* = k(z_*,z_*) - K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}K_{X*}.$$5, outperforming stationary FGP (RMSE $$\sigma^2_* = k(z_*,z_*) - K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}K_{X*}.$$6) and GPP (RMSE $$\sigma^2_* = k(z_*,z_*) - K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}K_{X*}.$$7). On the NARCCAP temperature dataset with $$\sigma^2_* = k(z_*,z_*) - K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}K_{X*}.$$8 daily temperatures over 92 days and $$\sigma^2_* = k(z_*,z_*) - K_{*X}[K_{XX}+\sigma_n^2 I]^{-1}K_{X*}.$$9 training points, spectral truncation to approximately 5% of $$\log p(y\mid \theta) = -\frac{1}{2}(y-m)^\top K^{-1}(y-m)-\frac{1}{2}\log|K|-\frac{n}{2}\log 2\pi,$$0 used about 53,000 frequencies, and 30,000 MCMC iterations completed in approximately 27 hours per run; NS-FGP achieved RMSE $$\log p(y\mid \theta) = -\frac{1}{2}(y-m)^\top K^{-1}(y-m)-\frac{1}{2}\log|K|-\frac{n}{2}\log 2\pi,$$1 and MAD $$\log p(y\mid \theta) = -\frac{1}{2}(y-m)^\top K^{-1}(y-m)-\frac{1}{2}\log|K|-\frac{n}{2}\log 2\pi,$$2, versus stationary FGP RMSE $$\log p(y\mid \theta) = -\frac{1}{2}(y-m)^\top K^{-1}(y-m)-\frac{1}{2}\log|K|-\frac{n}{2}\log 2\pi,$$3 and MAD $$\log p(y\mid \theta) = -\frac{1}{2}(y-m)^\top K^{-1}(y-m)-\frac{1}{2}\log|K|-\frac{n}{2}\log 2\pi,$$4 (Duan et al., 2015).

State-space implementations show a different empirical profile. On Colorado precipitation, the proposed $$\log p(y\mid \theta) = -\frac{1}{2}(y-m)^\top K^{-1}(y-m)-\frac{1}{2}\log|K|-\frac{n}{2}\log 2\pi,$$5 implementation reduced average MLM runtime to $$\log p(y\mid \theta) = -\frac{1}{2}(y-m)^\top K^{-1}(y-m)-\frac{1}{2}\log|K|-\frac{n}{2}\log 2\pi,$$6 seconds versus $$\log p(y\mid \theta) = -\frac{1}{2}(y-m)^\top K^{-1}(y-m)-\frac{1}{2}\log|K|-\frac{n}{2}\log 2\pi,$$7 seconds for the prior $$\log p(y\mid \theta) = -\frac{1}{2}(y-m)^\top K^{-1}(y-m)-\frac{1}{2}\log|K|-\frac{n}{2}\log 2\pi,$$8 approach, and analogous speedups were reported for GCV and SURE. On GHCN temperature, average fits with the proposed seasonal kernel plus Matérn component were $$\log p(y\mid \theta) = -\frac{1}{2}(y-m)^\top K^{-1}(y-m)-\frac{1}{2}\log|K|-\frac{n}{2}\log 2\pi,$$9 for MLM, $f(x,t)$00 for GCV, and $f(x,t)$01 for SURE, while the periodically decaying baseline gave strongly negative fits. These experiments emphasize that computational structure and kernel design are tightly coupled in large spatiotemporal GPR (Zhang et al., 2022).

For sparse irregular time series, the secondary autoregressive evaluation reported MAPE-AR values below 1% across all sparsity levels: at 3% sparsity, OU $f(x,t)$02 and Fractional $f(x,t)$03; at 5%, OU $f(x,t)$04 and Fractional $f(x,t)$05; at 7%, OU $f(x,t)$06 and Fractional $f(x,t)$07. The reported interpretation was that both GP reconstructions supported effective short-horizon forecasting under sparsity, with fractional dynamics appearing more robust at the most extreme sparsity level (Süzen et al., 2016).

In count-data forecasting, the ambulance-demand spatiotemporal LGCP improved both MAE and inhomogeneous Poisson log-likelihood over MEDIC on Cape Town data: STGP achieved MAE $f(x,t)$08 and log-likelihood $f(x,t)$09, versus MEDIC MAE $f(x,t)$10 and log-likelihood $f(x,t)$11. The paper attributes the particularly large log-likelihood gain to smooth spatial predictions that reduce wasted density in low-demand regions and avoid zero predictions at event locations (Nabarro et al., 2018).

In covariance-learning rather than mean-only prediction, gSTGP Model II was reported to estimate TESD far more accurately than separable, stationary nonseparable, and nonstationary baselines. With $f(x,t)$12 trials, TESD MSE was $f(x,t)$13, versus $f(x,t)$14–$f(x,t)$15 for the comparison methods; with $f(x,t)$16, TESD MSE was $f(x,t)$17, versus approximately $f(x,t)$18–$f(x,t)$19. In longitudinal PET imaging, the learned TESD revealed evolving ROI–POI correlations and decreasing graph connectivity in later stages for MCI and AD groups (Lan, 2019).

Physics-informed and manifold-aware models show similarly domain-specific behavior. In the systems-biology PDE example, simulated data with true $f(x,t)$20 yielded posterior means $f(x,t)$21, $f(x,t)$22, and $f(x,t)$23 after 7000 HMC iterations, while the Drosophila application gave $f(x,t)$24, $f(x,t)$25, and $f(x,t)$26 after 10000 iterations. On $f(x,t)$27, the manifold FGPR study reported that logarithmic truncation $f(x,t)$28 was well-suited under sparse spatial sampling, while power-law truncation improved accuracy in the more singular, long-range dependent subfamily for small $f(x,t)$29 (Niu et al., 2016, Ruiz-Medina et al., 22 Mar 2026).

A common misconception is that nonseparability necessarily requires treating location and time together as a joint variable. Several lines of work reject that premise for distinct reasons: gSTGP introduces time dependence through spatial eigenvalues rather than a generic joint kernel; NS-FGP maintains global Gaussianity through a shared spectral vector rather than conditional independence across local regimes; and state-space methods exploit separability not as a modeling limitation alone, but as a route to exact scalable inference when the temporal spectrum is rational (Lan, 2019, Duan et al., 2015, Todescato et al., 2017).