Covariance–Based Rational Approximations of Fractional SPDEs for Computationally Efficient Bayesian Inference

ABSTRACT

The stochastic partial differential equation (SPDE) approach is widely used for modeling large spatial datasets. It is based on representing a Gaussian random field u on ${\mathbb{R}}^d$ as the solution of an elliptic SPDE $L^{\beta }u=\mathcal{W}$ where L is a second-order differential operator, $2\beta \in \mathbb{N}$ is a positive parameter that controls the smoothness of u and $\mathcal{W}$ is Gaussian white noise. A few approaches have been suggested in the literature to extend the approach to allow for any smoothness parameter satisfying $\beta >d/4$. Even though those approaches work well for simulating SPDEs with general smoothness, they are less suitable for Bayesian inference since they do not provide approximations which are Gaussian Markov random fields (GMRFs) as in the original SPDE approach. We address this issue by proposing a new method based on approximating the covariance operator $L^{-2\beta }$ of the Gaussian field u by a finite element method combined with a rational approximation of the fractional power. This results in a numerically stable GMRF approximation which can be combined with the integrated nested Laplace approximation (INLA) method for fast Bayesian inference. A rigorous convergence analysis of the method is performed and the accuracy of the method is investigated with simulated data. Finally, we illustrate the approach and corresponding implementation in the R package rSPDE via an application to precipitation data which is analyzed by combining the rSPDE package with the R-INLA software for full Bayesian inference. Supplementary materials for this article are available online.

KEYWORDS:

• Gaussian Markov random field
• Gaussian process
• Latent Gaussian model
• Spatial statistics
• SPDE
• R-INLA

1 Introduction

Handling many observations from a Gaussian random field in spatial statistics can be challenging since the related computational tasks involve factorizations of large covariance matrices which are usually dense. This is often referred to as the “big N problem” (Banerjee, Carlin, and Gelfand 2015), and various approaches have been suggested to handle the computational issues (see, e.g., Heaton et al. 2019, for a recent comparison). One of the most widely used methods is the SPDE approach by Lindgren et al. (2011). This is based on the fact that a centered Gaussian random field u on the spatial domain $\mathcal{D}={\mathbb{R}}^d$ with an isotropic Matérn covariance function (Matérn 1960), (1) $$\varrho (\mathit{s})=\frac{{\sigma }^2}{2^{\nu -1}\Gamma (\nu )}{(\kappa \Vert \mathit{s}\Vert )}^{\nu }{K}_{\nu }(\kappa \Vert \mathit{s}\Vert ),\hspace{1em}\mathit{s}\in {\mathbb{R}}^d,$$(1) can be represented as a solution to the stochastic partial differential equation (SPDE) (2) $${({\kappa }^2-\Delta )}^{\beta }(\tau u)=\mathcal{W}\hspace{1em}\text{in}\mathcal{D}.$$(2)

Here, $\Gamma (·)$ is the Gamma function, $K_{\nu }$ is a modified Bessel function of the second kind, Δ is the Laplace operator, and $\mathcal{W}$ is Gaussian white noise. The parameter $\kappa >0$ controls the practical correlation range, ${\sigma }^2$ is the variance of u, ${\tau }^2=\Gamma (\nu )/({\sigma }^2{\kappa }^{2\nu }{(4\pi )}^{d/2}\Gamma (\nu +d/2))$, and the fractional power β is related to the smoothness parameter $\nu >0$ via the relation $2\beta =\nu +d/2$ (Whittle 1963). Lindgren et al. (2011) used this representation to construct computationally efficient Gaussian Markov Random Field (GMRF) approximations of Gaussian Matérn fields by considering the SPDE on a bounded domain $\mathcal{D}$, restricting the smoothness to $2\beta \in \mathbb{N}$, and then performing a finite element method (FEM) discretization.

The SPDE approach has become widely used in applications, and has initiated a great number of extensions and generalizations (Lindgren, Bolin, and Rue 2022). The reason for this is not only the computational benefits, but also that it provides a flexible framework for defining more sophisticated models for spatial data. It, for example, facilitates the construction of nonstationary Gaussian random fields by allowing the parameters κ and τ to be spatially varying (Lindgren et al. 2011; Fuglstad et al. 2015), and allows for the construction of Matérn-like random fields on more general manifolds by defining such fields via the SPDE (2) posed on the manifold (Lindgren et al. 2011; Bolin and Lindgren 2011).

One of the main criticisms of the SPDE approach is the requirement $2\beta \in \mathbb{N}$, which restricts the possible values of the corresponding smoothness parameter ν of the Matérn covariance function. Given the importance of ν when performing prediction, as shown by Stein (1999) and Bolin and Kirchner (2023), several methods for removing the restriction of $2\beta \in \mathbb{N}$ have been proposed. Lindgren et al. (2011, Author’s response) proposed to construct a GMRF approximation by approximating the spectrum of a Gaussian Matérn field by a spectrum that is a reciprocal of a polynomial. This method is applicable for stationary models but it can not be applied to nonstationary models, and it has a fixed accuracy which may not be sufficient for certain applications. Bolin, Kirchner, and Kovács (2020) proposed combining the FEM approximation of Lindgren et al. (2011) with a quadrature approximation of the fractional operator to obtain a numerical method that works for any $\beta >d/4$ and can be made arbitrarily accurate. That work also provided a theoretical convergence analysis of the method, which was extended in Bolin, Kirchner, and Kovács (2018) and Herrmann, Kirchner, and Schwab (2020). Bolin and Kirchner (2020) later proposed a different type of approximation referred to as the rational SPDE approach, which has a lower computational cost.

Even though the methods that work for nonstationary models with general smoothness are computationally efficient, they are much less used than the standard SPDE approach for statistical applications. The reason for this is that nonfractional SPDE models work in combination with the integrated nested Laplace approximation (INLA) method (Rue, Martino, and Chopin 2009) and are implemented in the widely used R-INLA (Lindgren and Rue 2015) R (R Core Team 2022) package. This software facilitates including SPDE-based models in general Bayesian latent Gaussian models, and the great majority of all applications of the SPDE approach have been done via this software.

Unfortunately, the methods of Bolin, Kirchner, and Kovács (2020) and Bolin and Kirchner (2020) provide approximations which are not compatible with R-INLA. The reason is that the methods do not yield a Markov approximation, so the precision matrix obtained from the approximations are not sparse. The covariance matrix of the approximations are of the form $\mathbf{P}{\mathbf{Q}}^{-1}\mathbf{P}$, where both P and Q are sparse matrices that depend on the parameters of the model. To achieve a sparse precision matrix, which is necessary for R-INLA, Bolin and Kirchner (2020) showed that one can work with a latent model with sparse precision matrix Q if the projection matrix A, which connects the locations of the mesh for the FEM approximation and the observation locations (see Section 2 for details), is adjusted to $\widehat{\mathbf{A}}=\mathbf{\text{AP}}$. This matrix, however, depends on the model parameters, which is not allowed in R-INLA.

