Neural Conditional Simulation for Complex Spatial Processes

TL;DR

This paper introduces neural conditional simulation (NCS), a neural diffusion-based method for efficient and accurate spatial process simulation conditioned on observations, applicable to complex models where traditional methods are intractable.

Contribution

The paper presents a novel neural diffusion model for conditional spatial simulation that is amortized and flexible, outperforming traditional methods like MCMC in efficiency and accuracy.

Findings

NCS accurately replicates true Gaussian process conditional distributions.

NCS outperforms MCMC in efficiency and accuracy for spatial extremes.

The method is broadly applicable to various spatial process models.

Abstract

A key objective in spatial statistics is to simulate from the distribution of a spatial process at a selection of unobserved locations conditional on observations (i.e., a predictive distribution) to enable spatial prediction and uncertainty quantification. However, exact conditional simulation from this predictive distribution is intractable or inefficient for many spatial process models. In this paper, we propose neural conditional simulation (NCS), a general method for spatial conditional simulation that is based on neural diffusion models. Specifically, using spatial masks, we implement a conditional score-based diffusion model that evolves Gaussian noise into samples from a predictive distribution when given a partially observed spatial field and spatial process parameters as inputs. The diffusion model relies on a neural network that only requires unconditional samples from the spatial process for training. Once trained, the diffusion model is amortized with respect to the observations in the partially observed field, the number and locations of those observations, and the spatial process parameters, and can therefore be used to conditionally simulate from a broad class of predictive distributions without retraining the neural network. We assess the NCS-generated simulations against simulations from the true conditional distribution of a Gaussian process model, and against Markov chain Monte Carlo (MCMC) simulations from a Brown--Resnick process model for spatial extremes. In the latter case, we show that it is more efficient and accurate to conditionally simulate using NCS than classical MCMC techniques implemented in standard software. We conclude that NCS enables efficient and accurate conditional simulation from spatial predictive distributions that are challenging to sample from using traditional methods.

Equations

$$\text{penalty}10000h(\cdot )\sim \text{SpatialProcess}(\mathbf{\theta }),$$

$$\text{penalty}10000{\mathring{\mathbf{x}}}_0:=h(\mathring{\mathcal{S}}),$$

$$\text{penalty}10000{\tilde{\mathbf{x}}}_0\sim p({\tilde{\mathbf{x}}}_0\mid {\mathring{\mathbf{x}}}_0;\mathbf{\theta }),$$

$$\text{penalty}10000{\mathbf{x}}_0\sim p({\mathbf{x}}_0),$$

$$\text{penalty}10000\mathrm{d}{\mathbf{x}}_t=\mathbf{f}({\mathbf{x}}_t,t)\mathrm{d}t+g(t)\mathrm{d}{\mathbf{B}}_t,t\in [0,T],$$

$$\text{penalty}10000{\mathbf{x}}_T\sim p_{\text{ref}}(\mathbf{x}),$$

$$\text{penalty}10000\mathrm{d}{\mathbf{x}}_t=(\mathbf{f}({\mathbf{x}}_t,t)-g(t{)}^2{\nabla }_{{\mathbf{x}}_t}\log (p({\mathbf{x}}_t)))\mathrm{d}t+g(t)\mathrm{d}{\mathbf{R}}_t,t\in [0,T],$$

$${\mathbf{\varphi }}^{\ast }$$

$$=\underset{\mathbf{\varphi }}{\mathrm{argmin}}{\mathbb{E}}_{t\sim U(0,T)}(\lambda (t){\mathbb{E}}_{{\mathbf{x}}_0}{\mathbb{E}}_{{\mathbf{x}}_t\mid {\mathbf{x}}_0}[{\Vert {\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t)-{\nabla }_{{\mathbf{x}}_t}\log (p({\mathbf{x}}_t\mid {\mathbf{x}}_0))\Vert }_2^2\left]\right),$$

$$\text{penalty}10000{\tilde{\mathbf{x}}}_0\sim p({\tilde{\mathbf{x}}}_0\mid {\mathring{\mathbf{x}}}_0),{\mathring{\mathbf{x}}}_0:=h(\mathring{\mathcal{S}}),$$

$$\text{penalty}10000\mathrm{d}{\tilde{\mathbf{x}}}_t=\tilde{\mathbf{f}}({\tilde{\mathbf{x}}}_t,t)\mathrm{d}t+g(t)\mathrm{d}{\tilde{\mathbf{B}}}_t,t\in [0,T],$$

$$\text{penalty}10000{\tilde{\mathbf{x}}}_T\sim p_{\text{ref}}(\tilde{\mathbf{x}}),{\mathring{\mathbf{x}}}_0:=h(\mathring{\mathcal{S}}),$$

$$\text{penalty}10000\mathrm{d}{\tilde{\mathbf{x}}}_t=(\tilde{\mathbf{f}}({\tilde{\mathbf{x}}}_t,t)-g(t{)}^2{\nabla }_{{\tilde{\mathbf{x}}}_t}\log (p({\tilde{\mathbf{x}}}_t\mid {\mathring{\mathbf{x}}}_0)))\mathrm{d}t+g(t)\mathrm{d}{\tilde{\mathbf{R}}}_t,t\in [0,T],$$

$${\mathbf{\varphi }}^{\ast }$$

$$=\underset{\mathbf{\varphi }}{\mathrm{argmin}}{\mathbb{E}}_{t\sim U(0,T)}(\lambda (t){\mathbb{E}}_{\mathring{\mathcal{S}}}{\mathbb{E}}_{{\mathring{\mathbf{x}}}_0}{\mathbb{E}}_{{\tilde{\mathbf{x}}}_0\mid {\mathring{\mathbf{x}}}_0}{\mathbb{E}}_{{\tilde{\mathbf{x}}}_t\mid {\mathring{\mathbf{x}}}_0,{\tilde{\mathbf{x}}}_0}[{\Vert {\tilde{\mathbf{s}}}_{\mathbf{\varphi }}({\mathbf{x}}_t,\mathring{\mathcal{S}},t)-{\nabla }_{{\tilde{\mathbf{x}}}_t}\log (p({\tilde{\mathbf{x}}}_t\mid {\mathring{\mathbf{x}}}_0,{\tilde{\mathbf{x}}}_0))\Vert }_2^2\left]\right)$$

$$=\underset{\mathbf{\varphi }}{\mathrm{argmin}}{\mathbb{E}}_{t\sim U(0,T)}(\lambda (t){\mathbb{E}}_{\mathring{\mathcal{S}}}{\mathbb{E}}_{{\mathbf{x}}_0}{\mathbb{E}}_{{\tilde{\mathbf{x}}}_t\mid {\mathbf{x}}_0}[{\Vert {\tilde{\mathbf{s}}}_{\mathbf{\varphi }}({\mathbf{x}}_t,\mathring{\mathcal{S}},t)-{\nabla }_{{\tilde{\mathbf{x}}}_t}\log (p({\tilde{\mathbf{x}}}_t\mid {\mathbf{x}}_0))\Vert }_2^2\left]\right),$$

$${\tilde{\mathbf{x}}}_{t+1}=\sqrt{1-{\beta }_t}{\tilde{\mathbf{x}}}_t+\sqrt{{\beta }_t}{\widetilde{\mathbf{ϵ}}}_t,{\widetilde{\mathbf{ϵ}}}_t\sim \mathcal{N}(\tilde{0},\tilde{\mathbf{I}}),t=0,\dots ,T-1,$$

$${\tilde{\mathbf{x}}}_{t-1}=(1-{\beta }_t{)}^{-\frac{1}{2}}({\tilde{\mathbf{x}}}_t+{\beta }_t{\nabla }_{{\tilde{\mathbf{x}}}_t}\log (p({\tilde{\mathbf{x}}}_t\mid {\mathring{\mathbf{x}}}_0)))+\sqrt{{\beta }_t}{\widetilde{\mathbf{ϵ}}}_t,{\widetilde{\mathbf{ϵ}}}_t\sim \mathcal{N}(\tilde{0},\tilde{\mathbf{I}}),t=T,\dots ,1,$$

$$p({\tilde{\mathbf{x}}}_t\mid {\mathbf{x}}_0):=\mathcal{N}(\sqrt{{\underline{\alpha }}_t}{\tilde{\mathbf{x}}}_0,{\underline{\sigma }}_t^2\tilde{\mathbf{I}}),$$

$${\mathbf{s}}_{\mathbf{\varphi }}(\cdot ,\mathring{\mathcal{S}},t)=({\tilde{\mathbf{s}}}_{\mathbf{\varphi }}(\cdot ,\mathring{\mathcal{S}},t{)}^{\prime },{\mathring{0}}_{\mathring{\mathcal{S}}}^{\prime }{)}^{\prime }\in {\mathbb{R}}^{|\mathcal{S}|},$$

