📊 statistics

Non-stationary Spatial Modeling Using Fractional SPDEs

This paper introduces a flexible Gaussian random field model based on fractional stochastic partial differential equations with spatially varying parameters, demonstrating through simulations and case studies that appropriate penalization prevents overfitting while the relative benefits of capturing non-stationarity versus fractional smoothness depend on the specific application.

Original authors: Elling Svee, Geir-Arne Fuglstad

Published 2026-06-26

📖 5 min read🧠 Deep dive

CC BY 4.0

Original authors: Elling Svee, Geir-Arne Fuglstad

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you are trying to draw a map of the weather or the ocean based on measurements taken at specific points. Usually, scientists use a standard "rulebook" (a mathematical model) to guess what the conditions are in the spaces between those points. This rulebook assumes the world is somewhat uniform: that the wind blows the same way everywhere and that the temperature changes smoothly in the same pattern across the entire map.

However, the real world is messy. In one part of the ocean, the water might be smooth and calm, while in another, it's choppy and turbulent. In one area, rain might fall in a gentle, predictable drizzle, while nearby, it might be a chaotic storm. The old rulebooks struggle with this because they force the whole map to follow one single pattern.

This paper introduces a new, smarter "rulebook" called Fractional SPDEs that can handle this messiness. Here is how it works, broken down into simple concepts:

1. The Flexible Rubber Sheet (Non-Stationarity)

Think of the old models as a rigid, flat sheet of plastic. If you stretch it, it stretches the same way everywhere. The new model is like a smart, stretchy rubber sheet.

The Problem: In the ocean, the "range" (how far one point influences its neighbor) might be short in a narrow fjord but long in the open sea. The old sheet can't change its stretchiness locally.
The Solution: This new model allows the "stretchiness" (anisotropy) and the "roughness" of the map to change from place to place. It can be tight and detailed in one corner and loose and smooth in another, just like the real environment.

2. The "Fractional" Smoothness (The Texture)

Imagine you are feeling a surface.

Smooth: Like a sheet of glass (easy to predict).
Rough: Like sandpaper (hard to predict).
The Old Way: The old models could only choose between "Glass" or "Sandpaper." They couldn't handle "Fine Grit" or "Medium Grit."
The New Way: This paper introduces Fractional Smoothness. It allows the model to choose any texture in between. It can say, "This part of the ocean is like fine sandpaper, but that part is like velvet." This is crucial because nature often has textures that aren't perfectly smooth or perfectly rough.

3. The "Overfitting" Trap (The Pen)

When you give a model too much freedom (like letting it change the stretchiness and texture everywhere), it can get confused. It might start memorizing the noise (random errors) in your data instead of learning the real pattern. This is called overfitting.

The Analogy: Imagine a student who memorizes the exact answers to a practice test. They get 100% on the practice, but fail the real exam because they didn't learn the concepts, just the specific numbers.
The Fix: The authors built a "penalty system" (a mathematical leash). If the model tries to get too crazy or too specific with its changes, the penalty pulls it back to a more reasonable, general shape. They tested this and found that even with many data points, the model didn't get confused; it stayed focused on the big picture.

4. How They Tested It (The Simulation)

The authors ran computer simulations to see how much data you need to make this new model work.

The Finding: You need a decent amount of data to teach the model these new tricks. They found that you need at least 500 measurement points in a single map to reliably figure out both the changing textures and the changing smoothness. If you have fewer than that, the model might get lost.
The Result: Once it had enough data, the new model was better at predicting the future than the old rigid models.

5. Real-World Tests (The Case Studies)

They tested this new tool on two very different real-world problems:

Ocean Salinity (Trondheim Fjord): Here, the water conditions change drastically depending on where you are (non-stationarity). The new model, with its ability to change its "stretchiness" locally, was much better at predicting the salt levels.
Rainfall (USA): Here, the texture of the rain patterns was the key. The new model, with its ability to choose the perfect "roughness" (fractional smoothness), predicted rainfall better than the old models.

The Bottom Line

This paper doesn't just say "we have a new math formula." It says: "Nature is too complex for one-size-fits-all maps."

They built a flexible, smart tool that can:

Change its shape locally (Non-stationarity).
Choose the perfect level of roughness (Fractional smoothness).
Avoid getting confused by too much freedom (Penalization).

They proved that if you have enough data (around 500 points), this tool gives you a much clearer, more accurate picture of the world, whether you are looking at the ocean or the weather. They also made the code for this tool available so others can use it.

Technical Summary: Non-stationary Spatial Modeling Using Fractional SPDEs

Problem Statement
Gaussian Random Fields (GRFs) are standard tools for spatial modeling, typically relying on the Matérn covariance family which assumes stationarity and isotropy. However, many real-world applications, such as ocean salinity and precipitation modeling, violate these assumptions. Existing methods for non-stationary modeling, such as deformation or convolution approaches, often incur high computational costs due to dense covariance matrices. Furthermore, while the Stochastic Partial Differential Equation (SPDE) approach offers computational efficiency via sparse precision matrices, previous extensions have been limited: some incorporate spatially varying anisotropy but assume integer smoothness, while others allow fractional smoothness but assume stationarity or isotropy. Additionally, flexible non-stationary models are prone to overfitting without careful penalization, and the relative importance of fractional smoothness versus non-stationarity in specific applications remains under-explored.