The main goal of this work is to solve this problem by proposing an alternative rational approximation. The main idea is to approximate the covariance operator of the random field directly, instead of first approximating the solution u and then deriving the corresponding covariance operator. This provides an approximation suitable for R-INLA, which is more numerically stable than that of the original rational SPDE approach. The proposed method is implemented in the R package rSPDE (Bolin and Simas 2023), which is available on CRAN and has an interface to R-INLA. Using the package, we show that the proposed method facilitates full Bayesian inference of all model parameters, including β, for latent Gaussian models based on fractional SPDEs.

The outline of the article is as follows. In Section 2, we give an overview of the model structure of the proposed approximation and show how it can be used for computationally efficient inference. The mathematical details and justifications of the method are provided in Sections 3–5. Specifically, in Section 3, we introduce the generalized Whittle–Matérn fields, which contain most of the previously proposed nonstationary SPDE-based Gaussian random fields as special cases, and for which our proposed method is applicable. In that section, we also provide the details of the FEM approximations. The new covariance-based rational approximation is introduced in Section 4, where we also prove that it provides an approximation of the covariance function of the generalized Whittle–Matérn field with an explicit rate of convergence in the L₂-norm. In Section 5, we show that the covariance-based rational approximation can be represented as a GMRF, and illustrate how this can be used for statistical inference. Some of the details of the rSPDE implementation are discussed in Section 6, and a comparison in terms of the accuracy of approximating covariance function by our method and some other methods is provided in Section 7. An application to modeling of precipitation data is presented in Section 8 and the article concludes with a discussion in Section 9. Finally, the supplementary materials contain seven appendices which provide further technical details and proofs.

2 Overview of the Approximation Strategy

As mentioned in the introduction, the main idea behind our strategy is to directly approximate the covariance operator of the random field. In this section we show the structure of the resulting approximation and also provide an illustration on how it can be used for inference in a simple application. More details will be given in later sections. The covariance-based rational approximation of the Whittle–Matérn field $u(\mathit{s})$ defined in (2), whose covariance operator is ${({\kappa }^2-\Delta )}^{-2\beta }$, uses a combination of the finite element method and rational approximations in order to approximate $u(\mathit{s})$ as $u_n(\mathit{s})={\sum }_{j=1}^n{w}_j{\phi }_j(\mathit{s})$, where ${\left\{w_j\right\}}_{j=1}^n$ are stochastic weights and ${\left\{{\phi }_j\right\}}_{j=1}^n$ are FEM basis functions. With our approximation, $\mathit{w}={[w_1,\dots ,w_n]}^{\top }$ can be expressed as a sum of m + 1 independent GMRFs ${\mathit{x}}_i={(x_{i1}\cdots x_{\mathit{\text{in}}})}^{\top }$ with sparse precision matrices: (3) $$\mathit{w}=\sum _{i=1}^{m+1}{\mathit{x}}_i,\hspace{1em}\text{where}\hspace{1em}{\mathit{x}}_i\sim N\left(\mathbf{0},{\mathit{Q}}_i^{-1}\right).$$(3)

Any linear predictor in R-INLA has this form, which means that we can perform inference in a computationally efficient manner based on the covariance-based rational approximation by using the same ideas as are used in R-INLA. For example, suppose that we observe $y_1,\dots ,y_N,N\in \mathbb{N}$, and spatial locations ${\mathit{s}}_1,\dots ,{\mathit{s}}_N\in {\mathbb{R}}^d,$ where (4) $$y_i=u({\mathit{s}}_i)+{ϵ}_i,\hspace{1em}i=1,\dots ,N,$$(4) and $\mathit{ϵ}={[{ϵ}_1,\dots ,{ϵ}_N]}^{\top }\sim N(\mathbf{0},{\mathit{Q}}_{\mathit{ϵ}}^{-1})$ for some sparse matrix ${\mathit{Q}}_{\mathit{ϵ}}$, such as ${\mathit{Q}}_{\mathit{ϵ}}=\frac{1}{{\sigma }_{ϵ}^2}{\mathit{I}}_N$ if we have independent measurement noise. Defining $\mathit{y}={[y_1,\dots ,y_N]}^{\top }$, (4) can be written in matrix form as $\mathit{y}=\mathbf{A}\mathit{w}+\mathit{ϵ}$, where A is the projector matrix with elements $A_{\mathit{\text{ij}}}={\phi }_j({\mathit{s}}_i)$. Let $\mathit{X}={[{\mathit{x}}_1^{\top },\dots ,{\mathit{x}}_{m+1}^{\top }]}^{\top }$. Then, the precision matrix of X is the block diagonal matrix (5) $$\mathit{Q}=\text{diag}({\mathit{Q}}_1,\dots ,{\mathit{Q}}_{m+1}).$$(5)

Writing the model in terms of the weights X allows us to equivalently write the model as $\mathit{y}=\overline{\mathit{A}}\mathit{X}+\mathit{ϵ}$ where $\overline{\mathit{A}}$ is a block matrix of size $N\times n(m+1)$ obtained by combining m + 1 copies of A as $\overline{\mathit{A}}=\left[\begin{array}{ccc}\mathit{A}& \cdots & \mathit{A}\end{array}\right]$. Thus, $\mathit{y}|\mathit{X}\sim N(\overline{\mathit{A}}\mathit{X},{\mathit{Q}}_{\mathit{ϵ}}^{-1})$ and $\mathit{X}\sim N(\mathbf{0},{\mathit{Q}}^{-1})$, where Q is given in (5). Standard results for latent Gaussian models then give us that the posterior distribution of X is $\mathit{X}|\mathit{y}\sim N({\mathit{\mu }}_{\mathit{X}|\mathit{y}},{\mathit{Q}}_{\mathit{X}|\mathit{y}}^{-1})$, where (6) $${\mathit{\mu }}_{\mathit{X}|\mathit{y}}={\mathit{Q}}_{\mathit{X}|\mathit{y}}^{-1}{\overline{\mathit{A}}}^{\top }{\mathit{Q}}_{\mathit{ϵ}}\mathit{y}\hspace{1em}\text{and}\hspace{1em}{\mathit{Q}}_{\mathit{X}|\mathit{y}}={\overline{\mathit{A}}}^{\top }{\mathit{Q}}_{\mathit{ϵ}}\overline{\mathit{A}}+\mathit{Q}.$$(6)

Finally, we can obtain the marginal likelihood, $\mathcal{l}(\mathit{y})$, of y as (7) $$\begin{array}{cc}2\mathcal{l}(\mathit{y})=\hfill & \hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left|\mathit{Q}\right|+\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left|{\mathit{Q}}_{\mathit{ϵ}}\right|-\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left|{\mathit{Q}}_{\mathit{X}|\mathit{y}}\right|-{\mathit{\mu }}_{\mathit{X}|\mathit{y}}^{\top }\mathit{Q}{\mathit{\mu }}_{\mathit{X}|\mathit{y}}\hfill \\ \hfill & -{(\mathit{y}-\overline{\mathit{A}}{\mathit{\mu }}_{\mathit{X}|\mathit{y}})}^{\top }{\mathit{Q}}_{\mathit{ϵ}}(\mathit{y}-\overline{\mathit{A}}{\mathit{\mu }}_{\mathit{X}|\mathit{y}})-n\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}(2\pi ).\hfill \end{array}$$(7)

