Fast simulation of Gaussian random fields with flexible correlation models in Euclidean spaces

Abstract

The efficient simulation of Gaussian random fields with flexible correlation structures is fundamental in spatial statistics, machine learning, and uncertainty quantification. In this work, we revisit the \emph{spectral turning-bands} (STB) method as a versatile and scalable framework for simulating isotropic Gaussian random fields with a broad range of covariance models. Beyond the classical Matérn family, we show that the STB approach can be extended to two recent and flexible correlation classes that generalize the Matérn model: the Bummer-Tricomi model, which allows for polynomially decaying correlations and long-range dependence, and the Gauss-Hypergeometric model, which admits compactly supported correlations, including the Generalized Wendland family as a special case. We derive exact stochastic representations for both families: a Beta-prime mixture formulation for the Kummer-Tricomi model and complementary Beta- and Gasper-mixture representations for the Gauss-Hypergeometric model. These formulations enable exact, numerically stable, and computationally efficient simulation with linear complexity in the number of spectral components. Numerical experiments confirm the accuracy and computational stability of the proposed algorithms across a wide range of parameter configurations, demonstrating their practical viability for large-scale spatial modeling. As an application, we use the proposed STB simulators to perform parametric bootstrap for standard error estimation and model selection under weighted pairwise composite likelihood in the analysis of a large climate dataset.

Figures (7)

• Figure 1: Median simulation time (in seconds) versus the number of spatial locations $n$, shown on a log--log scale, for SPDE--INLA sampling, Vecchia approximations with $m=5,10,20,30$, and the proposed STB algorithm. The left panel corresponds to $\nu=0.5$, and the right panel to $\nu=1$ . Black crosses connected by a dashed line indicate the one-time cost of mesh construction and precision-matrix assembly for INLA; this cost is not included in the per-realization simulation times.
• Figure 2: Simulated Gaussian RFs on $1{,}000{,}000$ points uniformly distributed on $[0,1]^2$ using the generalized Wendland model $\mathcal{GW}_{\nu,7,0.2}\equiv\widetilde{\mathcal{GH}}_{\nu,7,1/2,0.2}$, for $\nu=0,1,2$ from left to right.
• Figure 3: Validity (yellow and green) and simulatability (green) regions in the $(\nu,\mu)$ plane for $d=2$. The left panel corresponds to the $\mathcal{GW}$ model ($l=0.5$), while the right panel refers to the $\mathcal{H}$ model ($l=d/2+\nu$). The solid black line denotes the validity boundary, whereas the dashed black line marks the simulatability limit of the Beta--mixture construction.
• Figure 4: Empirical semivariograms from $1{,}000$ realizations of a Gaussian random field under different parameter configurations of the $\mathcal{GW}_{\nu,\mu,a}$ model (see Table \ref{['tab:tabla1']}), generated using the Beta--mixture algorithm. Gray lines show the empirical semivariograms for individual realizations, red lines their empirical mean over the $1{,}000$ realizations, and blue lines the corresponding theoretical semivariograms. Green lines show empirical semivariograms computed from realizations obtained via the Cholesky method.
• Figure 5: Empirical semivariograms from $1{,}000$ realizations of a Gaussian random field under different parameter configurations of the $\mathcal{GW}_{\nu,\mu,a}$ model (see Table \ref{['tab:tabla2']}), generated using the Gasper--mixture algorithm. Gray lines show the empirical semivariograms for individual realizations, red lines their empirical mean over the $1{,}000$ realizations, and blue lines the corresponding theoretical semivariograms. Green lines show empirical semivariograms computed from realizations obtained via the Cholesky method.