Dirichlet-Neumann Averaging: The DNA of Efficient Gaussian Process Simulation

TL;DR

This paper introduces Dirichlet-Neumann averaging (DNA), a novel, efficient method for generating Gaussian process and Gaussian random field realizations on high-resolution grids, with negligible error and applications to SPDE boundary conditions.

Contribution

The paper presents a new DNA sampling methodology for GPs and GRFs that is faster and maintains accuracy, and links this approach to SPDE boundary condition selection.

Findings

DNA sampling is computationally efficient with negligible covariance error.

Explicit error estimates are provided for Matérn covariances.

A link between DNA methodology and SPDE boundary conditions is established.

Abstract

Gaussian processes (GPs) and Gaussian random fields (GRFs) are essential for modelling spatially varying stochastic phenomena. Yet, the efficient generation of corresponding realisations on high-resolution grids remains challenging, particularly when a large number of realisations are required. This paper presents two novel contributions. First, we propose a new methodology based on Dirichlet-Neumann averaging (DNA) to generate GPs and GRFs with isotropic covariance on regularly spaced grids. The combination of discrete cosine and sine transforms in the DNA sampling approach allows for rapid evaluations without the need for modification or padding of the desired covariance function. While this introduces an error in the covariance, our numerical experiments show that this error is negligible for most relevant applications, representing a trade-off between efficiency and precision. We provide explicit error estimates for Mat\'ern covariances. The second contribution links our new methodology to the stochastic partial differential equation (SPDE) approach for sampling GRFs. We demonstrate that the concepts developed in our methodology can also guide the selection of boundary conditions in the SPDE framework. We prove that averaging specific GRFs sampled via the SPDE approach yields genuinely isotropic realisations without domain extension, with the error bounds established in the first part remaining valid.