The sparsity of Q is essential for computation. For instance, evaluating $\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}\left|\mathit{Q}\right|$ in the likelihood can be done efficiently based on sparse Cholesky decomposition (Rue and Held 2005). Sparsity of Q also facilitates computationally efficient sampling of w, and hence of u. See Appendix E for further details on the methods for sampling and likelihood evaluation.

3 Whittle–Matérn Fields and FEM Approximation

In this section we introduce the class of fractional-order SPDEs we are interested in as well as their FEM approximations. The model assumptions are presented in Section 3.1, and in Section 3.2 we introduce the FEM approximations and study their convergence.

Let us begin by introducing the notation that will be needed later on. Given a bounded domain $\mathcal{D}\subset {\mathbb{R}}^d$, $d\in \{1,2,3\}$, we denote by $L_2(\mathcal{D})$ the Lebesgue space of square-integrable real-valued functions endowed with the inner product ${(\varphi ,\psi )}_{L_2(\mathcal{D})}={\int }_{\mathcal{D}}\varphi (\mathit{x})\psi (\mathit{x})d\mathit{x}$. We denote the Sobolev space of order k by $H^k(\mathcal{D})$: $$H^k(\mathcal{D})=\left\{w\in L_2(\mathcal{D}):D^{\gamma }w\in L_2(\mathcal{D}),\forall \gamma \in {\mathbb{N}}^d,\left|\gamma \right|\le k\right\},$$where we are using the multiindex notation for the differential operator $D^{\gamma }$, and ${(·,·)}_{H^k(\mathcal{D})}$ is the Sobolev inner product: $${(u,v)}_{H^k(\mathcal{D})}=\sum _{\gamma \in {\mathbb{N}}^d:\left|\gamma \right|\le k}{(D^{\gamma }u,D^{\gamma }v)}_{L_2(\mathcal{D})}.$$

We denote by $H_0^1(\mathcal{D})$ the closure of $C_c^{\infty }(\mathcal{D})$ in $H^1(\mathcal{D})$, where $C_c^{\infty }(\mathcal{D})$ is the set of infinitely differentiable functions with compact support on $\mathcal{D}$. Additional notations needed for the theoretical analysis are given in Appendix A.

3.1. Model Assumptions

We are interested in the class of Gaussian random fields on $\mathcal{D}$ that can be represented as solutions to SPDEs of the form (8) $$L^{\beta }(\tau u)=\mathcal{W}\hspace{1em}\text{in}\mathcal{D},$$(8) where $L^{\beta }$ is a fractional power (in the spectral sense) of a second-order elliptic differential operator L which determines the covariance structure of u, $\tau >0$ is a constant parameter, and $\mathcal{W}$ is Gaussian white noise on $L_2(\mathcal{D})$. We have the following assumptions on $\mathcal{D}$:

Assumption 1.

The domain $\mathcal{D}$ is an open, bounded, convex polytope with closure $\overline{\mathcal{D}}$.

Under Assumption 1, we may define $$H_{\mathcal{N}}^2(\mathcal{D})=\{w\in H^2(\mathcal{D}):\partial w/\partial \nu =0\text{on}\partial \mathcal{D}\},$$where ν is the outward unit normal vector to $\partial \mathcal{D}$. Indeed, the expression $\partial w/\partial \nu =0$ on $\partial \mathcal{D}$ makes sense since the trace of Dw is well-defined in this case (see, e.g., Evans and Gariepy 2015, Theorem 4.6). Let us now describe the assumptions on the differential operator L:

Assumption 2.

The operator L is given in divergence form by $\mathit{\text{Lu}}=-\nabla ·(\mathit{H}\nabla u)+{\kappa }^{2}u$, and is equipped either with homogeneous Dirichlet or Neumann boundary conditions. Furthermore, the function $\mathit{H}:\overline{\mathcal{D}}\to {\mathbb{R}}^{d\times d}$ is symmetric, Lipschitz continuous and uniformly positive definite, and $\kappa :\mathcal{D}\to \mathbb{R}$ is an essentially bounded function, that is, $$\text{ess} {\text{sup}}_{\mathrm{x}\in \mathcal{D}}\kappa (x)=\mathrm{\text{inf}}\{a\in \mathbb{R}:\lambda (\{x:\kappa (x)>a\})=0\}<\infty .$$

Under Neumann boundary conditions, we additionally require that $$\text{ess} {\text{inf}}_{\mathrm{x}\in \mathcal{D}}\kappa (x)=\mathrm{\text{sup}}\{b\in \mathbb{R}:\lambda (\{x:\kappa (x)<b\})=0\}\ge {\kappa }_0>0,$$where λ is the Lebesgue measure on $\mathcal{D}$.

The SPDE (8) under Assumptions 1 and 2 defines a class of models that have previously been considered by Bolin, Kirchner, and Kovács (2020), Cox and Kirchner (2020), Herrmann, Kirchner, and Schwab (2020), and Bolin and Kirchner (2020) and is referred to as generalized Whittle–Matérn fields. It contains many previously proposed nonstationary SPDE-based spatial Gaussian random field models as special cases, such as those by Lindgren et al. (2011), Fuglstad et al. (2015), Fuglstad et al. (2019), and Hildeman, Bolin, and Rychlik (2021), and the method that we later introduce thus also applies to those models and their fractional extensions.

In the case of Dirichlet boundary conditions, define the space $V=H_0^1(\mathcal{D})\subset L_2(\mathcal{D})$, and in the case of Neumann boundary conditions let $V=H^1(\mathcal{D})\subset L_2(\mathcal{D})$. Then, under Assumptions 1 and 2, L induces the following continuous and coercive bilinear form on V: (9) $$a_L(v,u)={(\mathit{H}\nabla u,\nabla v)}_{L_2(\mathcal{D})}+{({\kappa }^{2}u,v)}_{L_2(\mathcal{D})},\hspace{1em}u,v\in V.$$(9)

Remark 1.

Under Assumptions 1 and 2, if $f\in L_2(\mathcal{D})$, then there exists a unique solution u of Lu = f and the operator L is $H^2(\mathcal{D})$-regular, that is, $u\in H^2(\mathcal{D})\cap H_0^1(\mathcal{D})$ under Dirichlet boundary conditions, whereas under Neumann boundary conditions, we have $u\in H_{\mathcal{N}}^2(\mathcal{D})$. See, for instance, (Grisvard 2011, Theorem 3.2.1.2) for Dirichlet boundary conditions or (Grisvard 2011, Theorem 3.2.1.3) for Neumann boundary conditions.

By Remark 1, specifically by the existence and uniqueness of the solution to the equation Lu = f, we can define the inverse operator $L^{-1}:L_2(\mathcal{D})\to L_2(\mathcal{D})$. By Rellich-Kondrachov theorem (Evans and Gariepy 2015, Theorem 4.11), $L^{-1}$ is a compact operator, and observe that $L^{-1}$ is self-adjoint, see Appendix A for a justification. Hence, by the spectral theorem for self-adjoint and compact operators, there exists an orthonormal basis ${\left\{e_j\right\}}_{j\in \mathbb{N}}$ in $L_2(\mathcal{D})$ formed by eigenvectors of L whose eigenvalues ${\left\{{\lambda }_j\right\}}_{j\in \mathbb{N}}$ are nonnegative and can be arranged in a nondecreasing order.

Remark 2.

