Whittle-Matérn Fields with Variable Smoothness

1. Problem Statement

Gaussian random fields (GRFs) are fundamental in statistics and machine learning, typically modeled using the Matérn covariance function. A primary limitation of the standard Matérn model is that its parameters (variance, smoothness, and correlation range) are constant across space. This assumption of stationarity fails to capture real-world phenomena where smoothness varies spatially (e.g., precipitation fields being smoother over oceans than over land).

While previous extensions using Stochastic Partial Differential Equations (SPDEs) allowed for spatially varying correlation ranges and variances, they generally retained a constant fractional order (and thus constant smoothness). Existing methods for non-stationary smoothness often rely on space deformations or process convolutions, which lack the rigorous spectral and regularity guarantees of the SPDE approach.

The Core Challenge: The authors aim to generalize the Whittle–Matérn field to allow for a spatially varying fractional order $s = s(x) \in [s, \bar{s}] \subset (0, 1)$. This introduces significant analytical difficulties because the resulting operator is nonlocal with a heterogeneous, singular kernel, departing from the classical constant-order pseudodifferential setting.

2. Methodology

A. Integral-Form Operator Formulation

Instead of using a spectral definition (which is difficult for variable orders on bounded domains), the authors define the model via a nonlocal integral operator.

• Kernel Construction: The operator is constructed using the modified Bessel function of the second kind, $K_\nu$.
• Variable Exponent: The local singularity is governed by a symmetric exponent $\beta(x, y) = \frac{s(x) + s(y)}{2}$.
• The Operator: For a function $u$ on a bounded domain $G$, the operator $\tilde{A}_s$ is defined as:
  $$\tilde{A}_s u(x) = \kappa^{2s(x)} u(x) + \int_{\mathbb{R}^d} \tilde{w}_{-2s}(x, y) \frac{u(x) - u(y)}{\|x-y\|^{d+2\beta(x,y)}} dy$$
  where $\tilde{w}_{-2s}$ involves the Bessel function $K_{\frac{d}{2}+\beta(x,y)}$.
• Domain Truncation: To make the problem computable, the exterior domain is truncated to a finite region $D_c^t$, and the problem is posed on $G = D \cup D_c^t$ with zero boundary conditions outside $D_c^t$.

B. Variational Framework

The authors develop a rigorous functional analytic framework:

1. Energy Space ($V_{\kappa, s}$): They define a Hilbert space adapted to the heterogeneous kernel, equipped with an energy norm combining an $L^2$ mass term and a nonlocal seminorm involving the kernel $\gamma(x,y)$.
2. Well-Posedness: They prove that $V_{\kappa, s}$ is complete, that smooth compactly supported functions are dense in it, and that the energy spaces for different $\kappa$ are equivalent.
3. Weak Formulation: The SPDE $\tilde{A}_s u = \mu^{-1} W$ (where $W$ is Gaussian white noise) is reformulated as a variational problem: Find $u \in V_{\kappa, s}$ such that $A_s(u, v) = \mu^{-1} W(v)$ for all test functions $v$.

C. Spectral and Regularity Analysis

• Operator Properties: Using the Lax-Milgram theorem and form representation theorems, they establish the existence of a unique, positive, self-adjoint operator $A$ with a compact resolvent.
• Spectral Bounds: They derive two-sided spectral bounds for the eigenvalues $\lambda_j$ of $A$, showing they behave asymptotically like $j^{2s/d}$ (where $s$ is the minimal local order).
• Sobolev Regularity: A key theoretical result is that the spectral solution $u$ belongs to the Sobolev space $H^r(G)$ almost surely for any $r < 2s - d/2$. This implies $u \in L^2(G)$ provided $s > d/4$.

D. Numerical Implementation

• Finite Element Method (FEM): A Galerkin FEM discretization is proposed. Due to the nonlocal nature, the stiffness matrix is dense ($O(N^2)$).
• Quadrature: Specialized quadrature rules (Duffy-type transformations) are developed to handle the singular near-field integrals arising when mesh elements share vertices or are identical.
• Sampling: The solution is sampled by solving the linear system $Au = b$, where $b$ is a realization of the white noise projected onto the finite element basis.

3. Key Contributions

1. Novel Integral Model: Introduction of a variable-order Whittle–Matérn field defined via a nonlocal integral operator with a Bessel-type kernel, explicitly handling spatially varying smoothness $s(x)$.
2. Rigorous Variational Theory: Proof of the completeness of the adapted energy space $V_{\kappa, s}$ and the well-posedness of the weak formulation.
3. Regularity Results: Establishment of Sobolev regularity for the solution controlled by the minimal local order $s = \min s(x)$, proving $u \in H^r(G)$ for $r < 2s - d/2$.
4. Uniqueness: Proof of the uniqueness of the Gaussian generalized solution within the class of generalized fields.
5. Numerical Algorithm: Development of a finite-element sampler with error estimates and specialized quadrature for singular kernels, demonstrated in 1D.

4. Key Results

• Theoretical:
  • The energy space $V_{\kappa, s}$ is a Hilbert space, and the embedding $V_{\kappa, s} \hookrightarrow L^2(G)$ is compact.
  • The eigenvalues $\lambda_j$ of the associated operator satisfy $C_1 j^{2s/d} \leq \lambda_j \leq C_2 j^{2s/d}$.
  • The solution is almost surely in $L^2(G)$ if $s > d/4$.
• Numerical (1D Experiments):
  • Convergence: The observed convergence rates for the FEM approximation align with the theoretical bound $O(h^{2s - d/2})$ (specifically $2s - 0.5$ in 1D).
  • Covariance Behavior: Numerical experiments with step, Gaussian bump, and oscillatory ramp functions for $s(x)$ demonstrate that the sample covariance adapts locally to the smoothness.
    • In regions with lower $s(x)$, the field is rougher (less smooth).
    • In regions with higher $s(x)$, the field is smoother.
    • The sample covariance $C_s(x, y)$ at a point $x$ is bounded by the Matérn covariances corresponding to the minimum and maximum smoothness values in the neighborhood, effectively capturing the non-stationarity.

5. Significance

• Bridging Theory and Application: This work provides a mathematically rigorous foundation for modeling non-stationary Gaussian fields with variable smoothness, a capability previously lacking in the SPDE framework.
• Physical Relevance: The model is directly applicable to physical phenomena where smoothness varies (e.g., geophysics, climate modeling, material science), offering a more realistic alternative to stationary models.
• Computational Pathway: While the current implementation is limited to 1D due to the $O(N^2)$ cost of the dense nonlocal matrix, the paper lays the groundwork for scalable implementations (e.g., using hierarchical matrices or wavelet methods) discussed as future work.
• Statistical Inference: The framework enables the potential for learning spatial heterogeneity directly from data, allowing for the calibration of $s(x)$ in inverse problems.

In summary, the paper successfully extends the popular Whittle–Matérn model to the variable-order regime, overcoming significant analytical hurdles to provide a robust, nonlocal framework for spatially varying smoothness in Gaussian random fields.