$$M_i(\mathring{\mathcal{S}})=\{\begin{array}{ll}1,& {\mathbf{s}}_i\in \mathring{\mathcal{S}},\\ 0,& {\mathbf{s}}_i\in \tilde{\mathcal{S}}.\end{array}$$

$$\ell (\mathbf{\varphi })$$

$$={\mathbb{E}}_{t\sim U(0,T)}(\lambda (t){\mathbb{E}}_{{\mathbf{x}}_t}(⟨{\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t)-{\nabla }_{{\mathbf{x}}_t}\log (p({\mathbf{x}}_t)),{\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t)-{\nabla }_{{\mathbf{x}}_t}\log (p({\mathbf{x}}_t))⟩\left)\right)$$

$$={\mathbb{E}}_{t\sim U(0,T)}(\lambda (t)({\mathbb{E}}_{{\mathbf{x}}_t}({\Vert {\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t)\Vert }^2)-2{\mathbb{E}}_{{\mathbf{x}}_t}(⟨{\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t),{\nabla }_{{\mathbf{x}}_t}\log (p({\mathbf{x}}_t))⟩)$$

$$+{\mathbb{E}}_{{\mathbf{x}}_t}({\Vert {\nabla }_{{\mathbf{x}}_t}\log (p({\mathbf{x}}_t))\Vert }^2)\left)\right),$$

$${\mathbb{E}}_{{\mathbf{x}}_t}\left(⟨{\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t),{\nabla }_{{\mathbf{x}}_t}\log (p({\mathbf{x}}_t))⟩\right)$$

$$=\int p({\mathbf{x}}_t)⟨{\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t),\frac{{\nabla }_{{\mathbf{x}}_t}p({\mathbf{x}}_t)}{p({\mathbf{x}}_t)}⟩\mathrm{d}{\mathbf{x}}_t$$

$$=\int ⟨{\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t),{\nabla }_{{\mathbf{x}}_t}p({\mathbf{x}}_t)⟩\mathrm{d}{\mathbf{x}}_t$$

$$=\int ⟨{\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t),{\nabla }_{{\mathbf{x}}_t}\int p({\mathbf{x}}_0)p({\mathbf{x}}_t\mid {\mathbf{x}}_0)\mathrm{d}{\mathbf{x}}_0⟩\mathrm{d}{\mathbf{x}}_t$$

$$=\int ⟨{\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t),\int p({\mathbf{x}}_0){\nabla }_{{\mathbf{x}}_t}p({\mathbf{x}}_t\mid {\mathbf{x}}_0)\mathrm{d}{\mathbf{x}}_0⟩\mathrm{d}{\mathbf{x}}_t$$

$$=\int ⟨{\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t),\int p({\mathbf{x}}_0)p({\mathbf{x}}_t\mid {\mathbf{x}}_0){\nabla }_{{\mathbf{x}}_t}\log (p({\mathbf{x}}_t\mid {\mathbf{x}}_0))\mathrm{d}{\mathbf{x}}_0⟩\mathrm{d}{\mathbf{x}}_t$$

$$=\iint p({\mathbf{x}}_0)p({\mathbf{x}}_t\mid {\mathbf{x}}_0)⟨{\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t),{\nabla }_{{\mathbf{x}}_t}\log (p({\mathbf{x}}_t\mid {\mathbf{x}}_0))⟩\mathrm{d}{\mathbf{x}}_0\mathrm{d}{\mathbf{x}}_t$$

$$=\iint p({\mathbf{x}}_0)p({\mathbf{x}}_t\mid {\mathbf{x}}_0)⟨{\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t),{\nabla }_{{\mathbf{x}}_t}\log (p({\mathbf{x}}_t\mid {\mathbf{x}}_0))⟩\mathrm{d}{\mathbf{x}}_t\mathrm{d}{\mathbf{x}}_0$$

$$={\mathbb{E}}_{{\mathbf{x}}_0\sim p({\mathbf{x}}_0)}{\mathbb{E}}_{{\mathbf{x}}_t\sim p({\mathbf{x}}_t\mid {\mathbf{x}}_0)}⟨{\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t),{\nabla }_{{\mathbf{x}}_t}\log (p({\mathbf{x}}_t\mid {\mathbf{x}}_0))⟩,$$

$${\mathbb{E}}_{{\mathbf{x}}_t\sim p({\mathbf{x}}_t)}({\Vert {\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t)\Vert }^2)={\mathbb{E}}_{{\mathbf{x}}_0\sim p({\mathbf{x}}_0)}{\mathbb{E}}_{{\mathbf{x}}_t\sim p({\mathbf{x}}_t\mid {\mathbf{x}}_0)}({\Vert {\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t)\Vert }^2).$$

$${\mathbf{\varphi }}^{\ast }$$

$$+{\mathbb{E}}_{{\mathbf{x}}_t}({\Vert {\nabla }_{{\mathbf{x}}_t}\log (p({\mathbf{x}}_t))\Vert }^2)\left)\right))$$

$$=\underset{\mathbf{\varphi }}{\mathrm{argmin}}{\mathbb{E}}_{t\sim U(0,T)}(\lambda (t)({\mathbb{E}}_{{\mathbf{x}}_0\sim p({\mathbf{x}}_0)}{\mathbb{E}}_{{\mathbf{x}}_t\sim p({\mathbf{x}}_t\mid {\mathbf{x}}_0)}({\Vert {\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t)\Vert }^2)$$

$$-2{\mathbb{E}}_{{\mathbf{x}}_0\sim p({\mathbf{x}}_0)}{\mathbb{E}}_{{\mathbf{x}}_t\sim p({\mathbf{x}}_t\mid {\mathbf{x}}_0)}⟨{\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t),{\nabla }_{{\mathbf{x}}_t}\log (p({\mathbf{x}}_t\mid {\mathbf{x}}_0))⟩+{\mathbb{E}}_{{\mathbf{x}}_t}({\Vert {\nabla }_{{\mathbf{x}}_t}\log (p({\mathbf{x}}_t))\Vert }^2)\left)\right)$$

$$=\underset{\mathbf{\varphi }}{\mathrm{argmin}}{\mathbb{E}}_{t\sim U(0,T)}(\lambda (t){\mathbb{E}}_{{\mathbf{x}}_0\sim p({\mathbf{x}}_0)}{\mathbb{E}}_{{\mathbf{x}}_t\sim p({\mathbf{x}}_t\mid {\mathbf{x}}_0)}({\Vert {\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t)\Vert }^2-2⟨{\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t),{\nabla }_{{\mathbf{x}}_t}\log (p({\mathbf{x}}_t\mid {\mathbf{x}}_0))⟩\left)\right)$$

$$=\underset{\mathbf{\varphi }}{\mathrm{argmin}}{\mathbb{E}}_{t\sim U(0,T)}(\lambda (t){\mathbb{E}}_{{\mathbf{x}}_0\sim p({\mathbf{x}}_0)}{\mathbb{E}}_{{\mathbf{x}}_t\sim p({\mathbf{x}}_t\mid {\mathbf{x}}_0)}({\Vert {\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t)\Vert }^2-2⟨{\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t),{\nabla }_{{\mathbf{x}}_t}\log (p({\mathbf{x}}_t\mid {\mathbf{x}}_0))⟩$$

$$+{\Vert {\nabla }_{{\mathbf{x}}_t}\log (p({\mathbf{x}}_t\mid {\mathbf{x}}_0))\Vert }^2\left)\right)$$

$$=\underset{\mathbf{\varphi }}{\mathrm{argmin}}{\mathbb{E}}_{t\sim U(0,T)}(\lambda (t){\mathbb{E}}_{{\mathbf{x}}_0\sim p({\mathbf{x}}_0)}{\mathbb{E}}_{{\mathbf{x}}_t\sim p({\mathbf{x}}_t\mid {\mathbf{x}}_0)}({\Vert {\mathbf{s}}_{\mathbf{\varphi }}({\mathbf{x}}_t,t)-{\nabla }_{{\mathbf{x}}_t}\log (p({\mathbf{x}}_t\mid {\mathbf{x}}_0))\Vert }^2\left)\right).$$

$${\mathbf{x}}_{t+1}$$

Full text

Keywords: approximate simulation, diffusion model, generative model, likelihood-free inference, simulation-based inference, spatial extremes

1 Introduction

Spatial process models are widely used across various disciplines, including environmental science, epidemiology, and social science. However, their practical implementation comes with several challenges, particularly in computation: fitting a spatial process model to data and using it to generate predictions and prediction intervals can be computationally demanding. While numerous approximations exist for Gaussian spatial process models when exact methods are unavailable, options remain more limited for complex non-Gaussian spatial models.

In recent years, neural networks have been increasingly used to tackle computational challenges in the first step of data analysis: model fitting through parameter inference (Louppe2016; Cramner; Dalmasso2020; Dalmasso2024). These techniques have been applied in various ways to the spatial setting. For instance, Walchessen2024 leverage neural networks to efficiently evaluate the likelihood function of a spatial process model, while Lenzi2021 and SD2024; SD2025 employ neural networks as Bayes estimators for parameter estimation. Both approaches rely on simulation-based inference (SBI) paired with machine learning—the practice of utilizing simulations from a statistical model to train a machine learning model for statistical inference. These methods preclude the need for an explicitly known or computationally tractable likelihood function (Cramner; Dalmasso2020). In contrast, less work has focused on the subsequent—yet equally important—spatial inference task: making predictions and quantifying uncertainty at unobserved locations when the predictive distributions are not available in closed form. Existing methods, which typically do not use SBI and machine learning, include Markov chain Monte Carlo (MCMC) (Gelfand) and variational inference (Jordan1999; Banerjee). However, these approaches often involve a trade-off between computational efficiency and statistical accuracy and usually become impractical when dealing with a large number of observations or prediction locations.

Here, we propose a new method using SBI to train a neural generative model that can simulate from predictive distributions in a spatial setting. Specifically, we employ state-of-the-art diffusion models that are used for image generation (Yang2024), because regularly gridded spatial data can effectively be viewed as image data. These models leverage both forward diffusion and backward evolution to generate samples from a target distribution in a computationally efficient manner. In the forward diffusion process, samples from a complex target distribution are progressively transformed into samples from a simpler reference distribution. This forward process, described by a stochastic differential equation (SDE), is designed to ensure the existence of a corresponding reverse SDE that evolves samples from the reference distribution back into samples from the target distribution (Song2021). This reverse SDE is characterized by a combination of user-defined terms from the forward SDE and an unknown score function, which is approximated using a neural network trained via a specifically designed loss function. Once the neural network is trained, a sample from the complex target distribution can be drawn efficiently by simulating from the reference distribution and applying the learned reversed SDE.

There are several approaches to incorporating conditional information into a diffusion process, many of which originate from the inpainting problem, a topic extensively studied in computer vision (e.g., Lugmayr2022; Quan2024; Simkus2025). Generally, in image inpainting, the objective is to produce realistic, yet diverse images from a partial image. These approaches usually make no assumptions about an underlying statistical model and, consequently, are not designed to approximate any true predictive distributions. In contrast, cardoso2025 adopt an approach that first trains a diffusion model to unconditionally sample from a spatial process and then performs conditional simulation from the predictive distribution via training-free guidance. However, this approach has not been applied to non-Gaussian processes and relies on several approximations during conditional simulation, which can lead to slower and potentially less accurate generation. A few inpainting methods, such as twisted diffusion (Wu2023) and controllable generation (Song2021), have asymptotic guarantees for exact conditional simulation. Yet, these methods have not thus far been validated in a statistically rigorous manner and may not scale easily to high-dimensional data.

A recent contribution by Gloeckler2024 provides a masking approach to incorporate conditional information into the score-based diffusion process that is designed to accurately simulate from conditional distributions of parameters and data including posterior distributions. However, their approach was not designed for high-dimensional conditioning, including spatial prediction and uncertainty quantification. A more explicit masking approach for handling missing or censored data is also well-established for estimating model parameters in SBI. Wang2024, Richards2024, and SD2025 include a mask labeling missing or censored data as input to the neural network for parameter inference tasks.