Under Assumptions 1 and 2, the operator L satisfies the Weyl’s law, that is, there exist $c,C>0$ such that for every $j\in \mathbb{N}$, $c{j}^{2/d}\le {\lambda }_j\le C{j}^{2/d}.$ See Davies (1995, Theorem 6.3.1) for the Dirichlet case. For the Neumann case, the Weyl’s law holds for the case in which H is a constant diagonal matrix (Fedosov 1963, 1964), in particular, it holds for the Neumann Laplacian. The result for a general H satisfying Assumption 2 is a direct consequence of the Weyl’s law for the Neumann Laplacian together with Proposition B.2 in Appendix B and the min–max principle.

Our goal is to obtain approximations of the covariance operator $L^{-2\beta }$ of the Gaussian random field u which solves equation (8)(8) $$L^{\beta }(\tau u)=\mathcal{W}\hspace{1em}\text{in}\mathcal{D},$$(8) . Let $${\varrho }^{\beta }(x,y)=\sum _{j=1}^{\infty }{\lambda }_j^{-2\beta }{e}_j(x)e_j(y).$$

Then, one can readily check, by Steinwart and Scovel (2012, Theorem 3.10), that the covariance operator $L^{-2\beta }$ is a kernel operator, with kernel ${\varrho }^{\beta }(·,·)$. That is, for any $f\in L_2(\mathcal{D})$, we have $(L^{-2\beta }f)(x)={\int }_{\mathcal{D}}{\varrho }^{\beta }(x,y)f(y)\mathrm{d}y$ for a.e. $x\in \mathcal{D}$. It is well-known that there exists a centered square-integrable Gaussian random field u that solves (8) if, and only if, its covariance operator, $L^{-2\beta }$, has finite trace (Lototsky and Rozovsky 2017, Theorem 3.2.5). Under Assumptions 1 and 2, one can use Weyl’s law (Remark 2) to show that $L^{-2\beta }$ has finite trace if, and only if, $\beta >d/4$. Hence, if $\beta >d/4$, then u is a centered square-integrable Gaussian random field with covariance function ${\varrho }^{\beta }(x,y)=E[u(x)u(y)],$ where the equality holds for a.e. $(x,y)\in \mathcal{D}\times \mathcal{D}.$

3.2. Finite Element Approximation

The goal is now to provide a convergence analysis for FEM approximations of the covariance operator $L^{-2\beta }$. Let us start by describing the setup we will use.

Assumption 3.

Let $V_h\subset V$ be a finite element space that is spanned by a set of continuous piecewise linear basis functions ${\left\{{\phi }_j\right\}}_{j=1}^{n_h}$ (see Appendix D), with $n_h\in \mathbb{N}$, defined with respect to a triangulation ${\mathcal{T}}_h$ of $\overline{\mathcal{D}}$ indexed by the mesh width $h:={\mathrm{\text{max}}}_{T\in {\mathcal{T}}_h}{h}_T$, where $h_T:=\text{diam}(T)$ is the diameter of the element $T\in {\mathcal{T}}_h$. We assume that the family ${({\mathcal{T}}_h)}_{h\in (0,1)}$ of triangulations inducing the finite-dimensional subspaces ${(V_h)}_{h\in (0,1)}$ of V is quasi-uniform, that is, there exist constants $K_1,K_2>0$ such that ${\rho }_T\ge K_1{h}_T$ and $h_T\ge K_{2}h$ for all $T\in {\mathcal{T}}_h$ and all $h\in (0,1)$. Here, ${\rho }_T>0$ is the radius of the largest ball inscribed in $T\in {\mathcal{T}}_h$.

We are now in a position to describe the FEM discretization of the model (8). Let $L_h:V_h\to V_h$ be defined in terms of the bilinear form a_L as its restriction to $V_h\times V_h$: $${(L_h{\varphi }_h,{\psi }_h)}_{L_2(\mathcal{D})}=a_L({\varphi }_h,{\psi }_h),\hspace{1em}{\varphi }_h,{\psi }_h\in V_h.$$

Note that L_h is a positive-definite, symmetric, linear operator on the finite-dimensional space V_h. Hence, we may arrange the eigenvalues of L_h as $0<{\lambda }_{1,h}\le {\lambda }_{2,h}\le \cdots \le {\lambda }_{n_h,h},$ with corresponding eigenvectors ${\left\{e_{j,h}\right\}}_{j=1}^{n_h}$ which are orthonormal in $L_2(\mathcal{D})$. Let ${\mathcal{W}}_h$ denote Gaussian white noise on V_h. That is, there exist independent standard Gaussian random variables ${\xi }_1,\dots ,{\xi }_{n_h}$ such that ${\mathcal{W}}_h={\sum }_{j=1}^{n_h}{\xi }_j{e}_{j,h}$. Then, we refer to the following SPDE on V_h as the discrete model of (8): (10) $$L_h^{\beta }{u}_h={\mathcal{W}}_h.$$(10)

Let u_h be a solution of (10), then the covariance operator of u_h is given by $L_h^{-2\beta }$, and $${\varrho }_h^{\beta }(x,y)=\sum _{j=1}^{n_h}{\lambda }_{j,h}^{-2\beta }{e}_{j,h}(x)e_{j,h}(y),\hspace{1em}\text{for a}.\mathrm{e}.(x,y)\in \mathcal{D}\times \mathcal{D},$$is the corresponding covariance function. We have the following result regarding the convergence of the FEM approximation ${\varrho }_h^{\beta }$ to the exact covariance function ${\varrho }^{\beta }$ in the $L_2(\mathcal{D}\times \mathcal{D})$-norm defined by $\left|\right|f|{|}_{L_2(\mathcal{D}\times \mathcal{D})}^2={\int }_{\mathcal{D}}{\int }_{\mathcal{D}}f{(x,y)}^{2}d\mathit{\text{xdy}}$. The proof is given in Appendix B.

Proposition 1.

Under Assumptions 1, 2, and 3, for each $\beta >d/4$ and each $\epsilon >0$, we have (11) $$\left|\right|{\varrho }^{\beta }-{\varrho }_h^{\beta }|{|}_{L_2(\mathcal{D}\times \mathcal{D})}{\lesssim }_{\epsilon ,\beta ,\mathit{H},\kappa ,\mathcal{D}}{h}^{\mathrm{\text{min}}\{4\beta -d/2-\epsilon ,2\}}.$$(11)

Here, and in the remainder of the article, the notation $A{\lesssim }_{{\theta }_1,\dots ,{\theta }_k}B$, where $k\in \mathbb{N}$, means that there exists a constant C depending on ${\theta }_1,\dots ,{\theta }_k$ (${\theta }_i,i=1,\dots ,k$, can be a parameter, a function, a domain, etc.) such that $A\le \mathit{\text{CB}}.$

Remark 3.

Cox and Kirchner (2020, Theorem 1) proved the bound (11) in the case of homogeneous Dirichlet boundary conditions. They did not provide a bound for the case of homogeneous Neumann boundary conditions. Proposition 1 arrives at the same bound for the Neumann case. For this, we additionally require that $\text{ess} {\text{inf}}_{\mathrm{x}\in \mathcal{D}}\kappa (\mathrm{x})\ge {\kappa }_0>0$ and that the domain $\mathcal{D}$ is a convex polytope in the Neumann case. As far as we know, this is a new result. The key step in the proof is to obtain an analogous result to Cox and Kirchner (2020, Lemma 2), which is given by Proposition B.2 in Appendix B.

