---
title: Spatiotemporal Gaussian Process Regression
url: https://www.emergentmind.com/topics/spatiotemporal-gaussian-process-regression-gpr
type: topic
---

# Spatiotemporal Gaussian Process Regression

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 1502.03042", max_results=5)
arxiv_search(query="spatiotemporal Gaussian Process Regression 1502.03042", 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 [1611.02978, 1502.03042, 2209.12565, 2603.21144, 2502.02450].

## 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 \(K=K_{XX}+\sigma_n^2 I\). 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 [1611.02978].

When the kernel is sufficiently differentiable, posterior inference extends to derivatives by differentiating the kernel. For profile inference, one may write
\[
m_*^{(x)}(x_*)=\mu'(x_*) + k_*^{(x)\top}(K+\Sigma_\varepsilon)^{-1}(y-\mu(X)),
\]
with analogous formulas for temporal and mixed derivatives. This is central in settings where gradients, rather than levels alone, are the inferential target [2507.05067].

## 2. Kernel design and dependence structure

A common construction is the separable spatiotemporal kernel
\[
k((x,t),(x',t')) = k_s(x,x')\,k_t(t,t'),
\]
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 [2209.12565, 1611.02978].

Nonseparable structure is introduced when the covariance depends jointly on space and time. One example is a composite ambulance-demand kernel
\[
k_{s,t}(x,x') = k_{t,m}(t,t') + k_{s,m}(s,s') + k_{t,i}(t,t')\,k_{s,i}(s,s'),
\]
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 \(S(\omega_s,\omega_t)\) depends jointly on spatial and temporal frequencies rather than factorizing as \(S_s(\omega_s)S_t(\omega_t)\) [1806.10873, 1502.03042].

A central limitation of classical separability appears in the temporal evolution of spatial dependence (TESD), defined as the zero time-lag spatial covariance
\[
C_{f\mid t}(s,s') = \mathrm{Cov}(f(s,t),f(s',t)).
\]
Under a separable kernel, TESD reduces to \(K_s(s,s')K_t(t,t)\), so time only rescales a fixed spatial pattern. The generalized spatiotemporal GP (gSTGP) addresses this by making the spatial Mercer eigenvalues time-dependent:
\[
K_s^{(t)}(s,s')=\sum_{i=1}^\infty \lambda_i^2(t)\phi_i(s)\phi_i(s').
\]
This yields a non-stationary, non-separable covariance targeted specifically at changing spatial dependence over time [1901.04030].

Anisotropy is treated in several equivalent ways. In Euclidean settings it can be encoded by axis-specific ranges, by inserting a spatial metric \(A\) into the phase of a spectral basis, or by using an elliptical norm \(\|\omega\|_M^2=\omega^\top M\omega\). In spatiotemporal models, different spatial axes may have distinct ranges \((\rho_1,\rho_2,\ldots)\), while the temporal axis has its own \(\rho_t\) [1502.03042].

## 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 \(\{\omega_k\}_{k=1}^m\), basis functions \(\phi_k(x)=\exp(i\omega_k^\top x)/\sqrt{n}\), and a shared complex Gaussian spectral vector \(Y(\omega)\), the latent field is
\[
f(x)=\sum_{k=1}^m g^{1/2}(\omega_k)Y(\omega_k)\phi_k(x)=\sum_{k=1}^m w_k\phi_k(x).
\]
For inputs \(X\), the covariance takes the form
\[
K_{XX}=Q_XGQ_X^*+\sigma^2 I.
\]
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 \(Q\) and \(Q^*\) via FFTs. The resulting weight-space posterior avoids \(n\times n\) inversion, with cost scaling as \(O(nm+m^3)\), while lattice and sparse-FFT variants give \(O(n\log n)\) or \(O(m\log n)\) transforms [1502.03042].

A different scalability route is the state-space realization of separable spatiotemporal kernels with rational temporal power spectral density. In that setting,
\[
\Sigma = K_t\otimes K_s + \sigma^2 I_{NM},
\]
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 \(t_k\) is a sufficient statistic to compute the minimum variance estimate of the process at any \(t\ge t_k\) over the entire spatial domain. Another exploits the SVD of \(K_s\) so that transformed innovation covariances become diagonal and transformed state covariances remain block-diagonal, reducing complexity from \(O(NM^3)\) to \(O(M^3+NM^2)\) under a fixed spatial grid and separable kernel [1705.01485, 2209.12565].

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 \(O(Bm^3)\), where \(B\) is the number of blocks and \(m\) a typical block size. Random projection methods instead project all data onto a lower-dimensional subspace through a sketch matrix \(\Phi\), approximating the covariance by
\[
Q^{RP}_{f,f}=K\Phi^\top(\Phi K\Phi^\top)^{-1}\Phi K.
\]
This yields \(O(n^2m+m^3)\) training and \(O(nm)\) prediction, and the paper derives a KL-divergence bound from the Frobenius approximation error [1511.00054, 1106.5779].

## 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 \(Y(\omega)\), while the local spectral modulation varies with location through stick-breaking weights. The process remains jointly Gaussian, and cross-component covariance stays nonzero:
\[
\mathrm{Cov}(Z_{x_1},Z_{x_2}\mid C_{x_1}=k_1,C_{x_2}=k_2)
=Q_{x_1}G_{k_1}^{1/2}G_{k_2}^{1/2}Q_{x_2}^*+\sigma^2 1[x_1=x_2].
\]
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 [1502.03042].

The gSTGP framework generalizes nonstationary space–time covariance in a different way. Rather than treating \((s,t)\) 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 [1901.04030].

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 \(M\), the covariance expansion is
\[
k(x,x',t,t')=\sum_{\ell=0}^{\infty}\sum_{m=1}^{N_\ell} S_\ell(t,t')\,\phi_{\ell m}(x)\phi_{\ell m}(x').
\]
The coefficient processes are Gaussian with diagonal cross-mode covariance, and truncation at \(\ell\le L(T)\) gives exact dimension reduction in the spectral domain. The paper studies logarithmic truncation \(L(T)=\lfloor \log T\rfloor\) and power-law truncation \(L(T)=\lfloor T^\varrho\rfloor\), with Empirical Bayes updates performed sequentially in time [2603.21144].

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 \(z=(x,t)\). 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 [2110.03712].

## 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
\[
\rho(x)=\exp(f(x)),\qquad y_i\mid f_i \sim \mathrm{Poisson}(\exp(o_i+f_i)),
\]
where \(o_i=\log(A_i\tau_i)\) is an exposure offset. With inducing points and stochastic variational inference, the ELBO contains the exact expectation
\[
\mathbb{E}[\log p(y_i\mid f_i)]
= y_i(o_i+\mu_i)-\exp(o_i)\exp(\mu_i+0.5\sigma_i^2)-\log(y_i!),
\]
which makes Poisson inference computationally practical [1806.10873].

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
\[
w_{\mathrm{IMQ}}(x,y)=\beta\left(1+\frac{(y-\gamma(x))^2}{c(x)^2}\right)^{-1/2}.
\]
The resulting filtering and smoothing updates remain Gaussian and conjugate. Adaptive choices
\[
\gamma_k = Hm_{k\mid k-1},\qquad c_k^2 = \mathrm{diag}(S_k),\qquad \beta=\sigma/\sqrt{2}
\]
remove manual tuning, and the paper proves bounded posterior influence under these weights [2502.02450].

A different prior-based strategy appears in spatio-temporal structured sparse regression. There the latent signal is sparse, with spike-and-slab variables \(\omega_{it}\), a spatial GP prior on support logits \(\gamma_t\), and a temporal GP prior on evolving spatial means \(\mu_t\):
\[
p(\gamma_t)=\mathcal{N}(\mu_t,\Sigma_0),\qquad \mu_t\sim \mathcal{N}(\mu_{t-1},W).
\]
Expectation propagation yields a Gaussian–Bernoulli approximation for \((x_t,\omega_t,\gamma_t,\mu_t)\), with \(O(N^3T)\) per iteration for the hierarchical model, compared with \(O(N^3T^3)\) for the one-level GP baseline [1807.05561].

Physics-informed spatiotemporal GPR uses linear operators to encode mechanistic knowledge directly in the kernel. In the systems-biology example,
\[
-D\frac{\partial^2 y}{\partial x^2}+\alpha\frac{\partial y}{\partial t}+\beta y = f,
\]
so if \(y\) has GP prior covariance \(k_{yy}\), then
\[
k_{ff}=\mathcal{L}_{x,t}\mathcal{L}_{x',t'}k_{yy},\qquad
k_{fy}=\mathcal{L}_{x,t}k_{yy},\qquad
k_{yf}=\mathcal{L}_{x',t'}k_{yy}.
\]
This yields a multi-output GP over protein concentration and mRNA expression without explicitly solving the PDE, and posterior inference over \((D,\alpha,\beta)\) is carried out by Hybrid Monte Carlo [1610.05163].

## 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
\[
\mathcal{F}[\bar m](\xi)\approx H(\xi)\mathcal{F}[f](\xi),\qquad
H(\xi)=\frac{\mathcal{F}[\kappa](\xi)}{\sigma_\varepsilon^2\Delta x+\mathcal{F}[\kappa](\xi)}.
\]
This defines an effective cutoff frequency \(\xi_*\) through \(H(\xi_*)=1/2\), and for Matérn kernels yields an explicit expression in terms of the signal-to-noise rate \(S=\sigma^2\ell/(\sigma_\varepsilon^2\Delta x)\). The paper also introduces
\[
N_{\mathrm{eff}}(x)=\frac{\left(\sum_i s_i\right)^2}{\sum_i s_i^2},
\]
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, \(N_{\mathrm{eff}}\) near 1–2 indicates high sensitivity to a small number of measurements [2507.05067].

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 [1611.02978].

## 7. Applications and empirical comparisons

The empirical literature is dominated by large-scale environmental, biomedical, and operational applications. In stationary simulations with \(n=2500\) irregular locations, FGP achieved the most accurate parameter estimates and fastest time for Matérn \((\kappa=1.5)\) and squared exponential kernels; for the squared exponential case, the reduced dimension \(m\approx 0.13n\) cut runtime by approximately \(6\times\) with negligible loss. In a nonstationary localized squared exponential experiment, NS-FGP converged to 3 clusters aligned with spatially varying ranges and achieved RMSE \(=1.81\) and MAD \(=1.62\), outperforming stationary FGP (RMSE \(=4.96\)) and GPP (RMSE \(=10.20\)). On the NARCCAP temperature dataset with \(1{,}343{,}752\) daily temperatures over 92 days and \(1{,}075{,}000\) training points, spectral truncation to approximately 5% of \(n\) used about 53,000 frequencies, and 30,000 MCMC iterations completed in approximately 27 hours per run; NS-FGP achieved RMSE \(=1.13\) and MAD \(=0.65\), versus stationary FGP RMSE \(=2.75\) and MAD \(=1.26\) [1502.03042].

State-space implementations show a different empirical profile. On Colorado precipitation, the proposed \(O(M^3+NM^2)\) implementation reduced average MLM runtime to \(0.65\) seconds versus \(197.37\) seconds for the prior \(O(NM^3)\) 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 \(60.34\) for MLM, \(60.30\) for GCV, and \(52.40\) 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 [2209.12565].

For sparse irregular time series, the secondary autoregressive evaluation reported MAPE-AR values below 1% across all sparsity levels: at 3% sparsity, OU \(0.0049\) and Fractional \(0.0042\); at 5%, OU \(0.0051\) and Fractional \(0.0056\); at 7%, OU \(0.0067\) and Fractional \(0.0057\). 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 [1611.02978].

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 \(=0.174\) and log-likelihood \(=-6.83\times 10^4\), versus MEDIC MAE \(=0.181\) and log-likelihood \(=-7.96\times 10^4\). 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 [1806.10873].

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 \(K=100\) trials, TESD MSE was \(5.1\times 10^{-3}\), versus \(0.11\)–\(0.15\) for the comparison methods; with \(K=1000\), TESD MSE was \(3.9\times 10^{-4}\), versus approximately \(0.10\)–\(0.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 [1901.04030].

Physics-informed and manifold-aware models show similarly domain-specific behavior. In the systems-biology PDE example, simulated data with true \(D=\alpha=\beta=1\) yielded posterior means \(D\approx 0.993\), \(\alpha\approx 1.002\), and \(\beta\approx 1.002\) after 7000 HMC iterations, while the Drosophila application gave \(D\approx 0.017\), \(\alpha\approx 13.04\), and \(\beta\approx 0.65\) after 10000 iterations. On \(\mathbb{S}^2\), the manifold FGPR study reported that logarithmic truncation \(L(T)=\lfloor\log T\rfloor\) was well-suited under sparse spatial sampling, while power-law truncation improved accuracy in the more singular, long-range dependent subfamily for small \(T\) [1610.05163, 2603.21144].

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 [1901.04030, 1502.03042, 1705.01485].

Source: https://www.emergentmind.com/topics/spatiotemporal-gaussian-process-regression-gpr