Methodology
The authors propose a unified framework that combines fractional smoothness with spatially varying anisotropy within the SPDE approach.

Model Definition: The GRF $u(\cdot)$ is defined as the solution to a non-stationary SPDE:
$(\kappa(\mathbf{s}) - \nabla \cdot \mathbf{H}(\mathbf{s})\nabla)^\beta u(\mathbf{s}) = \tau(\mathbf{s})\mathcal{W}(\mathbf{s})$
where $\mathcal{W}$ is Gaussian white noise. The parameters $\kappa(\cdot)$ (range), $\tau(\cdot)$ (variance), and the positive-definite matrix $\mathbf{H}(\cdot)$ (anisotropy) vary spatially. The smoothness parameter $\beta$ (related to the Matérn smoothness $\nu = 2\beta - 1$ ) is allowed to be non-integer.
Parameterization:
- Anisotropy: The matrix $\mathbf{H}(\cdot)$ is parameterized using a spatially varying vector field $\mathbf{v}(\cdot)$ via a smooth bijection involving hyperbolic functions, ensuring the determinant is 1 and removing identifiability issues between $\mathbf{v}$ and $-\mathbf{v}$ .
- Non-stationarity: The log-transforms of the parameters ( $\log \kappa, \log \sigma, v_x, v_y$ ) are decomposed into stationary components and non-stationary components modeled as latent GRFs.
- Fractional Smoothness: The rational SPDE approach (Bolin and Kirchner, 2020) is employed to approximate the fractional operator $(\mathcal{L})^{-\beta}$ using rational functions (polynomials in the numerator and denominator).
Discretization and Computation:
- The domain is discretized using a finite element method (FEM) on a triangular mesh, resulting in sparse precision matrices.
- Automatic Differentiation (AD): To handle the complex log-likelihood derivatives required for parameter estimation (especially with the rational approximation), the authors utilize AD. This avoids the poor scaling of finite-difference methods as the number of parameters increases.
- Sparse AD Implementation: A custom library (sparsejax) is developed to perform sparse matrix operations within the AD framework, preserving sparsity patterns during backpropagation to ensure computational feasibility on both CPUs and GPUs.
Inference and Priors:
- A hierarchical Bayesian model is constructed.
- Penalized Complexity (PC) Priors: Priors are constructed to penalize model flexibility and prevent overfitting. For non-stationary components, penalties are set such that the probability of the parameter deviating significantly from its stationary counterpart is limited (e.g., to 5%).
- Optimization: Parameters are estimated via Maximum A Posteriori (MAP) using a two-stage gradient-based optimization (Adam followed by L-BFGS-B).

Key Contributions

Unified Model: The paper extends the SPDE approach to simultaneously incorporate spatially varying anisotropy and fractional smoothness, creating a more flexible class of models for spatial processes.
Computational Efficiency: By integrating automatic differentiation with sparse matrix operations, the authors enable efficient gradient-based optimization for high-dimensional non-stationary fractional models, overcoming the scaling limitations of previous numerical differentiation approaches.
Robust Inference: The development of specific priors for both stationary and non-stationary parameters ensures that the flexible model does not overfit, even with varying numbers of observation locations.
Empirical Assessment: The study systematically investigates the data requirements for reliably estimating fractional smoothness and non-stationarity.

Results

Simulation Study:
- Data Requirements: The study finds that a single realization containing at least 500 observations is required to reliably estimate fractional smoothness and non-stationarity in the considered scenarios.
- Overfitting: The proposed penalization successfully prevents overfitting across varying numbers of observation locations.
- Basis Functions: Increasing the number of basis functions for non-stationary components (e.g., from 8 to 24) did not degrade performance, suggesting that overspecifying the structure is not inherently detrimental if penalized correctly.
- Smoothness Estimation: Estimating the smoothness parameter $\nu$ provided only marginal improvements in predictive performance when the true data was generated from a stationary model.
Case Studies:
- Ocean Salinity (Trondheimsfjorden): Non-stationarity was found to be the dominant factor improving predictive performance (measured by RMSE and CRPS), while fractional smoothness had less impact.
- Precipitation (Conterminous U.S.): Conversely, fractional smoothness significantly improved predictions, whereas non-stationarity offered less relative benefit.
- Conclusion: The relative importance of fractional smoothness versus non-stationarity is application-dependent.

Significance and Claims
The authors claim that their approach provides a computationally robust tool for exploring the presence of fractional smoothness and non-stationarity in spatial data. The paper emphasizes that while the model is highly flexible, its practical utility relies on careful prior specification to avoid overfitting. The results suggest that practitioners should not assume one feature (e.g., non-stationarity) is universally superior; rather, the choice to include fractional smoothness or spatially varying anisotropy should be driven by the specific characteristics of the dataset and the application. The proposed method serves as both a predictive tool and an exploratory mechanism to diagnose the underlying spatial structure of environmental processes.