4 Rational Approximation

Having introduced the FEM approximation, we are now ready to define the complete approximation of the covariance operator of the generalized Whittle–Matérn fields. The approximation is obtained by combining a rational approximation of the fractional power of the covariance operator with the FEM approximation. We begin by introducing the method and then provide a theoretical justification by showing an explicit rate of convergence of the approximate covariance function to the correct one in the $L_2(\mathcal{D}\times \mathcal{D})$-norm.

In Bolin and Kirchner (2020), the authors obtained an approximation of the solution to (8), which also implicitly defines an approximation of the corresponding covariance operator. However, as we have previously mentioned, this results in an approximation that is not implementable in R-INLA. Also, for statistical applications there is usually no need to have an approximation of the solution itself, since only the corresponding distribution matters for inference. With this in mind, we propose to directly approximate the covariance operator $L^{-2\beta }$. To this end, we first split $L_h^{-2\beta }=L_h^{-\left\{2\beta \right\}}{L}_h^{-\lfloor 2\beta \rfloor }$, where $\left\{x\right\}=x-\lfloor x\rfloor$ is the fractional part of x. Then, we approximate $L_h^{-\left\{2\beta \right\}}$ with a rational approximation. This yields an approximation (12) $$L_h^{-2\beta }\approx L_{h,m}^{-2\beta }:=L_h^{-\lfloor 2\beta \rfloor }p(L_h^{-1})q{(L_h^{-1})}^{-1}.$$(12)

Here, $p(L_h^{-1})={\sum }_{i=0}^m{a}_i{L}_h^{m-i}$ and $q(L_h^{-1})={\sum }_{j=0}^m{b}_j{L}_h^{m-j}$ are polynomials obtained from a rational approximation of order m of the real-valued function $f(x)=x^{\left\{2\beta \right\}}$. That is, $$x^{\left\{2\beta \right\}}\approx \frac{\sum _{i=0}^m{a}_i{x}^i}{\sum _{i=0}^m{b}_i{x}^i}.$$

Specifically, to obtain ${\left\{a_i\right\}}_{i=0}^m$ and ${\left\{b_i\right\}}_{i=0}^m$, we approximate the function $f(x)=x^{\left\{2\beta \right\}}$ on the interval $[{\lambda }_{n_h,h}^{-1},{\lambda }_{1,h}^{-1}]$, which covers the spectrum of $L_h^{-1}$. The coefficients are computed as the best rational approximation in the $L_{\infty }$-norm, which, for example, can be obtained via the second Remez algorithm (Remez 1934) or by the recent, and more stable, BRASIL algorithm (Hofreither 2021). See Appendix F for details about this algorithm and a justification for the choice of using the best rational approximation in $L_{\infty }$-norm.

By defining the covariance function $${\varrho }_{h,m}^{\beta }(x,y)=\sum _{j=1}^{n_h}{\lambda }_{j,h}^{-\lfloor 2\beta \rfloor }p({\lambda }_{j,h}^{-1})q{({\lambda }_{j,h}^{-1})}^{-1}{e}_{j,h}(x)e_{j,h}(y),$$for a.e. $(x,y)\in \mathcal{D}$, we have that ${\varrho }_{h,m}^{\beta }$ is the kernel of the covariance operator $L_{h,m}^{-2\beta }$. There are two sources of errors when we consider ${\varrho }_{h,m}^{\beta }$ as an approximation of the true covariance function ${\varrho }^{\beta }$ of the generalized Whittle–Matérn field: the FEM approximation and the rational approximation. The following proposition, whose proof is given in Appendix B, shows that we have control of these two sources of errors via the FEM mesh width h and the order of the rational approximation m.

Proposition 2.

Let $\beta >d/4$. Under Assumptions 1, 2, and 3, for every $\epsilon >0$ and for sufficiently small h, we have: (13) $$\left|\right|{\varrho }_{h,m}^{\beta }-{\varrho }^{\beta }|{|}_{L_2(\mathcal{D}\times \mathcal{D})}{\lesssim }_{\epsilon ,\beta ,\mathit{H},\kappa ,\mathcal{D}}{h}^{\mathrm{\text{min}}\{4\beta -d/2-\epsilon ,2\}}+{\mathbb{1}}_{2\beta \notin \mathbb{N}}{h}^{-d/2}{e}^{-2\pi \sqrt{\left\{2\beta \right\}m}}.$$(13)

Remark 4.

We can calibrate the accuracy of the rational approximation with the finite element error by choosing $m\in \mathbb{N}$ such that $m=\lceil {(\mathrm{\text{min}}\{4\beta -d/2-\epsilon ,2\}+d/2)}^2\frac{{(\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}h)}^2}{4{\pi }^2\left\{2\beta \right\}}\rceil .$ This ensures that the rate of convergence in (13) is $\mathrm{\text{min}}\{4\beta -d/2-\epsilon ,2\}$. See Section 7 for further details on the choice of m.

5 GMRF Representation

The goal of this section is to obtain a sparse matrix representation of the precision operator of the rational approximation from the previous section, so that the methods in Section 2 can be used for computationally efficient sampling and likelihood evaluation.

The solution u_h in (10) at spatial location s can be represented as $u_h(\mathit{s})={\sum }_{j=1}^{n_h}{w}_j{\phi }_j(\mathit{s})$, where ${\left\{w_j\right\}}_{j=1}^{n_h}$ are stochastic weights and ${\left\{{\phi }_j\right\}}_{j=1}^{n_h}$ are the piecewise linear finite element basis functions. We will now show how to represent $\mathit{w}={[w_1,\dots ,w_{n_h}]}^{\top }$ as a sum of independent GMRFs, each with a sparse precision matrix. The key step is to apply a partial fraction decomposition in (12): (14) $$L_{h,m}^{-2\beta }=L_h^{-\lfloor 2\beta \rfloor }\left(\sum _{i=1}^m{r}_i{(L_h-p_i{I}_{V_h})}^{-1}+k{I}_{V_h}\right).$$(14)

Here, ${\left\{r_i\right\}}_{i=1}^m,{\left\{p_i\right\}}_{i=1}^m$ and k are real numbers, and $I_{V_h}$ is the identity operator mapping the finite element space to itself. Let C be the mass matrix with elements ${\mathit{C}}_{i,j}={({\phi }_i,{\phi }_j)}_{L_2(\mathcal{D})}$, and let L be the matrix obtained by the bilinear form $a_L(·,·)$ induced by the differential operator L, which has elements ${\mathit{L}}_{i,j}={(\mathit{H}\nabla {\phi }_i,\nabla {\phi }_j)}_{L_2(\mathcal{D})}+{({\kappa }^2{\phi }_i,{\phi }_j)}_{L_2(\mathcal{D})}$. Then, we can use (14) to obtain the covariance matrix of w as (see Appendix C for a derivation): (15) $${\mathbf{\Sigma }}_{\mathit{w}}={({\mathit{L}}^{-1}\mathit{C})}^{\lfloor 2\beta \rfloor }\sum _{i=1}^m{r}_i{(\mathit{L}-p_i\mathit{C})}^{-1}+{\mathit{K}}_{\lfloor 2\beta \rfloor },$$(15) where ${\mathit{K}}_0=k\mathit{C}$ and ${\mathit{K}}_n=k{({\mathit{L}}^{-1}\mathit{C})}^{n-1}{\mathit{L}}^{-1}$ when $n\ge 1,n\in \mathbb{N}$. In the Matérn case, that is, when κ is a constant and H is an identity matrix, we simply have $\mathit{L}=\mathit{G}+{\kappa }^2\mathit{C}$, where G is the stiffness matrix with elements ${\mathit{G}}_{i,j}={(\nabla {\phi }_i,\nabla {\phi }_j)}_{L_2(\mathcal{D})}$.