Here, we combine the masking approach of Gloeckler2024 with the more explicit masking approach of Richards2024 for the spatial process model setting. The result is a conditional score-based diffusion model within the SDE framework which can quickly and accurately conditionally simulate from high-dimensional predictive distributions. Specifically, we train a neural network that approximates the time-evolving conditional score function, taking as inputs the observed locations, the observed values, the spatial process model parameters, and the diffusion timestep. Importantly, training the neural network only requires unconditional simulation from the spatial process model, which is feasible for a broad class of models. Once trained, the neural network can be used to evaluate the conditional score function in the conditional reverse process, enabling conditional sampling for any observations—regardless of value, number, and location of these observations—and any model parameters without retraining the neural network. Our approach to conditional simulation, which we refer to as neural conditional simulation (NCS), is said to be amortized with respect to these quantities.

We carry out experiments comparing NCS to classical methods for conditional simulation with two spatial processes: (i) a Gaussian process and (ii) a Brown–Resnick process (BrownResnick1977; Kabluchko2009), which is a non-Gaussian max-stable process commonly used to model spatial extremes (Huser2014; Engelke2015; Davison2015). In the former case, exact conditional simulation is both tractable and computationally efficient; this case study serves to validate NCS and verify that it correctly samples from the true predictive distribution. For the latter case, exact conditional simulation is computationally and analytically intractable in moderate to high dimensions (Dombry2012); the standard approximation involving Markov Chain Monte Carlo (MCMC) has been implemented in the SpatialExtremes package (SpatialExtremes). We show that NCS is both more accurate and computationally efficient than this standard approximation across a range of observed patterns and spatial process model parameters. From these case studies, we conclude that NCS is a viable methodology for sampling from intractable predictive distributions of spatial processes with both high accuracy and computational efficiency.

The paper is organized as follows. Section 2 provides an overview of score-based diffusion models for generative modeling, and describes how they can be adapted for conditional simulation of spatial processes. In Section 3, we evaluate NCS on simulations from a Brown–Resnick spatial process and show that the NCS-generated simulations are representative of the true predictive distributions via a variety of validation metrics. Finally, in Section 4, we discuss the limitations of NCS as well as possible extensions.

2 Methodology

2.1 General Setting

Consider a finite-variance spatial process $h(\cdot)$ dependent on parameter $\bm{\theta}$ and defined on a spatial domain $\mathcal{D}\subset\mathbb{R}^{d}$ (typically $d=2,3)$. Denote by $\mathcal{S}\mathrel{\mathop{:}}=\{\bm{s}_{i}\in\mathcal{D}:i\in[n]\}$ a collection of $n$ locations, where $[n]$ is shorthand notation for the set $\{1,\dots,n\}$. Further, let $\mathring{\mathcal{S}}\subset\mathcal{S}$ be the set of spatial locations at which we observe $h(\cdot)$, and let $\mathcal{\tilde{S}}\mathrel{\mathop{:}}=\mathcal{S}\setminus\mathring{\mathcal{S}}$ be the set of spatial locations at which we wish to predict $h(\cdot)$. The task of conditional simulation can be summarized as follows:

[TABLE]

where the target conditional distribution $p(\tilde{\bm{x}}_{0}\mid\bm{\mathring{x}}_{0};\bm{\theta})$ in (1c), which we also refer to as the predictive distribution, is induced by the stochastic process $h(\cdot)$. Note that there are multiple predictive distributions because $\bm{\theta}$, $\mathring{\mathcal{S}}$, and $\bm{\mathring{x}}_{0}$ can vary. We assume $\bm{\theta}$ is known when simulating from a predictive distribution; in practice, one often replaces $\bm{\theta}$ with an estimate $\bm{\hat{\theta}}$. We assume that, for any $\bm{\theta}$ and $\mathring{\mathcal{S}}$, it is straightforward to unconditionally simulate observations from the spatial process according to (1b) in which $\mathring{\mathcal{S}}=\mathcal{S}$. Note that although we focus on spatial processes, our methodology can be used for any continuous stochastic process defined on a fixed, finite, Euclidean domain.

In Section 2.2, we describe score-based diffusion models for simulating unconditionally from $h(\mathcal{S})$; in Section 2.3, we modify the score-based diffusion model to conditionally simulate from (1c) for any $\mathring{\mathcal{S}}\subset{\mathcal{S}}$. In Section 2.4, we describe the discretizations of the forward and reverse SDEs characterizing the conditional score-based diffusion. In Section 2.5, we describe the neural network architecture we use for approximating the conditional score function of a spatial process. Finally, in Section 2.6, we discuss validation methods. For ease of exposition, we omit the dependence of the spatial process and the induced quantities on the parameter $\bm{\theta}$ in the remainder of Section 2.

2.2 Unconditional Score-Based Diffusion

Here, we first consider the case in which we want to sample $\bm{x}_{0}:=h(\mathcal{S})\sim p(\bm{x}_{0})$, ${\bm{x}_{0}\in\mathbb{R}^{\lvert\mathcal{S}\rvert}}$, where $\lvert\mathcal{S}\rvert=n$ is the number of locations, via an unconditional score-based diffusion model. In score-based diffusion, one first defines an SDE that diffuses $\bm{x}_{0}$ into a sample $\bm{x}_{T}\sim p(\bm{x}_{T})$, where $p(\bm{x}_{T})$ closely approximates a fixed, predetermined reference distribution $p_{\textrm{ref}}(\bm{x})$ selected for ease of sampling (Pedrotti2024). Provided the reverse process exists, one can then approximately sample from $p(\bm{x}_{0})$ by first sampling $\bm{x}_{T}\sim p_{\textrm{ref}}(\bm{x})$ and then evolving the sample through the corresponding reverse SDE (Song2021).

In our approach, we specify the forward process as a sequence of Gaussian perturbations, such that the reference distribution $p_{\textrm{ref}}(\bm{x})$ is Gaussian. One can verify that, as $T\to\infty$, the distribution $p(\bm{x}_{T})$ converges to the reference distribution. The forward SDE we use has the following form, with $t\in[0,T]$ denoting time and $T$ the time horizon:

[TABLE]

where $\bm{B}_{t}$ is an $n$-dimensional Brownian motion, and where the drift and diffusion coefficients $\bm{f}(\cdot,\cdot):\mathbb{R}^{\lvert\mathcal{S}\rvert}\times[0,T]\to\mathbb{R}^{\lvert\mathcal{S}\rvert}$ and $g(\cdot):[0,T]\to\mathbb{R}^{+}$ are user-specified and fixed. Note that time here is artificial and distinct from any temporal dimension in the variables being diffused. The reverse SDE is also a diffusion process given by

[TABLE]

where $\bm{R}_{t}$ is an $n$-dimensional reverse Brownian motion (i.e., a Brownian motion that evolves backwards in time) and $p(\bm{x}_{t})$ is the distribution of the samples $\bm{x}_{t}$ from the perturbed target distribution at time $t$ (Anderson1982). Note that the drift coefficient in the reverse SDE differs from that in the forward SDE as it includes the score function $\nabla_{\bm{x}_{t}}\log(p(\bm{x}_{t}))$ evaluated at the sample $\bm{x}_{t}$, for time $t\in[0,T]$. If the target distribution $p(\bm{x}_{0})$ is non-Gaussian, $p(\bm{x}_{t})$ cannot be adequately approximated by a Gaussian density for $t$ sufficiently smaller than $T$, and hence the score function is generally analytically or computationally intractable.

Indirect simulation from $p(\bm{x}_{0})$ via the reverse process is possible by replacing the score function in (3b) with an accurate surrogate $\bm{s}_{\bm{\phi}}(\bm{x}_{t},t)$—a neural network with trainable parameters $\bm{\phi}$ in this case. One way to ensure the surrogate sufficiently approximates the score function for all $\bm{x}_{t}\in\mathbb{R}^{\lvert\mathcal{S}\rvert}$ and $t\in[0,T]$ is by score matching, which involves solving the following minimization problem (Song2021):

[TABLE]

where $\lambda:[0,T]\to\mathbb{R}^{+}$ is a time-varying weighting function, and where the loss function includes an expectation over time to ensure the surrogate $\bm{s}_{\bm{\phi}}(\bm{x}_{t},\cdot)$ learns the score function for any $t\in(0,T)$ in the reverse SDE. For completeness, we provide in Supplementary Materials S1.1 a proof of the second equality in (4) derived from that in Vincent2011.

Since the forward SDE (2) only involves perturbations of $\bm{x}_{0}$ with Gaussian noise, the transition kernel $p(\bm{x}_{t}\mid\bm{x}_{0})$ is Gaussian provided that the forward drift coefficient $\bm{f}(\bm{x}_{t},t)$ is an affine function of $\bm{x}_{t}$ for all $t\in[0,T]$. Consequently, the score function of the transition kernel $\nabla_{\bm{x}_{t}}\log(p(\bm{x}_{t}\mid\bm{x}_{0}))$ can be easily computed, allowing us to learn the unconditional score $\nabla_{\bm{x}_{t}}\log p(\bm{x}_{t})$ in (4) (Song2021). Further, for affine $f(\cdot,t)$, as $T\to\infty$, $p(\bm{x}_{T})$ approaches a multivariate Gaussian distribution, which we can use as a reference distribution. In this paper, we specifically use a variance-preserving SDE (VPSDE) in which the forward process employs an affine drift coefficient to diffuse $\bm{x}_{0}$ into $\bm{x}_{T}\penalty 10000\ \dot{\sim}\penalty 10000\ \mathcal{N}(\bm{0},\bm{I})$ (where $\dot{\sim}$ denotes “approximately distributed as”), and we therefore set $p_{\textrm{ref}}(\bm{x}):=\mathcal{N}(\bm{0},\bm{I})$ such that $\bm{I}$ is an $n\times n$ identity matrix. See Section 2.4 for more details.

2.3 Conditional Score-Based Diffusion

