rSPDE

spatial statistics applied statistics Gaussian random fields SPDEs bayesian inference geostatistics

rSPDE: A Computational Framework for Rational Approximations of Fractional Stochastic Partial Differential Equations

rSPDE is an R package that implements rational approximations of fractional stochastic partial differential equations (SPDEs). It builds on the SPDE and finite element framework, originally developed to deliver a computationally efficient representation of spatial Gaussian fields whose dependence mimics the Matérn class. The Whittle–Matérn Gaussian fields are defined indirectly via a SPDE rather than by explicitly specifying a covariance function. The differential operator in that formulation encodes both how fast correlations decay over space and how smooth the field is. The name “Whittle–Matérn” is motivated by a seminal result of Peter Whittle (1963, Bulletin of the International Statistical Institute), who showed that when the spatial domain is the entire Euclidean space, the stationary solution of that SPDE coincides exactly with a Gaussian field having the classical Matérn covariance function. This equivalence forms the theoretical basis and provides the motivation for representing Matérn-like Gaussian fields through the SPDE approach, as discussed in modern treatments by Lindgren, Rue and Lindström (2011, Journal of the Royal Statistical Society, Series B) and, for a recent review, by Lindgren, Bolin and Rue (2021, Spatial Statistics). In settings where the domain is not the full Euclidean space, the same SPDE construction still defines a Whittle–Matérn field: whether the domain is a bounded subset of Euclidean space, a manifold, or a metric graph, the SPDE formulation generalizes naturally and retains many desirable properties. Importantly, by leveraging the finite element method, one obtains sparse precision (inverse covariance) matrices, which leads to computationally efficient and memory-efficient approximations suitable for large-scale inference. The rSPDE package implements the methods described in Bolin & Simas, enabling simulation, prediction, and Bayesian inference, including estimation of smoothness, within this general SPDE framework.

To illustrate the practical applicability of the rSPDE framework, we consider a real‐data example involving precipitation observations in southern Brazil. Figure 1 displays the spatial layout of the monitoring stations used in this example (hundreds of gauges across Paraná state). From these observations, we perform spatial prediction of precipitation fields using the rational approximation methods implemented in rSPDE, and the resulting posterior mean surface is visualized in Figure 2. Further technical details and implementation notes can be found in the vignette "inlabru implementation of the rational SPDE approach".

Figure 1: Observation locations for the precipitation dataset in Southern Brazil. These are used in the rSPDE case study for model fitting and spatial prediction.

Figure 2: Posterior mean of precipitation measurements in Southern Brazil using the rSPDE package with rational approximation and FEM discretization.

At its core, rSPDE implements rational approximations for fractional stochastic partial differential equation models in which a differential operator defines the spatial dependence and a parameter controls the smoothness of the resulting random field. Traditional SPDE formulations are limited to specific discrete smoothness levels, whereas rSPDE extends these models to any admissible degree of smoothness, thus enabling direct statistical inference on this parameter from data. As discussed by Stein in his book Interpolation of Spatial Data: Some Theory for Kriging, one of the key features of Gaussian processes with Matérn covariance functions is precisely the possibility of estimating the smoothness from observations. The rSPDE framework makes it possible to achieve the same goal for models formulated through the SPDE approach. Moreover, Bolin and Kirchner (2023, Bernoulli) show that in Whittle–Matérn models, having the correct smoothness value is sufficient to obtain asymptotically optimal estimators, emphasizing the practical and theoretical importance of inferring this parameter accurately.

The package implements two complementary methodological approaches developed through cutting-edge research. The operator-based rational approximation, introduced by Bolin and Kirchner (2020, Journal of Computational and Graphical Statistics), provides the foundational framework by combining finite element methods with rational approximations of fractional operators. This approach is particularly useful when considering fractional SPDE models with non-Gaussian noise. For Gaussian noise, it offers explicit convergence guarantees. Figure 3 illustrates a simulation example where a Whittle–Matérn field is recovered from noisy spatial observations using the rSPDE package. The rational approximation enables accurate interpolation, capturing the underlying spatial structure efficiently.

Figure 3: Simulated dataset and prediction from a Whittle–Matérn field with general smoothness. The left panel shows noisy observations, and the right panel displays the posterior mean prediction using rSPDE with rational approximation.=

The covariance-based rational approximation, developed by Bolin, Simas, and Xiong (2023, Journal of Computational and Graphical Statistics), represents a significant methodological advancement by directly approximating the covariance operator L^(-2β) rather than the solution itself. This innovation produces Gaussian Markov Random Field (GMRF) representations with sparse precision matrices, enabling seamless integration with R-INLA for fast and stable Bayesian inference. The method supports multiple rational approximation algorithms, including the BRASIL algorithm and Clenshaw-Lord Chebyshev-Padé approximations, providing flexibility for different computational requirements. Furthermore, this method is more numerically stable than the operator-based approach.