Since we have the same degree for numerator and denominator in the rational approximation, we can use the BRASIL algorithm (Hofreither 2021) to compute the coefficients ${\left\{a_i\right\}}_{i=0}^m$ and ${\left\{b_i\right\}}_{i=0}^m$ in (12) and thus the coefficients ${\left\{r_i\right\}}_{i=0}^m,{\left\{p_i\right\}}_{i=0}^m$ and k in (14). Another option, commonly used in practice, is to use a “near best” rational approximation. One such option, which was used in Bolin and Kirchner (2020), and which is also implemented in the rSPDE package, is the Clenshaw–Lord Chebyshev–Padé algorithm (Baker and Graves-Morris 1996). See Appendix F for details about this algorithm. Also, observe that the interval $[{\lambda }_{n_h,h}^{-1},{\lambda }_{1,h}^{-1}]$ where one should compute the rational approximation may vary with the parameters κ and H, and that recomputing the coefficients ${\left\{a_i\right\}}_{i=0}^m$ and ${\left\{b_i\right\}}_{i=0}^m$ for different values of these parameters is not practical for implementations. To avoid this, recall from Assumption 2 that ${\kappa }_0^2$ is a lower bound for the eigenvalues of L in the case of Neumann boundary conditions and that ${\lambda }_1\le {\lambda }_{1,h}$ (see Proposition B.1 in Appendix B). We can, then, re-scale the operator L_h as $L_h/{\kappa }_0^2$ so that we can replace the interval $[{\lambda }_{n_h,h}^{-1},{\lambda }_{1,h}^{-1}]$ by $[\delta ,1]$, where, ideally, δ is chosen in such way that $\delta \le {\kappa }_0^2/{\lambda }_{n_h,h}$ for all considered mesh sizes h. In the rSPDE package, the choices δ = 0 and $\delta ={10}^{-(5+m)/2}$ are implemented. However, the difference in accuracy with respect to approximating the covariance function is negligible between these two choices.

For these options, we verified empirically that if $f_{\beta }(x)=x^{\left\{2\beta \right\}},\left\{2\beta \right\}=2\beta -\lfloor 2\beta \rfloor$, and ${\widehat{f}}_{\beta ,m}$ is the rational approximation of $f_{\beta }$ where the numerator and denominator have same degree m, then ${\widehat{f}}_{\beta ,m}=x^{\lfloor 2\beta \rfloor }{\sum }_{j=1}^m{r}_i{(x-p_i)}^{-1}+k,$ where ${\left\{p_i\right\}}_{i=1}^m$ are negative real numbers and ${\left\{r_i\right\}}_{i=1}^m$ and k are positive real numbers. This, together with the fact that the BRASIL algorithm is only implemented for rational approximations with numerator and denominator having the same degree, are the main reasons we chose to consider the numerator and denominator having the same degree m. Bolin and Kirchner (2020) instead considered a rational approximation where the numerator has degree m and the denominator has degree m + 1. However, with this choice the partial fractions would not yield a decomposition into positive-definite operators in our case.