In this section, we modify the unconditional score-based diffusion approach of Section 2.2 to generate samples from the predictive distribution at $\tilde{\mathcal{S}}$ using observations at $\mathring{\mathcal{S}}$, such that $\mathcal{S}={\tilde{\mathcal{S}}}\cup{\mathring{\mathcal{S}}}$ and $\tilde{\mathcal{S}}\cap\mathring{\mathcal{S}}=\emptyset$. We adapt the approach of Gloeckler2024 to allow conditional simulation from predictive distributions involving arbitrary partitionings of $\mathcal{S}$ in a spatial setting. We modify the forward and reverse SDEs in (2) and (3), respectively, so that only $\tilde{\bm{x}}_{t}$, the unobserved component of $\bm{x}_{t}=(\tilde{\bm{x}}^{\prime}_{t},\bm{\mathring{x}}^{\prime}_{t})^{\prime}$, is diffused. Specifically, we use the following conditional forward SDE:

[TABLE]

and the following conditional reverse-time SDE:

[TABLE]

where $\tilde{\bm{f}}(\cdot,\cdot)$ is the component of $\bm{f}(\cdot,\cdot)$ corresponding to $\tilde{\mathcal{S}}$, and where $\tilde{\bm{B}}_{t}$ and $\tilde{\bm{R}}_{t}$ are forward- and reverse-time Brownian motions defined on $\mathbb{R}^{\lvert\tilde{\mathcal{S}}\rvert}$, respectively. Note that there are generally different conditional forward and reverse SDEs for different initial conditions, $\mathring{\mathcal{S}}$ and $\bm{\mathring{x}}_{0}$, in (5a) and (6a).

Consider a distribution over the observed locations, $p(\mathring{\mathcal{S}})$, that has positive mass over the set of all possible partitions of $\mathcal{S}$. By the average optimality property for non-negative loss functions (Zammit2024), the modified score matching algorithm is given by:

[TABLE]

where $\tilde{\bm{s}}_{\bm{\phi}}(\cdot,\cdot,\cdot)$ is the output of the surrogate conditional score function corresponding to $\tilde{\mathcal{S}}$, $\bm{x}_{t}=(\bm{\mathring{x}}_{0}^{\prime},\tilde{\bm{x}}_{t}^{\prime})^{\prime}$, and the second line follows from the same argument as that of (4).

Crucially, (7a) reveals that one does not need to sample from or have knowledge of the predictive distributions $p(\tilde{\bm{x}}_{0}\mid\bm{\mathring{x}}_{0})$ to train the conditional score-matching surrogate—one only needs to simulate from the user-specified $p(\mathring{\mathcal{S}})$ and the unconditional distribution of the spatial process model before running the forward evolution (5). Thus, our diffusion model can be used with a wide class of process models that are straightforward to sample from unconditionally, but difficult to sample from conditionally. Note that the output dimension of the surrogate $\tilde{\bm{s}}_{\bm{\phi}}(\cdot,\cdot,\cdot)$ changes with $\lvert\tilde{\mathcal{S}}\rvert$; Section 2.5 describes how to accommodate this in practice. As in Section 2.2, we use a VPSDE, as detailed in the next section.

2.4 Discretization

A popular SDE used for score-based diffusion models is the VPSDE, which, in the conditional setting, sets $\tilde{\bm{f}}(\tilde{\bm{x}}_{t},t)=\frac{-\beta(t)}{2}\tilde{\bm{x}}_{t}$ and $g(t)=\sqrt{\beta(t)}$ in (5b), where $\beta(t):[0,T]\to\mathbb{R}^{+}$ is selected in such a manner that the process has bounded variance as $T\to\infty$ (Song2021). To evolve $\tilde{\bm{x}}_{t}$ with the conditional forward and reverse SDEs, we discretize time. The discretized versions of (5b) and (6b) are given by

[TABLE]

and

[TABLE]

where $\beta_{t}=\beta(t)$ and $\tilde{\bm{\epsilon}}_{t}$ is defined on $\tilde{\mathcal{S}}$. See Supplementary Materials S1.2 and Appendix B in Song2021 for derivation details. The corresponding transition kernels in the discretized conditional forward SDE have the form

[TABLE]

where $\underline{\sigma}_{t}=1-\underline{\alpha}_{t}$, $\underline{\alpha}_{t}=\prod_{s=1}^{t-1}\alpha_{s}$ and $\alpha_{s}=1-\beta_{s}$. To sample from the target distribution $p(\tilde{\bm{x}}_{0}\mid\bm{\mathring{x}}_{0})$ using (9), we replace the true score function $\nabla_{\tilde{\bm{x}}_{t}}\log(p(\tilde{\bm{x}}_{t}\mid\bm{\mathring{x}}_{0}))$ with the approximation $\tilde{\bm{s}}_{\bm{\phi}}(\bm{x}_{t},\mathring{\mathcal{S}},t)$. Note that in (7a), $t$ now needs to be equipped with a probability mass function on $\{1,\dots,T\}$ instead of a continuous one. The VPSDE framework for score-based diffusion in (8) and (9) is sometimes also referred to as a denoising diffusion probabilistic model (DDPM); Chapter 4 in Chan2025 clarifies the connection between these two models.

2.5 Architecture for Conditional Score Neural Network

A key consideration in our methodology is the architecture of the neural network surrogate for the conditional score function $\tilde{\bm{s}}_{\bm{\phi}}(\cdot,\cdot,\cdot)$ which must be flexible enough to accurately approximate highly non-linear score functions. Here, we restrict $\mathcal{S}$ to points on a regular grid, which enables the use of convolutional neural network (CNN) architectures (Lecun1998). In what follows, we regard each realization of the spatial process on $\mathcal{S}$ as an image, and $\tilde{\mathcal{S}}$ corresponds to pixels within the image that are treated as unobserved. To cater to score functions of variable dimension, we train a neural network whose input and output at time $t\in[0,T]$ has dimension $\mathbb{R}^{\lvert\mathcal{S}\rvert}$ and which outputs [math] for all observed locations. That is, our neural network,

[TABLE]

yields the appropriate subvector $\tilde{\bm{s}}_{\bm{\phi}}(\cdot,\mathring{\mathcal{S}},t)$ for use in (7a) and (9) as part of its output. For brevity, we assume (11) places the values at the correct locations to avoid cumbersome notation.

For $\bm{s}_{\bm{\phi}}(\cdot,\cdot,\cdot)$, we use a specific type of CNN, referred to as U-Net, which consists of two main parts: a contracting path (or encoder) which reduces the image dimension and yields a large number of feature maps through multiple layers, and an expanding path (or decoder), that uses several layers of transposed convolutions, skip connections, and upsampling stages to restore the original spatial resolution; see Ronneberger2015 for more details. There are many variants of U-Nets. In our implementation, we employ the noise conditional score network++ model (NCSN++) of Song2021 with minor modifications to reduce the U-Net size and thus, decrease training and evaluation time. Specifically, we reduced the number of blocks in both the contracting and expanding paths of the U-Net from $4$ to $2$ and the number of Gaussian Fourier projections determining the sinusoidal time embeddings from $128$ to $64$.

Recall from (7a) that the conditional score depends on the observed locations, and therefore $\mathring{\mathcal{S}}$ needs to be an input to the score matching neural network. Since CNNs can only handle inputs of fixed dimension, we instead input a mask $\bm{M}(\mathring{\mathcal{S}})\in\{0,1\}^{\lvert S\rvert}$ where

[TABLE]

Specifically, we input $\bm{x}_{t}\in\mathbb{R}^{\lvert\mathcal{S}\rvert}$ and ${\bm{M}}(\mathring{\mathcal{S}})\in\{0,1\}^{\lvert\mathcal{S}\rvert}$ as two input channels, enabling the U-Net to extract separate spatial information from both inputs.

2.6 Validation

Validating whether the neural conditional simulations come from the true predictive distributions is difficult due to the high-dimensional setting, the assumed intractability of the predictive distributions, as well as the variable conditioning set. To address these issues, we propose a sample-based validation method derived from Linhart2023, which compares the NCS-approximated unconditional distribution (described below) to the true unconditional distribution. This approach has multiple benefits: First, this method is tractable because simulation via NCS and from the true unconditional distribution is computationally efficient. Second, the NCS-approximated unconditional distribution does not depend on a single conditioning set, and thus, we can sidestep the issue of verifying the predictive distribution for each and every possible conditioning set. Finally, this validation technique enables the comparison of multivariate metrics, such as the spatial minimum and maximum of each unconditional simulation, which will only follow the correct distribution if the corresponding high-dimensional characteristics of the true unconditional distributions are accurately captured. However, as with any validation technique, passing these checks is a necessary but not sufficient condition for correct conditional simulation.

In addition to these multivariate metrics, we also use low-dimensional metrics. For the low-dimensional metrics, we created a validation data set involving replicates generated via NCS and either the true predictive distribution or an MCMC approximation of the true predictive distribution. Since these low-dimensional metrics are straightforward summaries, here we only discuss how to generate the NCS-approximated unconditional distribution in depth, and defer the technical details for the low-dimensional validation metrics to Supplementary Materials S2.3.1.

NCS-Approximated Unconditional Distribution

Consider the setting in which parameters are fixed and the conditioning sets are drawn from a fixed distribution $p(\mathring{\mathcal{S}})$. We simulate $m$ replicates $\{\bm{x}_{0,i}\}_{i\in[m]}$ from the true unconditional distribution and $m$ corresponding conditioning sets $\{\mathring{\mathcal{S}}_{i}\}_{i\in[m]}$ and masks $\{\bm{M}(\mathring{\mathcal{S}}_{i})\}_{i\in[m]}$ according to $p(\mathring{\mathcal{S}})$. For each of the $m$ replicates $\bm{x}_{0,i}$ and masks $\bm{M}(\mathring{\mathcal{S}}_{0,i})$, we obtain the corresponding partially observed spatial fields $\bm{\mathring{x}}_{0,i}$, generate a single NCS simulation $\tilde{\bm{x}}_{0,i}$ conditional on $\bm{\mathring{x}}_{0,i}$, and concatenate the NCS simulation with the partially observed field to produce an “unconditional” simulation $(\bm{\mathring{x}}_{0,i}^{\prime},\tilde{\bm{x}}_{0,i}^{\prime})^{\prime}$. If NCS accurately samples from $p(\tilde{\bm{x}}_{0,i}\mid\bm{\mathring{x}}_{0,i})$, then $({\bm{\mathring{x}}}_{0,i}^{\prime},\tilde{\bm{x}}_{0,i}^{\prime})^{\prime}$ is a sample from $p(\bm{x}_{0})$, and can therefore be validated against other samples drawn from $p(\bm{x}_{0})$. See Algorithms S3 and S4 in the Supplementary Materials for pseudocode of this procedure.

