We introduce and analyze a nonlocal generalization of Whittle--Matérn Gaussian fields in which the smoothness parameter varies in space through the fractional order, $s=s(x)\in[\underline{s}\,,\bar{s}]\subset(0,1)$. The model is defined via an integral-form operator whose kernel is constructed from the modified Bessel function of the second kind and whose local singularity is governed by the symmetric exponent $β(x,y)=(s(x)+s(y))/2$. This variable-order nonlocal formulation departs from the classical constant-order pseudodifferential setting and raises new analytic and numerical challenges. We develop a novel variational framework adapted to the kernel, prove existence and uniqueness of weak solutions on truncated bounded domains, and derive Sobolev regularity of the Gaussian (spectral) solution controlled by the minimal local order: realizations lie in $H^r(G)$ for every $r<2\underline{s}-\tfrac{d}{2}$ (here $H^r(G)$ denotes the Sobolev space on the bounded domain $G$), hence in $L_2(G)$ when $\underline s>d/4$. We also present a finite-element sampling method for the integral model, derive error estimates, and provide numerical experiments in one dimension that illustrate the impact of spatially varying smoothness on samples covariances. Computational aspects and directions for scalable implementations are discussed.
1. Introduction. Gaussian random fields are fundamental models in statistics and machine learning. Gaussian random fields can be specified through their covariance function, and the most commonly used covariance function is the Matérn covariance function [23],
where Γ denotes the gamma function, ∥ • ∥ the Euclidean norm on R d and K ν is the modified Bessel function of the second kind of order ν. The parameters σ, ν, κ > 0 control the variance, smoothness, and correlation range of the random field, respectively.
One drawback of this model is that the parameters are constant over space, which means that it cannot capture non-stationarities, such as a spatially varying practical correlation range. Because non-stationary covariance structures often are observed in applications in statistics, there has been considerable interest over the past 30 years to construct flexible non-stationary extensions of the Matérn covariance function. Some popular approaches are space-deformations [30], process convolutions [18], and more explicit constructions which directly formulate valid non-stationary covariance functions [25]. One of the most popular approaches is the stochastic partial differential equation (SPDE) approach [22], which is based on the fact that a stationary solution u : R d ×Ω → R (Ω is a sample space) to the fractional-order SPDE
has covariance (1.1) [36], with range parameter κ, smoothness ν = 2s -d/2, and variance
Here ∆ is the Laplacian, µ > 0 is a scaling factor and W is a spatial Gaussian white noise in R d . Motivated by the SPDE formulation, [22] proposed extending the Matérn model to accommodate non-stationary behavior and more general spatial domains. Their work sparked a vibrant line of research focused on spatial models constructed via SPDEs; see, for example, [4,7,14,19,9,5,8]. In particular, the authors derived a connection between Gaussian fields and Gaussian Markov random fields on a bounded domain D ⊂ R d by using an approximate weak solution to (1.2), posed on a bounded domain D ⊂ R d , where the operator was equipped with Neumann boundary conditions. Besides reducing computational cost compared to the covariance-based approximations, [22] also used this idea to extend the Matérn fields to non-stationary models by allowing the parameters κ and µ to be spatially varying functions. This idea was later used to construct the so-called generalized Whittle-Matérn fields, which are formulated by considering the SPDE
where the operator is given as L = κ 2 -∇ • (H∇), for bounded and non-negative functions κ, µ and a (sufficiently nice) matrix-valued function H [6,10]. The term Whittle-Matérn fields is used to denote, in a broader sense, the class of solutions to (1.4), encompassing not only the stationary fields defined on R d , but also solutions formulated on bounded domains with boundary conditions and on general manifold domains.
A drawback of the finite element (FE) approximation proposed by [22] is that it is computable only if 2s ∈ N. This restriction has later been removed by combining a FE approximation with a rational approximation of either the fractional operator L -s [6] or of the corresponding covariance operator L -2s [9]. Either option results in a computationally efficient approximation which is applicable for s > d/4. In contrast to the approach in [22], where Neumann boundary conditions are employed, the authors in [6] consider homogeneous Dirichlet boundary conditions. However, the choice of boundary conditions typically has little practical impact, as the error in the covariance of the field decays rapidly away from the boundary [26].
Although the SPDE approach can be used to construct non-stationary Gaussian processes, it is limited to models with constant smoothness. Similarly, most previously proposed methods for constructing non-stationary covariance functions, such as the deformation approach, do not allow for spatially varying smoothness. This is a significant limitation, as many physical phenomena naturally exhibit variability in smoothness across space. For example if precipitation or temperature is modeled on a global scale, it is likely that these fields are smoother over oceans compared to over land. Motivated by this, we aim in this work to extend the SPDE approach to allow a spatially varying fractional order s(x), which in turn controls the smoothness of the field. To accomplish this, we consider an integral-based formulation of the fractional operator instead of a spectral definition of it as used, for example, in [22] and [6].
The main contributions of this paper are summarized below.
• We formulate a variable-order Whittle-Matérn integral operator with a heterogeneous Bessel-type kernel on a bounded domain G. • We develop a rigorous variational framework, including completeness of the adapted energy space, denoted V κ,s , and a well-posed weak formulation. • We establish spectral and regularity results: compact embedding V κ,s → L 2 (G), compact resolv
This content is AI-processed based on open access ArXiv data.