Since ${\left\{p_i\right\}}_{i=1}^m$ are negative real numbers and ${\left\{r_i\right\}}_{i=1}^m$ and k are positive real numbers, we have that $r_i{({\mathit{L}}^{-1}\mathit{C})}^{\lfloor 2\beta \rfloor }{(\mathit{L}-p_i\mathit{C})}^{-1}$, for $i=1,\dots m$, and ${\mathit{K}}_{\lfloor 2\beta \rfloor }$ are valid covariance matrices. Thus, w can be expressed as a sum of m + 1 independent random vectors ${\mathit{x}}_i$ as in (3), where ${\mathit{Q}}_i$ is the precision matrix of ${\mathit{x}}_i$. By (15), we obtain that (16) $${\mathit{Q}}_i=\{\begin{array}{ccc}{r}_i^{-1}(\mathit{L}-p_i\mathit{C}){({\mathit{C}}^{-1}\mathit{L})}^{\lfloor 2\beta \rfloor }& & i=1,\dots ,m,\\ {\mathit{K}}_{\lfloor 2\beta \rfloor }^{-1}& & i=m+1.\end{array}$$(16)

Let $\mathit{X}={[{\mathit{x}}_1^{\top },\dots ,{\mathit{x}}_{m+1}^{\top }]}^{\top }$. Then, the precision matrix of X is the block diagonal matrix shown in (5). The final step in order to obtain a GMRF representation is to use the mass lumping technique as for the standard SPDE approach, see Appendix C.5 in Lindgren et al. (2011). Thus, the mass matrix C in (16) is replaced by a lumped mass matrix $\tilde{\mathit{C}}$, where $\tilde{\mathit{C}}$ is a diagonal matrix with ${\tilde{\mathit{C}}}_{\mathit{\text{ii}}}={\sum }_{j=1}^{n_h}{\mathit{C}}_{\mathit{\text{ij}}}$, for $i=1,\dots ,n_h$. With this adjustment, Q in (5) is sparse and we thus have obtained a GMRF representation.

6 Implementation and the rSPDE Package

The proposed covariance-based rational approximation method has been implemented in the R package rSPDE. In the following sections, we will use this package to illustrate the performance of the method. In this section, we give a brief introduction to the package and how it can be used in combination with R-INLA for computationally efficient Bayesian inference of latent Gaussian models involving the generalized Whittle–Matérn fields.

The usual workflow of fitting standard SPDE models in R-INLA can be divided into six steps. Namely, constructing the FEM mesh, defining SPDE model, creating a projector matrix, building the INLA stack, specifying the model formula, and finally calling the function inla to fit the model. Details about this can be found in Lindgren and Rue (2015). To fit a model with a fractional SPDE, this procedure remains the same. The only difference is that when defining the SPDE model, creating the projector matrix and building the index for INLA stack, we use functions from the rSPDE package. These functions are very similar to the corresponding R-INLA functions in terms of functionality. For example, a fractional SPDE model can be created with the command:

model <- rspde.matern(mesh = mesh) where mesh is a FEM mesh that can be obtained by inla.me sh.2d function from R-INLA. The default order of the rational approximation in this function is m = 2, which provides a good trade off between computational cost and accuracy, see Figure 1. As for the corresponding inla.spde2.matern function that can be used to define non-fractional SPDE models in INLA, one can also set priors for κ and τ in rspde.matern. Further, we can also define a prior for the smoothness parameter ν or specify ν so that a SPDE model with a fixed smoothness parameter can be generated. This feature can be used, for example, in the case that one already knows what ν is or wants to compare two different models with different ν, as we will do in Section 8.

Figure 1. Errors in $L_2(\mathcal{D}\times \mathcal{D})$-norm (top) and supremum norm ($L_{\infty }(\mathcal{D}\times \mathcal{D})$) (bottom) on $\mathcal{D}={[0,1]}^2$ for different practical ranges ρ for different values of ν. All methods use the same FEM mesh, with 100 equally spaced nodes in each direction.

The projector matrix $\overline{\mathit{A}}$ for a given mesh and observation locations loc is computed as

A <- rspde.make.A(mesh = mesh,

loc = loc)

As for the creation of the model, the default order of the rational approximation when creating the projector matrix is $m=2$, which can be changed by the user. The other arguments of the function are the same as those in the corresponding R-INLA function inla.spde.make.A. In the step of building the INLA stack, usually an index set is needed. The index can be computed with the function rspde.make.index, which replaces the R-INLA function inla.spde.make.index and has the same arguments. With these functions, the fractional models can be used as any other random effect in R-INLA. After fitting the model with the R-INLA function inla, posterior samples from a latent field and hyperparameters can be obtained by using inla.posterior.sample, and posterior distributions of the model parameters can be extracted via the rSPDE function rspde.result.

Besides the INLA-related functions, the rSPDE package also provides various utility functions. For example, once a fractional SPDE model, model, has been created with the rspde.matern function, one can simulate from it by calling simulate(model) to obtain a prior sample from a given choice of parameters, and the marginal log-likelihood from Section 2 can be computed by

l <- rSPDE.matern.loglike(model,

y, A, sigma.e)

Here, y is the observed data and sigma.e is the standard deviation of the measurement noise. In addition, if a model is fitted with this approach, then kriging and posterior sampling can be obtained by using the predict function. For further details and examples, we refer the reader to the vignettes at https://davidbolin.github.io/rSPDE.

Finally, rSPDE also provides an interface to the inlabru package (Bachl et al. 2019), which simplifies the construction of spatial models. This was used in the application in Section 8, where the entire code for defining and fitting the fractional model is

mesh <- inla.mesh.2d(loc = loc,

max.edge = c(0.5, 10),

cutoff = 0.35)

spde <- rspde.matern(mesh = mesh,

nu.upper.bound = 1)

res <- bru(z ∼ −1 + field(coordinates,

model = spde), data = data)

Here loc are the measurement locations and data is a data frame with the locations and observations.

7 Numerical Experiments

In this section, we compare the accuracy of the covariance–based rational approximation with the operator-based method from Bolin and Kirchner (2020), and with the “parsimonious” method from Lindgren et al. (2011). Since the latter method is implemented in R-INLA, we refer to it as the INLA approximation. We also note that the INLA method constructs a covariance-based Markov approximation (see also Bolin and Kirchner 2020, sec. 2), so it can be viewed as a 0th order covariance-based rational approximation.

For the comparison, we consider the SPDE model (2) with homogeneous zero Neumann boundary conditions on the unit square $\mathcal{D}={[0,1]}^2$, with τ chosen such that ${\sigma }^2$ in the Matérn covariance is one. The reason we consider the square domain, is that we have an explicit expression for the covariance function of the solution u. Indeed, we have, from Khristenko et al. (2019, eq. (2.13)), that the covariance function of u is given by (17) $$\begin{array}{c}{\varrho }_u^{\beta }\left(\mathit{x},\mathit{y}\right)=\sum _{\mathit{k}\in {\mathbb{Z}}^2}[\varrho (\left|\right|\mathit{x}+2\mathit{k}-\mathit{y}\left|\right|)\\ +\varrho (||(x_1+2k_1-y_1,x_2+2k_2+y_2)||)\\ +\varrho (||(x_1+2k_1+y_1,x_2+2k_2-y_2)||)\\ +\varrho (||\mathit{x}+2\mathit{k}+\mathit{y}||)],\end{array}$$(17) where $\left|\right|·\left|\right|$ is the Euclidean norm on ${\mathbb{R}}^2$ and $\varrho (·)$ is the Matérn covariance function in (1) with σ = 1 and $\nu =2\beta -1$. To compare the accuracy of the covariance approximations, we evaluate the true and approximate covariance functions on a regular mesh on ${[0,1]}^2$ with $N=100$ equally spaced nodes on each axis. We will compare these approximations with respect to the $L_2({[0,1]}^2\times {[0,1]}^2)$-norm and the supremum norm on ${[0,1]}^2\times {[0,1]}^2$. Appendix G shows how we approximate the errors in these two norms in detail.

For the operator-based and covariance-based rational approximations, we consider the orders of rational approximation as $m=1,2,3,4$. We choose smoothness parameters ranging from 0.1 to 3.1 with steps of size 0.05. Further, we test three possible values of κ. These values of κ, say ${\kappa }_1(\nu ),{\kappa }_2(\nu )$ and ${\kappa }_3(\nu )$ are chosen in such a way that the practical range $\rho =\sqrt{8\nu }/\kappa$ is fixed as $0.1,0.5$, and 1, respectively, for all values of ν. The resulting errors for the different methods are shown in Figure 1.

We begin by observing that for smoothness parameters $\nu =1,2$, or 3, there is no rational approximation and the errors come from the FEM approximation. With this in mind, one should note that for smaller range parameters most of the approximation error comes from the FEM approximation, thus, yielding a small difference of errors across the different methods. However, for larger ranges, such as, in this case, practical range equal to 1, the errors have different orders of magnitude as the order of the rational approximation increases, with the errors from the operator-based and covariance-based approximations of same rational approximation order having approximately the same order of magnitude. Furthermore, we can observe numerical instabilities of the operator-based approximations of order 3 and 4 as ν increases for both practical ranges 0.5 and 1, whereas the covariance-based method is stable for all orders of approximation.

In order to further illustrate the effect of the FEM error on the rational approximation of the covariance operator we repeated the analysis from above but with a coarser FEM mesh, consisting of 50 equally spaced nodes on each axis over the domain ${[0,1]}^2$. The results are shown in Figure G.1 in Appendix G. We now observe that for practical range 0.1, there is no visible difference between the covariance-based or operator-based rational approximations of orders 1 to 4, with a very small difference between the “parsimonious” INLA approximation and the remaining rational approximations. Further, for practical ranges 0.5 and 1, we hardly see any differences between the rational approximations of orders 2, 3, and 4. The only noticeable difference being that for large values of ν, the operator-based rational approximation becomes numerically unstable. On the other hand, it is noteworthy that for practical ranges 0.5 and 1, there is a significant difference (difference in orders of magnitude) between the rational approximations of order 0, 1 and the remaining orders.

To summarize, the results indicate that the covariance-based method generally has a similar accuracy as the operator-based method, which is higher than the accuracy provided by INLA’s method. The results also show that the covariance-based method is more numerically stable, especially for larger values of m, the order of the rational approximation.

It is important to remember that the INLA method only provides a fixed approximation, furthermore it only works in this case of stationary parameters, whereas the other methods are applicable also for nonstationary models and can be made arbitrarily precise by increasing the order m. As previously mentioned, the operator-based method is not suitable for inference in R-INLA, but the covariance-based method is. Thus, in conclusion, the covariance-based method provides a method that facilitates inference for stationary and nonstationary fractional SPDE-based models in R-INLA, which is also more accurate than the current INLA method for stationary models.

The numerical experiments in this section were implemented using the rSPDE package. All plots in this section, along with several more, for different choices of all the parameters involved, can be found in a shiny (Chang et al. 2021) app available at https://github.com/davidbolin/rSPDE. The results above were obtained by the Clenshaw–Lord Chebyshev–Padé algorithm with δ = 0 (see Section 5). The shiny app also contains the results by the BRASIL algorithm and the Clenshaw–Lord Chebyshev–Padé algorithm with $\delta ={10}^{-(5+m)/2}$. We also include results on likelihood errors in Appendix G.

8 Application

In this section, we illustrate the usage of the covariance-based rational approximation method through an application to a spatial dataset of precipitation observations. The dataset, available at https://www.image.ucar.edu/Data/precip_tapering/ contains annual precipitation anomalies observed by weather stations in the United States (standardized by the long-run mean and standard deviation for each station). We study the data from the year 1962, which contains observations from 7352 stations throughout the contiguous United States. We chose this dataset because it is simple enough to use a stationary model, which allows us to highlight the advantages of the fractional model without having to construct a complicated hierarchical model. Kaufman, Schervish, and Nychka (2008) also studied this data as an illustration for the covariance tapering method.

We model the data by (4) where $\mathit{ϵ}$ is independent Gaussian measurement noise with ${\mathit{Q}}_{\mathit{ϵ}}={\sigma }_{ϵ}^{-2}\mathbf{I}$ and u is a Whittle-Matérn field obtained as a solution to (2), where $\mathcal{D}$ is a bounded region (see Figure 2). The field is discretized using a finite element mesh that covers the contiguous United States with 9485 nodes. Figure 2 shows the mesh and the 7352 stations. Our interest is to compare the stationary SPDE models with either a fractional smoothness parameter ν (referred to as the fractional model) or a fixed parameter ν = 1 (referred to as the integer model) in terms of predictive power. In order to more easily interpret the parameters, we consider a parameterization of the Whittle–Matérn field in terms the standard deviation $\sigma =\sqrt{\Gamma (\nu )}/(\tau {\kappa }^{\nu }\sqrt{(4\pi )\Gamma (\nu +1)})$, the practical correlation range $\rho =\sqrt{8\nu }/\kappa$, and the smoothness ν. The prior distributions for the parameters are chosen as the default choices from the rSPDE package. That is, the priors of $\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}(\rho )$ and $\hspace{0.17em}\mathrm{\text{log}}\hspace{0.17em}(\sigma )$ are independent Gaussian distributions with variance 10 and the mean values are chosen based on size of the domain. Further, the prior of ν is a Beta distribution on the interval (0, 1) with mean 1/2 and variance 1/16. The choice of prior for ν is motivated by the fact that we do not believe that this should be a very smooth field. We also tested with Beta distributions on larger intervals and found that this did not affect the parameter estimates or the predictive performance of the model much.