3 Simulation Experiments

To demonstrate the accuracy of NCS, we use it to conditionally simulate from two spatial processes: a Gaussian process and a Brown–Resnick process. The first Gaussian case study verifies that our method is able to generate samples from true (known) predictive distributions; the second Brown–Resnick case study serves to demonstrate the utility of our method in settings where fast and efficient unconditional simulation is possible, but from which the predictive distributions are difficult or intractable to sample. For brevity, here we only present the results for the Brown–Resnick case study; see Supplementary Materials S3 for the Gaussian case study.

3.1 Experimental Setup

We let the spatial locations in $\mathcal{S}$ be a 32 $\times$ 32 gridding on $\mathcal{D}=[-10,10]\times[-10,10]$. The conditioning sets $\mathring{\mathcal{S}}$ and corresponding masks $\bm{M}(\mathring{\mathcal{S}})$ were constructed either by sampling at random a fixed number of locations from the $32^{2}$ grid points (Section 3.2), or through $32^{2}$ independent Bernoulli draws with probability of success (corresponding to whether a pixel is observed or not) equal to $\rho$ (Section 3.3). For the diffusion processes in (8) and (9), we set $T=1000$, $\beta_{0}=0.0001$, $\beta_{T}=0.02$ and we let $\beta_{t}=\beta(t)$ increase linearly with $t$ (specifically, we let $\beta(t)=0.0001+0.0199t/T$). We used a weighting function $\lambda(t)=1-\underline{\alpha}_{t}$ for $\underline{\alpha}_{t}=\prod_{s=1}^{t}(1-\beta_{s})$ in the conditional loss function (7a); see Song2021 and Supplementary Materials S2.2 and S4.2 for details. Training and validation data were generated “on the fly” (Chan2018).

The Brown–Resnick process that we use in this paper is a two-parameter max-stable process with $\bm{\theta}\equiv(\lambda,\nu)^{\prime}$, with range parameter $\lambda\in\mathbb{R}^{+}$ and smoothness parameter $\nu\in(0,2]$. It admits a likelihood function that is extremely costly to evaluate, and predictive distributions that are computationally challenging to sample from (Padoan2010; Davison; Dombry2012; Huser2019). See Supplementary Materials S2.1 for greater detail on the Brown–Resnick process.

For the purpose of assessing NCS, we consider the Gibbs sampler (Dombry2012) implemented in the SpatialExtremes R package (SpatialExtremes) as an alternative conditional simulation method. We found that full conditional simulation (FCS) via the SpatialExtremes package is only computationally feasible when the conditioning set contains seven or fewer observations. For conditioning sets with more than seven observed locations, we resort to local conditional simulation (LCS) by only using the the seven observed locations closest to a prediction location for conditional simulation. In contrast to FCS and NCS, LCS targets a low-dimensional approximation of a predictive distribution and will only be used here for univariate and bivariate approximations of the predictive distributions. The Gibbs sampler is computationally intensive, and requires a computing wall time that depends on several factors, such as the number of observed locations and model parameters. See Supplementary Materials S4.1 for more details.

In contrast, as we show below in Section 3.2, NCS allows for consistent computational efficiency, regardless of the number of observed locations or process-model parameters. We implemented NCS and LCS/FCS on an Nvidia Tesla T4 GPU and an AMD EPYC 7V12 processor (with 64 cores, no multithreading) respectively. For the Brown–Resnick process, we log-transformed the simulated spatial fields from the unit Fréchet scale to the standard Gumbel scale before training the U-Nets for estimating the score function.

3.2 Results: Small Number of Observed Locations

In this section, we compare FCS and NCS for the case when the number of observed locations $|\mathring{\mathcal{S}}|\in\{1,\dots,7\}$ (where FCS is computationally feasible) and where the range parameter $\lambda\in\{1,\dots,5\}$ (weak to strong spatial dependence) varies. For each range value, we trained and implemented separate U-Nets that are each amortized for conditioning sets of $1$ to $7$ locations. See Supplementary Materials S2.2 for details. To assess NCS, we examine the extremal correlation function $\chi(h)=2-\zeta(h)$, where the spatial extremal coefficient function $\zeta(h):\mathbb{R}^{+}\to[1,2]$ quantifies the extremal dependence between pairs of locations a fixed distance $h$ apart (deHaan1985; Smith1990). We compute the empirical extremal correlation functions for the true unconditional distributions, the FCS-approximated unconditional distributions, and the NCS-approximated unconditional distributions using the approach described in Section 2.6 (Linhart2023). In addition to the extremal correlation functions, we also empirically estimate the unconditional densities of the minimum, maximum, and the absolute sum of the process over $\mathcal{S}$. See Supplementary Materials S4 for training and validation details as well as additional results.

In Figure 1, we show the true, NCS, and FCS extremal correlation functions for the case in which only one to seven locations are observed, when the range parameter is $\lambda=1$ (left panel) and $\lambda=5$ (right panel). When $\lambda=1$, the true and NCS extremal correlation functions match exactly, but the FCS extremal correlation function is reflective of a spatial dependence that is too long-range. As shown in Figure S5, the FCS simulations are negatively biased at unobserved locations far from the conditioning sites when the range parameter is small. In contrast, when $\lambda=5$, the true and FCS extremal correlation functions match exactly, whereas the NCS extremal correlation function is slightly inaccurate in the case of seven observed locations. In this case, the NCS simulations appear more diffuse than the FCS simulations. However, this discrepancy is not as stark as the negative bias in the FCS simulations when the range parameter is small.

The empirical distributions of the spatial minimum and maximum computed from the true, FCS-approximated, and NCS-approximated unconditional simulations are shown in the second and third rows of Figure 1. We see that the distributions of the minima and maxima from NCS match the true ones (except for some slight inaccuracy in the $\lambda=5$ case with seven observations), while those from FCS are clearly biased for the $\lambda=1$ case. In the fourth row of Figure 1, we show the distribution of the absolute sum of the process over $\mathcal{S}$. Again, we see that the true and NCS distributions of this quantity match almost exactly. As with the distributions of the spatial minimum and maximum, FCS does not capture the absolute summation distributions well for $\lambda=1$.

To highlight the computational inefficiency and ultimately infeasibility of FCS when compared to NCS, we record the average time (in seconds) to generate one conditional simulation for varying number of observations, and divide it by the time needed to conditionally simulate with $|\mathring{\mathcal{S}}|=1$. As shown through these time ratios in Figure 2, the average time to generate an FCS simulation when seven locations are observed is five times greater than the time needed when one location is observed. Since the Gibbs sampler underlying FCS involves evaluating multivariate Gaussian cumulative distribution functions, we expect FCS to scale exponentially with the number of observed locations $|\mathring{\mathcal{S}}|$ (Hinrich2014; Balakrishnan2014; SpatialExtremes) and use this to show the likely scaling of FCS beyond seven observations—the maximum number of locations for which our implementation of the Gibbs sampler in the SpatialExtremes R package remained stable. We conclude that FCS is difficult to implement using available software, and its computational complexity quickly increases with the number of observed locations.

On the other hand, the time needed to generate an NCS simulation does not change with the number of observed locations. The average time to generate a single FCS and NCS simulation from seven observed locations was $28$ and $27$ seconds, respectively, with our setup. We currently evaluate the U-Net at every time step in the discretized conditional reverse SDE, but there are methods designed to reduce the number of steps needed to simulate from diffusion models (Bond-Taylor2022). Once implemented, these may further reduce the time needed for a single NCS simulation.

3.3 Results: Amortizing with Respect to Observed Proportion and Process Parameters

In this section, we evaluate NCS’s performance on the Brown–Resnick process with larger numbers of observed locations, and assess its adaptability to varying parameter values. Specifically, we consider a scenario where $\bm{\theta}$ is fixed but where $\rho$, the proportion of observed values, is allowed to vary, and another where $\rho$ is fixed and where the range parameter is allowed to vary. We trained the U-Net in the first scenario using simulations from a Brown–Resnick process with a fixed parameter vector $\bm{\theta}=(\lambda,\nu)^{\prime}=(3,1.5)$ and with $\rho$ sampled from a uniform distribution on $[0.01,0.5]$. This U-Net, which we refer to as the “proportion U-Net”, is amortized with respect to the observed proportion of grid locations since it does not require re-training for different proportions $\rho\in[0.01,0.5]$. We trained the second U-Net using simulations from a Brown–Resnick process with fixed $\nu=1.5$ and $\rho=0.05$, but with $\lambda$ allowed to vary uniformly on $[1,5]$. This U-Net, which we refer to as the “parameter U-Net”, is amortized with respect to the range parameter since it does not require re-training for different range values $\lambda\in[1,5]$. While in principle it is possible to train a single U-Net that amortizes with respect to both observed proportion and spatial-process parameters, we anticipate that larger U-Nets and computing capacity would be required for such capability.

When an element of $\bm{\theta}$, $\theta_{1}$ say, is allowed to vary, the score function estimate must also be allowed to vary with $\theta_{1}$. To accommodate this, in the parameter U-Net, we include $\theta_{1}$ as an extra input, i.e., we re-define the approximate score function to be $\bm{s}_{\bm{\phi}}(\bm{x}_{t},\mathring{\mathcal{S}},\theta_{1},t)=\big(\tilde{\bm{s}}_{\bm{\phi}}(\bm{x}_{t},\mathring{\mathcal{S}},\theta_{1},t)^{\prime},\bm{0}_{\mathring{\mathcal{S}}}^{\prime}\big)^{\prime}$, and modify (7a) accordingly to include an expectation with respect to a distribution over $\theta_{1}$ (here, a uniform distribution on $[1,5]$). Through experimentation, we found that information about the parameter $\theta_{1}$ needs to be included at each step in the contracting and expanding paths of the U-Net to get good results. We implemented this by replacing the second input channel to the U-Net, $\bm{M}(\mathring{\mathcal{S}})$, with $\theta_{1}\bm{M}(\mathring{\mathcal{S}})$.

In this section, we validate the resulting NCS simulations with various metrics including empirical conditional univariate and marginal densities as well as empirical unconditional densities of the minimum, maximum, and the absolute sum of the process over $\mathcal{S}$ via the approach described in Section 2.6. See Supplementary Materials S4.2 and S4.3 for training and validation details.