In addition to the FEM-based formulations, rSPDE also supports one-dimensional models that do not require finite element discretization. This feature is based on the rational approximation framework developed by Bolin, Mehandiratta, and Simas (2025, Journal of Machine Learning Research), which provides an efficient method for approximating stationary Gaussian models with Matérn covariance on intervals, that is, one-dimensional domains. This approach is particularly well suited for temporal or other one-dimensional spatial models, where the structure of the domain allows for accurate and computationally efficient inference without mesh construction.

A distinguishing feature of rSPDE is its comprehensive integration with Bayesian spatial modeling ecosystems. The package provides native interfaces to R-INLA and inlabru, enabling users to perform full Bayesian inference on all model parameters, including the smoothness parameter. Moreover, rSPDE integrates with the MetricGraph package to support spatial models defined on metric graphs (e.g. networks). In the MetricGraph vignette “Whittle–Matérn fields with general smoothness” the combination of finite element and rational approximation methods (as in rSPDE) is used to define Whittle–Matérn Gaussian random fields with general smoothness on metric graphs. In the vignette “Log-Gaussian Cox processes on metric graphs” the authors illustrate how rSPDE is used within the MetricGraph / INLA / inlabru pipeline for point process modeling over graphs. Through these integrations, rSPDE allows full Bayesian inference, including uncertainty quantification, posterior sampling, and model comparison, not only in Euclidean domains but also over complex network geometries.

The computational efficiency of rSPDE stems from its sophisticated approximation strategies. By constructing sparse precision matrices and leveraging finite element discretizations, the package achieves scalability to large spatial datasets while maintaining numerical stability. The covariance-based approach particularly excels in this regard, avoiding parameter-dependent projection matrices that would otherwise complicate inference procedures. Multiple rational approximation orders can be selected to balance computational cost with approximation accuracy, providing users with flexible control over the speed-accuracy trade-off.

The theoretical foundations of rSPDE are rigorously established through convergence analyses and explicit error bounds for both approximation schemes, which provide formal justification for the methods and guidance for practical implementation. These results ensure that the numerical approximations behave predictably as the discretization is refined, supporting reliable inference across different spatial settings. The framework’s flexibility enables applications in a wide range of domains that demand sophisticated spatial modeling. In environmental monitoring, it allows the analysis of pollution fields, climate variables, and ecological processes under flexible smoothness assumptions. In geostatistics and mining, it supports kriging with general Matérn covariance functions, enabling optimal interpolation of mineral concentrations and geological properties. In epidemiology, the same methodology can be used to study disease spread and spatially varying health outcomes, while in urban analytics it facilitates modeling of traffic flows, infrastructure, and demographic patterns characterized by complex spatial dependencies.

Building on the theoretical results above, the empirical accuracy of the rational SPDE approximation can be evaluated by comparing the resulting covariance function with the exact Matérn covariance. This is achieved by increasing the order m of the rational approximation, which corresponds to the degree of the approximating rational function applied to the differential operator. As expected, the error decreases rapidly with increasing m. For instance, using the same mesh and evaluating the covariance at fixed locations, the norm of the difference between the true and approximate covariances drops from approximately 1.01 for m = 1 to just 0.017 for m = 4. This exponential decay in error highlights one of the key advantages of the rational SPDE framework: it provides high accuracy even for relatively low orders, which helps balance precision and computational efficiency. For further details see the vignette “Operator-based rational approximation”. Figure 4 provides a direct comparison between the exact Matérn covariance function and its rational approximation for order m = 1. Despite the low approximation order, the curves are nearly indistinguishable, illustrating the high accuracy of the rational SPDE method even with minimal computational effort.

Figure 4: Comparison between the exact Matérn covariance function and its rational approximation with m = 1. Even at the lowest approximation order, the rational method closely reproduces the target covariance structure.

The package supports both stationary and non-stationary spatial processes, spatio-temporal extensions, and point process models. These features make it a unified computational platform for a broad class of spatial statistical applications as well as applications in time series. Implementation follows best practices in scientific software development: the codebase is extensively tested and validated for consistency across systems, the documentation is comprehensive, and the vignettes cover both methodological background and applied case studies. Integration with other spatial modeling tools is straightforward, allowing researchers to incorporate rSPDE into existing workflows without substantial modification.

To facilitate learning and adoption, the package provides detailed documentation, reproducible examples, and vignettes illustrating different modeling approaches and applications, including case studies based on real datasets such as precipitation measurements from southern Brazil as shown in Figures 1 and 2. These materials guide users from introductory concepts to advanced modeling scenarios, making the methodology accessible to both applied practitioners and researchers. Further information and learning resources are available on the rSPDE homepage rSPDE homepage. Overall, rSPDE bridges theoretical developments in fractional SPDEs with efficient computational tools for Bayesian spatial inference. By enabling estimation of general smoothness parameters while maintaining numerical stability and scalability, it broadens the scope of SPDE-based modeling and provides a robust foundation for modern spatial statistical analysis.

rSPDE

David Bolin

Alexandre de Bustamante Simas

Share

Related Links

Stochastic Processes and Mathematical Statistics (StochProc)

rSPDE

Related People

David Bolin

Alexandre de Bustamante Simas

Share

Related Links