Figure 2. the finite element mesh over the contiguous U.S. And the stations shown as dots.

We fit the models using R (version 4.2.1) and the rSPDE package (version 2.2.0) combined with inlabru (version 2.7.0) and R-INLA (version 23.02.17) running on a machine with an Intel i9-12900KF CPU, 64GB RAM and an Ubuntu operating system. The complete code for the analysis can be found in the supplementary materials. The total time for fitting the fractional and integer models are 38.4 sec and 15.4 sec, respectively.

The posterior distributions of the parameters of the Gaussian field and standard deviation of measurement noise for the two models are shown in Figure 3. One can note that the posterior mode of ν for the fractional model is around 0.52, which indicates that a fractional smoothness is needed. Compared to the fractional model, the integer model has a smaller σ and a larger ${\sigma }_{ϵ}$, indicating that the latent field explains less of the variability of the data. Finally, the practical correlation range of the integer model is substantially smaller than that of the fractional model, which likely is caused by the fact that a small range is needed to better explain the short range behavior of the data if the smoothness parameter is forced to be an integer.

Figure 3. Posterior distributions of σ, ρ, ν, and ${\sigma }_{ϵ}$.

To further compare the models, we perform two leave-group-out pseudo cross-validation studies (Liu and Rue 2022). In the first, for each station, we predict the value of the station based on all data except that from stations that are closer than a certain distance D (referred to as the distance of removed data). We then vary this distance and compute the accuracy of the predictions as functions of D. In the second, for each station, we instead remove the data from the k nearest stations and compute the accuracy of the predictions as functions of k. According to the screening effect (Stein 2002), in both cases the removed observations are the most informative. The quality of the prediction is measured in terms of Mean Squared Errors (MSE) and the negative Log-Score (LS) (Good 1952). Both metrics are negatively oriented, which means that a lower value indicates a better result. The results of the two cross-validation studies are shown in Figure 4. We see that the fractional model outperforms the integer model in both cases. For example, the fractional model with the distance of removed data being 400 km achieves the same levels of MSE and negative LS as the integer model with the distance of removed data being 300 km (indicated by dashed lines). Also, the fractional model with 125 removed data points achieves the same levels of MSE and negative LS as the integer model with 100 removed data points. We can note that the two models have similar performance when the distance of removed data and number of removed data are close to zero. This is expected due to the mean-squared continuity of the latent fields, combined with the fact that the models have nugget effects. This means that both models will have an MSE close to the variance of measurement error when the distance of removed data or number of nearest removed data are close to zero.

Figure 4. MSE and negative Log-Score as functions of distance (in km) of removed data (top) and number of removed data (bottom).

9 Discussion

We have introduced a new rational SPDE approach which provides stable and computationally efficient approximations for the covariance structure of generalized Whittle–Matérn Gaussian random fields with general smoothness $\beta >d/4$. We further derived an explicit rate of convergence of the method, which provides a theoretical justification for the approach. Compared to the rational SPDE approach of Bolin and Kirchner (2020), the main advantage is that we obtain a GMRF representation of the approximation. This allowed us to implement the method so that fractional SPDE models now can be estimated in R-INLA, where we in particular can estimate the smoothness parameter from data.

The current version of rSPDE has truncated log-normal and beta priors for the smoothness parameter ν as possible choices. A natural question for future research is how this prior should be chosen in a more systematic way. A potential way to do this is following the idea of penalized complexity priors (PC-priors) (Simpson et al. 2017). Fuglstad et al. (2019) derived PC-priors for κ and τ of the Whittle–Matérn fields assuming a fixed value of ν. We plan to extend that work by deriving PC-priors for all three parameters. Another potential area of future work is to extend the proposed method to spatio-temporal SPDE models as those proposed by Lindgren et al. (2020).

Supplementary Materials

Appendices: The seven appendices of the manuscript (Appendix.pdf).

R code and data: Code for replicating results of the application considered in Section 8 with data associated with the code (code_for_application.zip).

Acknowledgments

Our sincere thanks to Elias T. Krainski and Håvard Rue for their help with explaining some details of the internal structure of the R-INLA software and to the anonymous reviewers for insightful comments and suggestions on the article.

Disclosure Statement

The authors report there are no competing interests to declare.