In Figure 3 we show observations, underlying true fields, and conditional simulations obtained using NCS. In the left panel, we show results from fixing $\rho=0.05$, varying the range parameter $\lambda$ between 1 and 5, and implementing NCS using the parameter U-Net. In the right panel we show results from fixing $\lambda=3$, varying the proportion parameter $\rho$ between $0.01$ and $0.5$, and implementing NCS using the proportion U-Net. We see that the NCS conditional simulations exhibit spatial patterns that are similar to those of the corresponding true (partially observed) fields, and that amortization with respect to either the range or the proportion of observed locations appears to be effective.

Figure 4 shows the LCS (purple) and NCS (orange) conditional univariate and bivariate densities for select locations in $\mathcal{S}$. We do not expect the NCS and LCS densities to match exactly since LCS is an approximation that utilizes less information than NCS. Yet, we observe that the two predictive densities at these select locations are very similar. In particular, the densities correctly capture the true values shown by the dashed red lines as either the range $\lambda$ or the proportion of observed locations $\rho$ vary.

Extremal correlation plots, as well as empirical distributions of the minimum, maximum, and absolute sum, are shown in Figure 5. The true and NCS extremal correlation functions match almost exactly, even as the range parameter $\lambda$ varies. The NCS-approximated unconditional distributions of the minimum, maximum, and absolute sum also match almost perfectly with the underlying true ones excepting the case when the proportion of observed locations is $\rho=0.01$—corresponding to approximately $10$ observed locations. In Section 3.2, the metrics for the NCS-approximated unconditional distribution with only seven observed locations similarly deviated from the truth.

3.4 Summary of Results

In this case study, we found that NCS simulations are at least comparable to and often outperform the standard approximations (FCS and LCS) in generating simulations representative of the intractable true predictive distributions of a Brown–Resnick process. Importantly, NCS has fewer limitations: FCS becomes computationally intractable when the number of observed spatial locations exceeds seven, and LCS is inherently restricted in the degree of high-dimensional patterns it can capture. In addition, the high-dimensional patterns for the FCS-approximated unconditional distributions differ from those of the true unconditional distribution more significantly than the equivalent for NCS in many cases except when the range is large and roughly one percent of locations are observed. It is unclear whether this dip in performance is inherent to NCS when applied to Brown–Resnick processes or whether improved training in this case may address the decreased performance.

4 Discussion and Conclusions

In this paper, we introduced neural conditional simulation (NCS) as a method to sample from predictive distributions of a spatial process using a conditional score-based diffusion model within an SDE framework. For each partially observed field and spatial process parameters, there is an individual conditional forward SDE which diffuses the values of the unobserved locations into white noise and a corresponding conditional reverse SDE which evolves white noise to simulations from the true predictive distribution. While the conditional forward SDEs are known, the corresponding conditional reverse SDEs are unknown due to an intractable conditional score function. To sample from the reverse SDEs, we first train a neural network that depends on the values at the observed locations and the locations themselves to approximate the score function. We also show how the score-function neural estimator can be amortized with respect to the proportion of observed locations, and the parameters of the underlying spatial process. Importantly, training the neural network to approximate the conditional score function only requires simulations from the unconditional distribution. NCS thus implicitly learns the true predictive distributions of spatial processes for which unconditional simulation is fast and efficient yet conditional simulation is not. With NCS, fast and accurate spatial prediction and uncertainty quantification using complex, non-Gaussian spatial processes and partially observed data becomes possible.

To demonstrate the accuracy of NCS, we evaluated how representative NCS-generated simulations are of the true predictive distributions for Gaussian processes (Supplementary Materials S3.3) and Brown–Resnick processes (Section 3), which have tractable and intractable predictive distributions, respectively. Through the use of diagnostics that depend on the joint properties of the predictive distribution, such as the minimum and maximum, we gave empirical evidence that the NCS-approximated conditional simulations are accurate. Finally, we showed that NCS remains computationally efficient regardless of the number of observed locations. We also demonstrated that NCS can be more accurate and comes with fewer limitations than MCMC in practical settings. NCS therefore appears to be ideally suited for efficiently simulating from intractable predictive distributions for a wide class of complex spatial processes.

In this paper, we have presented NCS in the context of spatial process models. However, the methodology presented in Section 2.3 is sufficiently flexible to also apply to spatial hierarchical models. Specifically, in the case where one has a data model and a latent process model, (5) and (6b) can be modified slightly so that the target distribution is that of the latent process given the data. Furthermore, we envision that NCS can enable accurate prediction and uncertainty quantification for a broad class of statistical and physical models beyond spatial statistics for which conditional simulation currently necessitates approximations or is computationally infeasible.

Data Availability Statement

The code used to produce these results is publicly available at https://github.com/jmwalchessen/ncs. The simulation experiments in this paper can be reproduced with the given code.

Acknowledgments

Part of this work was done when J.W. visited A.Z.-M. at University of Wollongong and R.H. at the King Abdullah University of Science and Technology (KAUST). We would like to thank Jordan Richards and Brian Trippe for helpful discussions about masking approaches for partially observed spatial data and conditional diffusion, respectively. We would like to acknowledge Microsoft for providing Azure computing resources for this work. J.W., A.Z.-M. and R.H. were supported by KAUST grant ORFS-2023-OFP-55020. J.W. and M.K. were supported in part by National Science Foundation grant DMS-2053804, NOAA grant NA21OAR4310258 and a grant from C3.ai Digital Transformation Institute.

Declaration of Interest

No competing interests are declared.

S1 Supplementary Materials

This section provides detailed justification for the loss function (4) in Section 2.2 and the discretized forward (8) and reverse SDEs (9) in Section 2.4. All the derivations presented here are for the unconditional case; they are applicable to the conditional case with minor adjustments.

S1 Methodology Details

S1.1 Loss Function for Unconditional Score-Based Diffusion

Here, we demonstrate how the seemingly intractable loss function in the first line of (4) can be rewritten as the tractable expression in the second line via a derivation from Vincent2011. First, expand the norm difference in the intractable loss function via matrix multiplication:

[TABLE]

where $\left\langle\cdot,\cdot\right\rangle$ is the inner product operator. Second, rewrite the second term in the expansion as follows,

[TABLE]

where (S2a) is a rearrangement of the gradient $\nabla_{\bm{x}_{t}}\log(p(\bm{x}_{t}\mid\bm{x}_{0}))$ and (S2b) follows from Fubini’s Theorem. Next, the first term in the expansion of $\ell(\bm{\phi})$ can be reworked via the Law of Total Expectation,

[TABLE]

Finally, substitute the reworked first and second terms into the expanded loss function in (S1) and minimize $\ell(\bm{\phi})$ with respect to the parameters $\bm{\phi}$ according to

[TABLE]

S1.2 Discretization of the Forward and Reverse SDEs

Here, we provide the details on the time-discretized forward and reverse SDEs for the unconditional case. Note that the derivation is applicable to the conditional case, as well, with minor adjustments. For simplicity, we assume $T$ is sufficiently large such that ${{\Delta{t}/T}\approx 0}$ for $\Delta t=1$.

To emphasize the connection between denoising diffusion probabilistic models (DDPM) and the VPSDE framework for score-based diffusion, we start by writing down our derivations from the forward and reverse processes of a DDPM (Ho2020; Chan2025). The discretized forward process of a DDPM is

[TABLE]

and therefore,

[TABLE]

where $\bm{f}(\bm{x}_{t},t)=-{\beta_{t}\bm{x}_{t}}/2$, $g(t)=\sqrt{\beta_{t}}$, and the first approximation follows due to a MacLaurin series expansion of $h(\beta_{t})=\sqrt{1-\beta_{t}}$ for $\beta_{t}$ chosen such that $\beta_{t}\approx 0$ for $t\in[0,T]$. The discretized reverse process of a DDPM is

[TABLE]

and therefore,

[TABLE]

where $\bm{f}(\bm{x}_{t},t)=-{\beta_{t}\bm{x}_{t}}/2$, $g(t)=\sqrt{\beta_{t}}$, and the first approximation follows due to a MacLaurin series expansion of $h(\beta_{t})=(1-\beta_{t})^{-\frac{1}{2}}$ for $\beta_{t}$ chosen such that $\beta_{t}\approx 0$ for $t\in[0,T]$.

S2 Case Study Methodology

S2.1 Brief Description of Spatial Processes

Gaussian Process

A Gaussian process is a stochastic process whose finite-dimensional distributions are Gaussian. Its predictive distributions are analytically available, and therefore we can use them to validate NCS. Results from applying NCS to Gaussian processes are shown in Section S3.3.

Brown–Resnick Process

Here, we give more details on the Brown–Resnick process employed in Section 3.1. Max-stable processes, such as Brown–Resnick processes, can be expressed as the pointwise maxima of an infinite collection of random spatial processes (Davison; Huser2019; Davison2019). More specifically, a max-stable process has the following construction (deHaan1985)

[TABLE]

in which $\{\eta_{i}\}_{i=1}^{\infty}$ is a Poisson point process on $(0,\infty)$ with intensity function $\operatorname{d}\!{\Lambda}(\eta)=\eta^{-2}\operatorname{d}\!{\eta}$ and $\{W_{i}(\bm{s})\}$ are independent copies of a non-negative stochastic process with mean one (Kabluchko2009). Different types of max-stable processes arise from the definition of $W_{i}(\bm{s})$. For a Brown–Resnick model, the process $W_{i}(\bm{s})$ is log-normal and has the form:

[TABLE]

where $\{\epsilon_{i}(\bm{s})\}_{i\in\mathbb{N}^{+}}$ are realizations of an intrinsically stationary Gaussian process, conventionally with $\epsilon_{i}(\bm{s})=0$ almost surely, and with semivariogram $\gamma(\bm{h})$, $i\in\mathbb{N}^{+}$, where $\mathbb{N}^{+}$ denotes the positive integers. A common choice, which we also adopt in this work, is to set $\gamma(\bm{h})=(\left.\left\lVert\bm{h}\right\rVert\middle/\lambda\right.)^{\nu}$, with $\bm{h}$ the spatial separation between two locations in spatial domain $\mathcal{D}$ , $\lambda\in\mathbb{R}^{+}$ a range parameter, and $\nu\in(0,2]$ a smoothness parameter (Castruccio).

S2.2 Training Details Common to the Gaussian Process and Brown–Resnick Process Case Studies

Weighting Function

In our loss function (7a), we use the same weighting function $\lambda(\cdot)$ as in Song2021, which is proportional to the inverse of the expectation of the conditional score function. This is equivalent to setting $\lambda(t)=\underline{\sigma}_{t}^{2}$ for a variance preserving SDE (VPSDE), since, from (10),

[TABLE]

Spatial Process Parameters

In the main text, we consider two differently amortized U-Nets, one amortized with respect to a spatial process parameter referred to as a parameter U-Net and the other amortized with respect to the proportion of missing observations referred to as a proportion U-Net. For the parameter U-Net, we consider the parameter space $\Theta=[1,5]\times\{1.5\}$, where $\nu=1.5$ is fixed for both processes described in Section S2.1, and where the length scale $\ell$ and range $\lambda$ vary from one to five for the Gaussian process and Brown–Resnick process, respectively. We define the training parameter space $\tilde{\Theta}$—the space from which we draw unconditional simulations during training—to extend beyond $\Theta=[1,5]\times\{1.5\}$ to avoid suboptimal behavior of the U-Net at the boundary. On the other hand, when training the proportion U-Net, we do not vary the spatial process parameters. Specifically, we fix the length scale $\ell=3$ and variance $\tau^{2}=1.5$ for the Gaussian process, and the range $\lambda=3$ and smoothness $\nu=1.5$ for the Brown–Resnick process. Therefore, for the parameter U-Net, $\Theta=\tilde{\Theta}=\{3\}\times\{1.5\}$.

Simulating Training and Validation Data

We employ data generation “on the fly.” That is, we simulate a dataset, use it to train and evaluate the diffusion model for a fixed number of epochs, and repeat the process with a new simulated dataset each time. This approach avoids overfitting and minimizes data storage (Chan2018).

We simulated training and validation data for the parameter and proportion U-Nets according to the following procedure: For each data draw, we first simulate $p$ parameters from the training parameter space $\tilde{\Theta}$, and for each of those $p$ parameters, we simulate $s$ full spatial field realizations via the unconditional simulator. Note that for the proportion U-Net, these $p$ parameter vectors are identical. Finally, we simulate $m$ masks and $m$ discrete timesteps uniformly from $\{0,1,\dots,T\}\subset\mathbb{Z}$ per each full spatial field realization. Each mask is simulated using the Bernoulli distribution with $\rho$—the probability of observing a single location—for all $n=32^{2}$ locations, and our convention is that [math] indicates missingness. For the parameter U-Net, $\rho=0.05$ is kept fixed. For the proportion U-Net, we sample $\rho$ uniformly from $[0.01,0.525]$ during training to avoid boundary issues during evaluation, where we instead restrict the space of proportions to $[0.01,0.5]$. Algorithm S1 illustrates this data generation process. Once finished with simulation, we have training data of the form:

[TABLE]

where $\bm{\mathring{x}}_{0}=(x_{0,j}:M_{j}(\mathring{\mathcal{S}})=1)^{\prime},\quad\tilde{\bm{x}}_{t}=\sqrt{\underline{\alpha}_{t}}\cdot\mathring{\bm{x}}_{0}+\underline{\sigma}_{t}\tilde{\bm{\epsilon}},\quad\tilde{\bm{\epsilon}}\sim\mathcal{N}(\tilde{\bm{0}},\tilde{\bm{I}})$, and where $M_{j}(\mathring{\mathcal{S}})$ is defined in (12).

S2.3 Validation Details Common to Both the Gaussian Process and Brown–Resnick Process Case Studies

S2.3.1 Generation of Validation Datasets

In this section, we give details on how we generate test data for the simulation experiments conducted in Section 3 and Supplementary Materials S3 and S4.

Conditional Validation Dataset (Proportion U-Net)

For each of the five observed proportions $\rho\in\{0.01,0.05,0.1,0.25,0.5\}$, we generated a conditioning set $\mathring{\mathcal{S}}$ and mask $\bm{M}(\mathring{\mathcal{S}})$ using the Bernoulli distribution with probability $\rho$, and a reference spatial field $\bm{x}_{0}$ from its unconditional distribution with either variance $\tau^{2}=1.5$ and length scale $\ell=3$ (Gaussian process), or smoothness $\nu=1.5$ and range $\lambda=3$ (Brown–Resnick process). For each of the resulting five partially observed spatial fields $\bm{\mathring{x}}_{0}=(x_{0,j}:M_{j}(\mathring{\mathcal{S}})=1)^{\prime}$, we generated $m=4000$ replicates of the unobserved spatial field $\{\tilde{\bm{x}}_{0,i}\}_{i\in[m]}$ with NCS and either the true predictive distribution (Gaussian process) or an MCMC approximation (Brown–Resnick process). We use this dataset in Section 3.3 to validate NCS when the percentage of observed locations varies.

Conditional Validation Dataset (Parameter U-Net)

In the Gaussian process case, we unconditionally simulated a reference spatial field $\bm{x}_{0}$ for each parameter length scale $\ell\in\{1,\ldots,5\}$ and variance $\tau^{2}=1.5$. In the Brown–Resnick case, we unconditionally simulated a reference spatial field $\bm{x}_{0}$ for each range parameter $\lambda\in\{1,\ldots,5\}$ and smoothness $\nu=1.5$. Then, for each reference spatial field, we generated a conditioning set $\mathring{\mathcal{S}}$ and corresponding mask $\bm{M}(\mathring{\mathcal{S}})$ using the Bernoulli distribution with probability $\rho=0.05$ of observing each location. In each case, using each of the five partially observed fields $\bm{\mathring{x}}_{0}=(x_{0,j}:M_{j}(\mathring{\mathcal{S}})=1)^{\prime}$, we generated $m=4000$ unobserved spatial field replicates $\{\tilde{\bm{x}}_{0,i}\}_{i\in[m]}$ with NCS and either the true predictive distribution (Gaussian process) or an MCMC approximation (Brown–Resnick process). We use this dataset in Section 3.3 to validate NCS when the spatial process parameter varies.

S2.3.2 Diagnostic and Validation Metrics

In this section, we describe diagnostics and metrics used to assess the quality of conditional simulation in Sections 3.2, 3.3, S3.2 and S4.4.

Conditional Univariate and Bivariate Densities

Conditional univariate and bivariate densities are estimated using the conditional validation datasets and kernel density estimation with a Gaussian kernel and bandwidth selected according to Scott’s Rule (Scott1992) on the conditional validation datasets. Specifically, in our results (e.g., Figure 4), we show a $5\times 3$ grid of plots in which the first column contains the partially observed spatial fields and the second and third columns contain the conditional univariate and bivariate kernel density estimates respectively for the selected unobserved locations.

Conditional Mean Fields

The empirical conditional mean fields are computed by averaging the $m=4000$ replicates at each of the unobserved locations in the conditional validation dataset. In our results (e.g., Figures S3 and S6), we display a $5\times 3$ grid of plots that contains the partially observed spatial fields in the first column, the true empirical (Gaussian) or an MCMC-approximated (Brown–Resnick) conditional mean field in the second column, and the empirical NCS conditional mean field in the third column.

S3 Gaussian Process: Tractable Case Study

We applied NCS to a Gaussian process with zero mean and exponential covariance function parametrized by the length scale $\ell>0$ and the variance $\tau^{2}>0$. Hence, the model is

[TABLE]

S3.1 Training Details

Parameter U-Net

For the parameter U-Net in the Gaussian process case, we fixed the variance to $\tau^{2}=1.5$ and the proportion of observed locations to $\rho=0.05$, but let the length scale $\ell\in[1.0,5.0]$ vary. To simulate training data, we follow the process described in Section S2.2. Specifically, for training, we simulated $p=1000$ parameters samples from $\tilde{\Theta}=[0.5,5.5]\times\{1.5\}$, $s$ = 1 full spatial field realizations per parameter sample, and $m=100$ masks and timesteps per full spatial field realization.

Proportion U-Net

For the proportion U-Net in the Gaussian process case, we fixed the variance to $\tau^{2}=1.5$, the length scale to $\ell=3$, while we allowed the proportion of observed locations $\rho\in[0.01,0.5]$ to vary. We follow the process in Section S2.2 for simulating training data. Specifically, for training, we sampled $r=50$ observed proportions $\rho$ uniformly from $[0.01,0.525]$, $p=1$ parameters uniformly from $\Theta=\tilde{\Theta}=\{3\}\times\{1.5\}$ for each $\rho$, $s=25$ full spatial field realizations for each $\rho$, and used $m=100$ masks and timesteps for each full spatial field realization.

S3.2 Validation Metrics Particular to the Gaussian Process Case Study

Conditional Correlation Metrics

To assess how well the NCS simulations capture the spatial correlations induced by the true predictive distributions in the Gaussian case study, we produce empirical conditional correlation heatmaps (Figure S4). For a fixed unobserved spatial location in a partially observed spatial field from the conditional validation dataset, we compute the Pearson correlation coefficient (PCC) between the unobserved spatial location and all other unobserved spatial locations using both NCS simulations and conditional simulations from the true predictive distribution. PCC measures the linear correlation between the process evaluated at two unobserved spatial locations $\bm{s}_{1},\bm{s}_{2}$ given the partially observed spatial field. If the PCC is $-1$ or $1$, there is perfect linear correlation given the partially observed spatial field, while a PCC value of [math] indicates no linear correlation given the partially observed field. Specifically, we display a $5\times 2$ grid of plots for both the parameter and proportion U-Nets, in which the first column contains the true empirical correlation heat maps and the second column contains the heatmaps derived from NCS.

S3.3 Results

Visualizations

In Figure S1, we show simulations from the true predictive distributions and NCS. These simulations display similar spatial patterns and seamlessly incorporate the partially observed field as intended. Additionally, amortization appears to be highly effective: as the length scale increases in the left panel, the spatial correlations increase in the NCS simulations, mirroring the true conditional simulations. Further, no matter the number of observed spatial locations, NCS seamlessly integrates the simulated unobserved values with the partially observed field as shown in the right panel. From these results, we conclude that NCS simulations are reasonable as one cannot discriminate between the truth and NCS by the naked eye.

Conditional Marginal and Bivariate Densities

The true predictive distributions of a Gaussian process are fully determined by its first and second moments. As such, if all the true and NCS conditional univariate and bivariate densities are equivalent up to some small approximation error due to finite sample estimation, then the NCS simulations are representative of the true predictive distributions. In Figure S2, the true and NCS conditional univariate and bivariate densities appear equivalent up to approximation error. While Figure S2 only shows densities for a select number of spatial locations and partially observed spatial fields, we compared conditional univariate and bivariate densities for several other spatial locations and observed similar behavior. Thus, we conclude that the first and second moments of the NCS and true predictive distributions are approximately equivalent, and that the NCS simulations are representative of the true predictive distributions.

Conditional Mean and Correlation Metrics

Figure S3 and Figure S4 contain the true and NCS conditional mean fields and correlation heatmaps, respectively. These serve as additional diagnostics to assess whether—and, if so, how—the first and second moments of the NCS and true predictive distributions differ. In both figures, there are little to no visual differences between the true and NCS empirical conditional mean fields and correlation heatmaps. For both the NCS and true predictive distributions in Figure S4, spatial (conditional) correlations generally increase with the length scale, unless there are observed locations (marked in white) near the spatial location of interest (marked in yellow). This behavior is best illustrated in the fourth row of the left panel; the conditional correlations sharply decrease above the spatial location of interest (in yellow) due to the presence of an observation two pixels above (in white). Yet, the conditional correlations are unaffected below the spatial location of interest due to the absence of observed data in that area. Such behavior is expected for the true predictive distributions; it is yet a striking indication how effective NCS is at correctly conditioning on the given partially observed spatial fields.

S4 Brown–Resnick: Intractable Case Study

S4.1 Details on Standard Approximation for the Brown–Resnick Process

Further Description

For both FCS and LCS, we utilize the function condrmaxstab in the SpatialExtremes package to approximate the predictive distributions (SpatialExtremes). This function applies a three-step procedure in which the first step involves identifying feasible partitions of the $k$ observed locations based on whether they share common extremal functions (Dombry2012). The extremal functions $\varphi_{\bm{s}_{i}}^{+}$ of a max-stable realization $\bm{z}$ are functions $\eta_{i}W_{i}(\bm{s})$ from the stochastic representation of deHaan1985 in (S9) over the set of spatial locations $\mathcal{S}$ for which the realization $\bm{z}$ and the extremal function $\varphi_{\bm{s}_{i}}^{+}$ evaluated at location $\bm{s}_{i}$ are equal. Extremal functions are not necessarily unique for all spatial locations $\bm{s}_{i}\in\mathcal{S}$ because an extremal function $\varphi_{\bm{s}_{i}}^{+}$ evaluated at another location $\bm{s}_{j}$ may equal the realization $\bm{z}$ at location $\bm{s}_{j}$ as well. As such, this first step requires an approximation via a MCMC algorithm, specifically a Gibbs sampler, to achieve computational tractability because the number of all possible partitions for $k$ locations is the Bell number $B_{k}$ which grows more than exponentially with respect to $k$ (Castruccio). The second and third steps involve simulating extremal and sub-extremal functions respectively from tractable, closed-form distributions and as such, are less computationally intensive. See Dombry2012 for full mathematical details and Huser2019 for another application and assessment of the above Gibbs sampler.

As shown in Figure S5 below, FCS suffers from negative bias for unobserved locations far away from the spatial domain center when the range parameter is small. This is probably due to the insufficient number of subextremal functions drawn by the condrmaxstab function, which affects mostly locations far away from the observed locations when the variogram increases to infinity.

Implementation

To implement FCS and LCS, we use the condrmaxstab function with the given set of observed locations (FCS) or the nearest seven neighbors to the select unobserved location or locations (LCS), the nugget set to $0.00001$, the burn-in period set to 1000, and the thinning factor set to $100$. Here, a nugget of $0.00001$ adds an insignificant amount of noise to the Brown–Resnick unconditional simulations used during the third step, but helps address computational problems which occur in the absence of a nugget.

We use univariate and bivariate LCS to produce conditional univariate and bivariate densities and conditional mean fields. In the univariate case, for LCS, we condition on the seven observations closest to the location of interest. In the bivariate case, we first identify the three nearest neighbors for each of the two locations of interest, and then randomly sample from the remaining unique seven nearest neighbors for each of the two select locations until a conditioning set of size seven is obtained. This process ensures there are neighbors for each of the two locations of interest in the conditioning set for bivariate LCS.

S4.2 Training Details

The major difference between the Gaussian and Brown–Resnick case studies is the log-transformation of the training data simulated from the Brown–Resnick process. As such, NCS produces conditional simulations on the (log-transformed) Gumbel scale, and all results are on the Gumbel scale.

Parameter U-Net

The details for simulating training are the same as for the Parameter U-Net in the Gaussian process case study (Section S3.1) except that the number of proportion values was set to $r=200$ for the training data.

Proportion U-Net

The details for simulating training data are the same as for the Proportion U-Net in the Gaussian process case study.

Small Conditioning Set U-Net

We train U-Nets for the small conditioning sets via the following data generation process: We first fix the parameter $\bm{\theta}=(\lambda,\nu)^{\prime}$ and select a minimum $o_{L}=1$ and maximum $o_{U}=10$ number of observations for the conditioning sets according to the experimental settings. For each data draw, we iterate through each number of observed locations $o_{j}$ starting from the minimum $o_{L}$ and ending with the maximum $o_{U}$. For each observed number $o_{j}$, we simulate $s$ full spatial field realizations via the unconditional simulator with fixed parameter $\bm{\theta}$. For each of these $s$ full spatial fields, we construct $m$ masks and simulate $m$ discrete timesteps uniformly from $\{0,1,\dots,T\}$. Each mask is constructed from a conditioning set $\mathring{\mathcal{S}}$ containing $o_{j}$ locations uniformly sampled without replacement from $\mathcal{S}$. Algorithm S2 illustrates this data process, and the resulting data have the form shown in (S12).

Specifically, we train one U-Net for each specific range parameter ${\lambda\in\{1,2,3,4,5\}}$. However, we fix the smoothness to $\nu=1.5$ and the minimum and maximum number of observations to one and ten respectively for all U-Nets. For validation, the number of observed locations does not exceed seven because FCS as implemented in the SpatialExtremes R package is fragile beyond seven. To simulate training and validation data, we follow the process described in Section S2.2. For simulating training data, the minimum and maximum number of observations were set to $o_{L}=1$ and $o_{U}=10$, $s=32$ full spatial field realizations were generated per number of observations between $o_{L}$ and $o_{U}$, and $m=100$ conditioning sets, masks and timesteps were simulated for each full spatial field realization.

S4.3 Validation Details

Conditional Validation Dataset (Small Conditioning Set)

For each range value $\lambda\in\{1,2,3,4,5\}$ and each conditioning set size $m$ ranging from one to seven, we unconditionally simulated a reference spatial field $\bm{x}_{0}$ from a Brown–Resnick process with the given range $\lambda$ and fixed smoothness $\nu=1.5$. For each reference spatial field, we generated a conditioning set $\mathring{\mathcal{S}}$ and corresponding mask $\bm{M}(\mathring{\mathcal{S}})$ by uniformly sampling from $\mathcal{S}$ without replacement, such that $\lvert\mathring{\mathcal{S}}\rvert=m$. Using each of the $35$ partially observed fields $\bm{\mathring{x}}_{0}=(x_{0,j}:M_{j}(\mathring{\mathcal{S}})=1)^{\prime}$, we generated $m=4000$ full spatial field replicates $\{(\bm{\mathring{x}}_{0}^{\prime},\tilde{\bm{x}}_{0,i}^{\prime})^{\prime}\}_{i\in[m]}$ with NCS and FCS (an MCMC approximation described in Section 3.1). We use this dataset in Section 3.2 to validate NCS when the number of observed locations does not exceed seven.

Unconditional Validation Dataset (Small Conditioning Set U-Net)

For the validation method described in Section 2.6, we need to generate unconditional samples from the conditional ones. As discussed in Section 2.6, we do this by generating unconditionally at observed locations and then conditionally simulating at unobserved locations; the resulting concatenated sample is then one from the unconditional distribution. Algorithm S3 summarizes how we do this, using $m=4000$ NCS-approximated or FCS-approximated unconditional simulations using varying range ${\lambda\in\{1,2,3,4,5\}}$, fixed smoothness $\nu=1.5$, and conditioning sets with varying sizes ranging from one to seven.

Unconditional Validation Dataset (Parameter and Proportion U-Net)

We apply the process described in Section 2.6 and the previous paragraph but with varying parameters and observed proportions as shown in Algorithm S4. We generate $m=4000$ NCS-approximated unconditional simulations using either varying range $\lambda\in\{1,2,3,4,5\}$, fixed smoothness $\nu=1.5$, and fixed observed proportion $\rho=.05$ for the parameter U-Net, or varying observed proportion $\rho\in\{0.01,0.05,0.1,0.25,0.5\}$ with fixed range $\lambda=3$ and smoothness $\nu=1.5$ for the proportion U-Net.

Multivariate Metrics

The metrics in Section S2.3.2 such as the conditional univariate and bivariate densities, conditional mean fields, and conditional correlation metrics only compare low-dimensional properties of NCS to those of the true predictive distributions or standard approximations. Here, we discuss metrics that implicitly validate the high-dimensional patterns of NCS-approximated predictive distributions via NCS-approximated unconditional simulations. These metrics, referred to as multivariate metrics, are evaluated using the aforementioned unconditional validation datasets. For the Brown–Resnick case study, where these multivariate metrics are especially relevant, we considered the extremal correlation function and the spatial minimum, maximum, and absolute summation in any given full realization. Distributions of these quantities are shown in our results (e.g., Figure 5) for different parameters and missingness cases.

S4.4 Results

Visualizations (Small Conditioning Set)

In Figure S5, we show results for both FCS (middle) and NCS (right) simulations for the case when three locations are observed. Clearly, FCS suffers from negative bias for unobserved locations far away from the spatial domain center when the range parameter is small. This is due to the way MCMC is implemented for this approximation method. Due to this negative bias, the NCS simulations appear more representative overall, although NCS suffers from slightly increased graininess as the range parameter increases.

Conditional Mean Fields (Small Conditioning Set)

In Figure S6, we show results for various spatial process parameter and proportion settings. We see no discernible differences between the univariate LCS and NCS empirical conditional mean fields.