The Rational SPDE Approach for Gaussian Random Fields With General Smoothness

Abstract

A popular approach for modeling and inference in spatial statistics is to represent Gaussian random fields as solutions to stochastic partial differential equations (SPDEs) of the form $L^{\beta }u=\mathcal{W}$, where $\mathcal{W}$ is Gaussian white noise, L is a second-order differential operator, and $\beta >0$ is a parameter that determines the smoothness of u. However, this approach has been limited to the case $2\beta \in \mathbb{N}$, which excludes several important models and makes it necessary to keep β fixed during inference. We propose a new method, the rational SPDE approach, which in spatial dimension $d\in \mathbb{N}$ is applicable for any $\beta >d/4$, and thus remedies the mentioned limitation. The presented scheme combines a finite element discretization with a rational approximation of the function $x^{-\beta }$ to approximate u. For the resulting approximation, an explicit rate of convergence to u in mean-square sense is derived. Furthermore, we show that our method has the same computational benefits as in the restricted case $2\beta \in \mathbb{N}$. Several numerical experiments and a statistical application are used to illustrate the accuracy of the method, and to show that it facilitates likelihood-based inference for all model parameters including β.

Supplementary materials for this article are available online.

Keywords:

• Fractional operators
• Matérn covariances
• Nonstationary Gaussian fields
• Spatial statistics
• Stochastic partial differential equations

1 Introduction

One of the main challenges in spatial statistics is to handle large datasets. A reason for this is that the computational cost for likelihood evaluation and spatial prediction is in general cubic in the number N of observations of a Gaussian random field. A tremendous amount of research has been devoted to coping with this problem and various methods have been suggested (see Heaton et al. 2019, for a recent review).

A common approach is to define an approximation u_h of a Gaussian random field u on a spatial domain $\mathcal{D}$ via a basis expansion, (1.1) $$u_h(\mathbf{s})=\sum _{j=1}^n{u}_j {\phi }_j(\mathbf{s}),\hspace{1em}\mathbf{s}\in \mathcal{D},$$(1.1) where ${\phi }_j: \mathcal{D}\to \mathbb{R}$ are fixed basis functions and the vector $\mathbf{u}={(u_1,\dots ,u_n)}^{\top }\sim \mathrm{N}(0,{\mathbf{\Sigma }}_{\mathbf{u}})$ are stochastic weights. The computational effort can then be reduced by choosing $n\ll N$. However, methods based on such low-rank approximations tend to remove fine-scale variations of the process. For this reason, methods which instead exploit sparsity for reducing the computational cost have gained popularity in recent years. One can construct sparse approximations either of the covariance matrix of the measurements (Furrer, Genton, and Nychka 2006), or of the inverse of the covariance matrix (Datta et al. 2016). Alternatively, one can let the precision matrix ${\mathbf{\Sigma }}_{\mathbf{u}}^{-1}$ of the weights in (1.1) be sparse, as in the stochastic partial differential equation (SPDE) approach by Lindgren, Rue, and Lindström (2011), where usually $n\approx N$. To increase the accuracy further, several combinations of the methods mentioned above have been considered (e.g., Sang and Huang 2012) and multiresolution approximations of the process have been exploited (Nychka et al. 2015; Katzfuss 2017). However, theoretical error bounds have not been derived for most of these methods, which necessitates tuning these approximations for each specific model.

In this work we propose a new class of approximations, whose members we refer to as rational SPDE approximations or rational approximations for short. Our approach is similar to some of the above methods in the sense that an expansion (1.1) with compactly supported basis functions is exploited. The main novelty is that neither the covariance matrix ${\mathbf{\Sigma }}_{\mathbf{u}}$ nor the precision matrix ${\mathbf{\Sigma }}_{\mathbf{u}}^{-1}$ of the weights u is assumed to be sparse. The covariance matrix is instead a product ${\mathbf{\Sigma }}_{\mathbf{u}}=\mathbf{P}{\mathbf{Q}}^{-1}{\mathbf{P}}^{\top }$, where P and Q are sparse matrices and the sparsity pattern of P is a subset of that of Q. We show that the resulting approximation facilitates inference and prediction at the same computational cost as a Markov approximation with ${\mathbf{\Sigma }}_{\mathbf{u}}^{-1}=\mathbf{Q}$, and at a higher accuracy.

For the theoretical framework of our approach, we consider a Gaussian random field on a bounded domain $\mathcal{D}\subset {\mathbb{R}}^d$ which can be represented as the solution u to the SPDE (1.2) $$L^{\beta }u=\mathcal{W}\hspace{1em}\text{in} \mathcal{D},$$(1.2) where $\mathcal{W}$ is Gaussian white noise on $\mathcal{D}$, and $L^{\beta }$ is a fractional power of a second-order differential operator L which determines the covariance structure of u. Our rational approximations are based on two components: (i) a finite element method (FEM) with continuous and piecewise polynomial basis functions ${\left\{{\phi }_j\right\}}_{j=1}^n$, and (ii) a rational approximation of the function $x^{-\beta }$. We explain how to perform these two steps in practice to explicitly compute the matrices P and Q. Furthermore, we derive an upper bound for the strong mean-square error of the rational approximation. This bound provides an explicit rate of convergence in terms of the mesh size of the finite element discretization, which facilitates tuning the approximation without empirical tests for each specific model.Examples of random fields which can be expressed as solutions to SPDEs of the form (1.2) include approximations of Gaussian Matérn fields (Matérn 1960). Specifically, if $\mathcal{D}:={\mathbb{R}}^d$ a zero-mean Gaussian Matérn field can be viewed as a solution u to (1.3) $${({\kappa }^2-\Delta )}^{\beta } (\tau u)=\mathcal{W},$$(1.3) where Δ denotes the Laplacian (Whittle 1963). The constant parameters $\kappa ,\tau >0$ determine the practical correlation range and the variance of u. The exponent β defines the smoothness parameter ν of the Matérn covariance function via the relation $2\beta =\nu +d/2$ and, thus, the differentiability of the field. For applications, variance, range and differentiability typically are the most important properties of a Gaussian field. For this reason, the Matérn model is highly popular in spatial statistics and has become the standard choice for Gaussian process priors in machine learning (Rasmussen and Williams 2006). Since (1.3) is a special case of (1.2) we believe that the outcomes of this work will be of great relevance for many statistical applications, see also Section 7.

In contrast to covariance-based models, the SPDE approach additionally has the advantage that it allows for a number of generalizations of stationary Matérn fields including (i) nonstationary fields generated by nonstationary differential operators (e.g., Fuglstad et al. 2015), (ii) fields on more general domains such as the sphere (e.g., Lindgren, Rue, and Lindström 2011), and (iii) non-Gaussian models (e.g., Wallin and Bolin 2015).

Lindgren, Rue, and Lindström (2011) showed that, in the case that $2\beta \in \mathbb{N}$, one can construct accurate approximations of the form (1.1) for Gaussian Matérn fields on bounded domains $\mathcal{D}⊊{\mathbb{R}}^d$, such that ${\mathbf{\Sigma }}_{\mathbf{u}}^{-1}$ is sparse. To this end, (1.3) is considered on $\mathcal{D}$ and the differential operator ${\kappa }^2-\Delta$ is augmented with appropriate boundary conditions. The resulting SPDE is then solved approximately by means of a FEM. Due to the implementation in the R-INLA software, this approach has become widely used, see Bakka et al. (2018) for a comprehensive list of recent applications.

However, the restriction $2\beta \in \mathbb{N}$ implies a significant limitation for the flexibility of the method. In particular, it is therefore not directly applicable to the important special case of exponential covariance ($\nu =1/2$) on ${\mathbb{R}}^2$, where $\beta =3/4$. In addition, restricting the value of β complicates estimating the smoothness of the process from data. In fact, β typically is fixed when the method is used in practice, since identifying the value of $2\beta \in \mathbb{N}$ with the highest likelihood requires a separate estimation of all the other parameters in the model for each value of β. A common justification for fixing β is to argue that it is not practicable to estimate the smoothness of a random field from data. However, there are certainly applications for which it is feasible to estimate the smoothness. We provide an example of this in Section 7. Furthermore, having the correct smoothness of the model is particularly important for interpolation, and the fact that the Matérn model allows for estimating the smoothness from data was the main reason for why Stein (1999) recommended the model.

The rational SPDE approach presented in this work facilitates an estimation of β from data by providing an approximation of u which is computable for all fractional powers $\beta >d/4$ (i.e., $\nu >0$), where $d\in \mathbb{N}$ is the dimension of the spatial domain $\mathcal{D}\subset {\mathbb{R}}^d$. It thus enables to include this parameter in likelihood-based (or Bayesian) parameter estimation for both stationary and nonstationary models. Although the SPDE approach has been considered in the nonfractional case also for nonstationary models, Lindgren, Rue, and Lindström (2011) showed convergence of the approximation only for the stationary model (1.3). Our analysis in Section 3 closes this gap since we consider the general model (1.2) which covers the nonstationary case and several other previously proposed generalizations of the Matérn model.

The structure of this article is as follows: We briefly review existing methods for the SPDE approach in the fractional case in Section 2. In Section 3, the rational SPDE approximation is introduced and a result on its strong convergence is stated. The procedure of applying the rational SPDE approach to statistical inference is addressed in Section 4. Section 5 contains numerical experiments which illustrate the accuracy of the proposed method. The identifiability of the parameters in the Matérn SPDE model (1.3) is discussed in Section 6, where we derive necessary and sufficient conditions for equivalence of the induced Gaussian measures. In Section 7, we present an application to climate data, where we consider fractional and nonfractional models for both stationary and nonstationary covariances. We conclude with a discussion in Section 8. Finally, the supplementary materials contain four appendices providing details about (A) the finite element discretization, (B) the convergence analysis, (C) a comparison with the quadrature method by Bolin, Kirchner, and Kovács (2018a), and (D) the equivalence of Gaussian measures. The method developed in this work has been implemented in the R (R Core Team 2017) package rSPDE (Bolin 2019).

2 The SPDE Approach in the Fractional Case Until Now

A reason for why the approach by Lindgren, Rue, and Lindström (2011) only works for integer values of $2\beta$ is given by Rozanov (1977), who showed that a Gaussian random field on ${\mathbb{R}}^d$ is Markov if and only if its spectral density can be written as the reciprocal of a polynomial, $\tilde{S}(\mathbf{k})={(2\pi )}^{-d}{({\sum }_{j=0}^m{b}_j|\left|\mathbf{k}\right|{|}^{2j})}^{-1}$. Since the spectrum of a Gaussian Matérn field is (2.1) $$S(\mathbf{k})=\frac{1}{{(2\pi )}^d}\frac{1}{{({\kappa }^2+|\left|\mathbf{k}\right|{|}^2)}^{2\beta }},\hspace{1em}\mathbf{k}\in {\mathbb{R}}^d,$$(2.1) the precision matrix Q will not be sparse unless $2\beta \in \mathbb{N}$. For $2\beta \notin \mathbb{N}$, Lindgren, Rue, and Lindström (2011) suggested to compute a Markov approximation by choosing $m=\lceil 2\beta \rceil$ and selecting the coefficients $\mathbf{b}={(b_1,\dots ,b_m)}^{\top }$ so that the deviation between the spectral densities $${\int }_{{\mathbb{R}}^d}w(\mathbf{k}){(S(\mathbf{k})-\tilde{S}(\mathbf{k}))}^2 \mathrm{d}\mathbf{k}$$ is minimized. For this measure of deviation, w is some suitable weight function which should be chosen to get a good approximation of the covariance function. For the method to be useful in practice, the coefficients b_j should be given explicitly in terms of the parameters κ and ν. Because of this, Lindgren, Rue, and Lindström (2011) proposed a weight function that enables an analytical evaluation of the integral, $${\int }_{{\kappa }^2}^{\infty }[z^{2\beta }-\sum _{j=0}^m{b}_j{(z-{\kappa }^2)}^j{]}^2{z}^{-2m-1-\theta } \mathrm{d}z,$$ where $\theta >0$ is a tuning parameter. By differentiating this integral with respect to the parameters and setting the differentials equal to zero, a system of linear equations is obtained, which can be solved to find the coefficients b. The resulting approximation depends strongly on θ, and one could use numerical optimization to find a good value of θ for a specific value of β, or use the choice $\theta =2\beta -\lfloor 2\beta \rfloor$, which approximately minimizes the maximal distance between the covariance functions (Lindgren, Rue, and Lindström 2011). This method was used for the comparison in Heaton et al. (2019), and we will use it as a baseline method when analyzing the accuracy of the rational SPDE approximations in later sections.

Another Markov approximation based on the spectral density was proposed by Roininen et al. (2018). These Markov approximations may be sufficient in certain applications; however, any approach based on the spectral density or the covariance function is difficult to generalize to models on more general domains than ${\mathbb{R}}^d$, nonstationary models, or non-Gaussian models. Thus, such methods cannot be used if the full potential of the SPDE approach should be kept for fractional values of β.

There is a rich literature on methods for solving deterministic fractional PDEs (e.g., Gavrilyuk, Hackbusch, and Khoromskij 2004; Bonito and Pasciak 2015; Jin et al. 2015; Nochetto, Otárola, and Salgado 2015), and some of the methods that have been proposed could be used to compute approximations of the solution to the SPDE (1.3). However, any deterministic problem becomes more sophisticated when randomness is included. Even methods developed specifically for sampling solutions to SPDEs like (1.3) may be difficult to use directly for statistical applications, when likelihood evaluations, spatial predictions or posterior sampling are needed. For instance, it has been unclear if the sampling approach by Bolin, Kirchner, and Kovács (2018a), which is based on a quadrature approximation for an integral representation of the fractional inverse $L^{-\beta }$, could be used for statistical inference. In Appendix C, we show that it can be viewed as a (computationally less efficient) version of the rational SPDE approximations developed in this work. Consequently, the results in Section 4 on how to use the rational SPDE approach for inference apply also to that method. In Section 5.1, we compare the performance of the two methods in practice within the scope of a numerical experiment.

3 Rational Approximations for Fractional SPDEs

In this section, we propose an explicit scheme for approximating solutions to a class of SPDEs including (1.3). Specifically, in Sections 3.1 and 3.2, we introduce the fractional order equation of interest as well as its finite element discretization. In Section 3.3, we propose a nonfractional equation, whose solution after specification of certain coefficients approximates the random field of interest. For this approximation, we provide a rigorous error bound in Section 3.4. Finally, in Section 3.5, we address the computation of the coefficients in the rational approximation.

3.1 The Fractional Order Equation

With the objective of allowing for more general Gaussian random fields than the Matérn class, we consider the fractional order Equation (1(1.1) $$u_h(\mathbf{s})=\sum _{j=1}^n{u}_j {\phi }_j(\mathbf{s}),\hspace{1em}\mathbf{s}\in \mathcal{D},$$(1.1) .2), where $\mathcal{D}\subset {\mathbb{R}}^d$, $d\in \{1,2,3\}$, is an open, bounded, convex polytope, with closure $\overline{\mathcal{D}}$, and $\mathcal{W}$ is Gaussian white noise in $L_2(\mathcal{D})$. Here and below, $L_2(\mathcal{D})$ is the Lebesgue space of square-integrable real-valued functions, which is equipped with the inner product ${(w,v)}_{L_2(\mathcal{D})}:={\int }_{\mathcal{D}}w(\mathbf{s})v(\mathbf{s}) \mathrm{d}\mathbf{s}$. The Sobolev space of order $k\in \mathbb{N}$ is denoted by $H^k(\mathcal{D}):=\left\{w\in L_2(\mathcal{D}):D^{\gamma }w\in L_2(\mathcal{D})\forall \left|\gamma \right|\le k\right\}$ and $H_0^1(\mathcal{D})$ is the subspace of $H^1(\mathcal{D})$ containing functions with vanishing trace.

We assume that the operator $L: \mathit{D}(L)\subset L_2(\mathcal{D})\to L_2(\mathcal{D})$ is a linear second-order differential operator in divergence form, (3.1) $$Lu=-\nabla ·(\mathbf{H}\nabla u)+{\kappa }^{2}u,$$(3.1) whose domain of definition $\mathit{D}(L)$ depends on the choice of boundary conditions on $\partial \mathcal{D}$. Specifically, we impose homogeneous Dirichlet or Neumann boundary conditions and set $V=H_0^1(\mathcal{D})$ or $V=H^1(\mathcal{D})$, respectively. Furthermore, we let the functions H and κ in (3.1) satisfy the following assumptions:

1. $\mathbf{H}: \overline{\mathcal{D}}\to {\mathbb{R}}^{d\times d}$ is symmetric, Lipschitz continuous on the closure $\overline{\mathcal{D}}$, that is, $$\exists C_{\mathrm{Lip}}>0:\hspace{1em}|H_{ij}(\mathbf{s})-H_{ij}(\mathbf{s}\prime )|\le C_{\mathrm{Lip}}\left|\right|\mathbf{s}-\mathbf{s}\prime \left|\right|$$ $$\forall \mathbf{s},\mathbf{s}\prime \in \overline{\mathcal{D}},\hspace{1em}\forall i,j\in \{1,\dots ,d\},$$

and uniformly positive definite, that is, $$\exists C_0>0:\hspace{1em}\mathrm{ess}\underset{\mathbf{s}\in \mathcal{D}}{\inf }{\mathit{\xi }}^{\top }\mathbf{H}(\mathbf{s})\mathit{\xi }\ge C_0\left|\right|\mathit{\xi }|{|}^2\hspace{1em}\forall \mathit{\xi }\in {\mathbb{R}}^d;$$

1. $\kappa : \mathcal{D}\to \mathbb{R}$ is essentially bounded, $\kappa \in L_{\infty }(\mathcal{D})$. Furthermore, if Neumann boundary conditions are imposed, then $\mathrm{ess}\underset{\mathbf{s}\in \mathcal{D}}{\inf }\kappa (\mathbf{s})\ge {\kappa }_0>0$ holds.

If I.–II. are satisfied, the differential operator L in (3.1) induces a symmetric, continuous and coercive bilinear form a_L on V, (3.2) $$a_L(u,v):={(\mathbf{H}\nabla u,\nabla v)}_{L_2(\mathcal{D})}+{({\kappa }^{2}u,v)}_{L_2(\mathcal{D})},\hspace{1em}u,v\in V,$$(3.2) and its domain is given by $\mathit{D}(L)=H^2(\mathcal{D})\cap V$. Furthermore, Weyl’s law (see, e.g., Davies 1995, Theorem 6.3.1) shows that the eigenvalues ${\left\{{\lambda }_j\right\}}_{j\in \mathbb{N}}$ of the elliptic differential operator L in (3.1), in nondecreasing order, satisfy the spectral asymptotics (3.3) $${\lambda }_j\eqsim j^{2/d}\hspace{1em}\text{as} j\to \infty .$$(3.3)

Thus, existence and uniqueness of the solution u to (1.2) readily follow from Lemma 2.1 and Proposition 2.3 of Bolin, Kirchner, and Kovács (2018a). We formulate this as a proposition.

Proposition 3.1

. Let L be given by (3.1) where H and κ satisfy the assumptions I.–II. above and assume $\beta >d/4$. Then (1.2) has a unique solution u in $L_2(\Omega ;L_2(\mathcal{D}))$.

The assumptions I.–II. on the differential operator L are satisfied, for example, by the Matérn operator $L={\kappa }^2-\Delta$, in which case the condition $\beta >d/4$ on the fractional exponent in (1.2) corresponds to a positive smoothness parameter ν, that is, to a nondegenerate field. Moreover, Equation (1(1.1) $$u_h(\mathbf{s})=\sum _{j=1}^n{u}_j {\phi }_j(\mathbf{s}),\hspace{1em}\mathbf{s}\in \mathcal{D},$$(1.1) .2) as considered in our work includes several previously proposed nonfractional nonstationary models as special cases, such as the nonstationary Matérn models by Lindgren, Rue, and Lindström (2011), the models with locally varying anisotropy by Fuglstad et al. (2015), and the barrier models by Bakka et al. (2019). Thus, Proposition 3.1 shows existence and uniqueness of the fractional versions of all these models, which can be treated in practice by using the results of the following sections.

3.2 The Discrete Model

To discretize the problem, we assume that $V_h\subset V$ is a finite element space with continuous piecewise linear basis functions ${\left\{{\phi }_j\right\}}_{j=1}^{n_h}$ defined with respect to a triangulation ${\mathcal{T}}_h$ of the domain $\overline{\mathcal{D}}$ indexed by the mesh width $h:={\max }_{T\in {\mathcal{T}}_h}{h}_T$, where $h_T:=\mathrm{diam}(T)$ is the diameter of the element $T\in {\mathcal{T}}_h$. Furthermore, 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 supposed to be quasi-uniform, that is, there exist constants $C_1,C_2>0$ such that ${\rho }_T\ge C_1{h}_T$ and $h_T\ge C_{2}h$ for all $T\in {\mathcal{T}}_h$ and $h\in (0,1)$. Here, ${\rho }_T>0$ is radius of largest ball inscribed in $T\in {\mathcal{T}}_h$.

The discrete operator $L_h: V_h\to V_h$ is defined in terms of the bilinear form a_L in (3.2) via the relation ${(L_h{\varphi }_h,{\psi }_h)}_{L_2(\mathcal{D})}=a_L({\varphi }_h,{\psi }_h)$ which holds for all ${\varphi }_h,{\psi }_h\in V_h$. We then consider the following SPDE on the finite-dimensional state space V_h, (3.4) $$L_h^{\beta }{u}_h={\mathcal{W}}_h\hspace{1em}\text{in} \mathcal{D},$$(3.4) where ${\mathcal{W}}_h$ is Gaussian white noise in V_h, that is, ${\mathcal{W}}_h={\sum }_{j=1}^{n_h}{\xi }_j{e}_{j,h}$ for a basis ${\left\{e_{j,h}\right\}}_{j=1}^{n_h}$ of V_h which is orthonormal in $L_2(\mathcal{D})$ and ${\xi }_j\sim \mathrm{N}(0,1)$ iid for $j=1,\dots ,n_h$.

We note that the assumptions I.–II. from Section 3.1 on the functions H and κ combined with the convexity of $\mathcal{D}$ imply that the operator L in (3.1) is $H^2(\mathcal{D})$-regular, that is, for a right-hand side $f\in L_2(\mathcal{D})$, the weak solution $u\in V$ to Lu = f satisfies $u\in H^2(\mathcal{D})\cap V$, see, for example, Grisvard (2011, Theorem 3.2.1.2) for the case of Dirichlet boundary conditions. By combining this observation with the spectral asymptotics (3.3) we see that the assumptions in Lemmas 3.1 and 3.2 of Bolin, Kirchner, and Kovács (2018a) are satisfied (since then, in their notation, $r=s=q=2$ and $\alpha =2/d$) and we obtain an error estimate for the finite element approximation $u_h=L_h^{-\beta }{\mathcal{W}}_h$ in (3.4) for all $\beta \in (d/4,1)$. Furthermore, since their derivation requires only that $\beta >d/4$, we can formulate this result for all such values of β in the following proposition.

Proposition 3.2

. Suppose that $\beta >d/4$ and that L is given by (3.1) where H and κ satisfy the assumptions I.–II. from Section 3.1. Let u, u_h be the solutions to (1.2) and (3.4), respectively. Then, there exists a constant C > 0 such that, for sufficiently small h, $$\left|\right|u-u_h|{|}_{L_2(\Omega ;L_2(\mathcal{D}))}\le C{h}^{\text{min}\{2\beta -d/2, 2\}}.$$

3.3 The Rational Approximation

Proposition 3.2 shows that the mean-square error between u and u_h in $L_2(\mathcal{D})$ converges to zero as $h\to 0$. It remains to describe how an approximation of the random field u_h with values in the finite-dimensional state space V_h can be constructed.

For $\beta \in \mathbb{N}$ one can use, for example, the iterated finite element method presented in Appendix A to compute u_h in (3.4) directly. In the following, we construct approximations of u_h if $\beta \notin \mathbb{N}$ is a fractional exponent. For this purpose, we aim at finding a nonfractional equation (3.5) $$P_{\ell ,h}{u}_{h,m}^R=P_{r,h}{\mathcal{W}}_h\hspace{1em}\text{in} \mathcal{D},$$(3.5) such that $u_{h,m}^R$ is a good approximation of u_h, and where the operator $P_{j,h}:=p_j(L_h)$ is defined in terms of a polynomial p_j of degree $m_j\in {\mathbb{N}}_0$, for $j\in \{\ell ,r\}$. Since the so-defined operators $P_{\ell ,h}$, $P_{r,h}$ commute, this will lead to a nested SPDE model of the form (3.6) $$\begin{array}{c}{P}_{\ell ,h}{x}_{h,m}={\mathcal{W}}_h\text{in} \mathcal{D},\hfill \\ u_{h,m}^R=P_{r,h}{x}_{h,m}\hspace{1em}\text{in} \mathcal{D},\hfill \end{array}$$(3.6) which facilitates efficient computations, see Section 4 and Appendix A.

Comparing the initial Equation (1(1.1) $$u_h(\mathbf{s})=\sum _{j=1}^n{u}_j {\phi }_j(\mathbf{s}),\hspace{1em}\mathbf{s}\in \mathcal{D},$$(1.1) .2) with (3.7) $$P_{\ell }{u}_m^R=P_r\mathcal{W}\hspace{1em}\text{in} \mathcal{D},$$(3.7) where $P_j:=p_j(L),j\in \{\ell ,r\}$, motivates the choice $m_{\ell }-m_r\approx \beta$ to obtain a similar smoothness of $u_m^R={(P_r^{-1}{P}_{\ell })}^{-1}\mathcal{W}$ and $u=L^{-\beta }\mathcal{W}$ in (1.2). In practice, we first choose a degree $m\in \mathbb{N}$ and then set (3.8) $$\begin{array}{c}{m}_r:=m\hspace{1em}\text{and}\hspace{1em}{m}_{\ell }:=m+m_{\beta },\hspace{1em}\text{where}\hfill \\ m_{\beta }:=\max \{1,\lfloor \beta \rfloor \}.\hfill \end{array}$$(3.8)

In this case, the solution $u_m^R$ of (3.7) has the same smoothness as the solution v of the nonfractional equation $L^{\lfloor \beta \rfloor }v=\mathcal{W}$, if $\beta \ge 1$, and as v in $Lv=\mathcal{W}$, if $\beta <1$. Furthermore, for fixed h, the degree m controls the accuracy of the approximation $u_{h,m}^R$.

We now turn to the problem of defining the nonfractional operators $P_{\ell ,h}$ and $P_{r,h}$ in (3.5). To compute u_h in (3.4) directly, one would have to apply the discrete fractional inverse $L_h^{-\beta }$ to the noise term ${\mathcal{W}}_h$ on the right-hand side. Therefore, a first idea would be to approximate the function $x^{-\beta }$ on the spectrum of L_h by a rational function $\tilde{r}$ and to use $\tilde{r}(L_h){\mathcal{W}}_h$ as an approximation of u_h. This is, in essence, the approach proposed by Harizanov et al. (2018) to find optimal solvers for the problem ${\mathbf{L}}^{\beta }\mathbf{x}=\mathbf{f}$, where L is a sparse symmetric positive definite matrix. However, the spectra of L and of L_h as $h\to 0$ (considered as operators on $L_2(\mathcal{D})$) are unbounded and, thus, it would be necessary to normalize the spectrum of L_h for every h, since it is not feasible to construct the rational approximation $\tilde{r}$ on an unbounded interval. We aim at an approximation $L_h^{-\beta }\approx p_{\ell }{(L_h)}^{-1}{p}_r(L_h)$, where in practice the choice of $p_{\ell }$ and p_r can be made independent of L_h and h. Thus, we pursue another idea.

In contrast to the differential operator L in (3.1), its inverse $L^{-1}: L_2(\mathcal{D})\to L_2(\mathcal{D})$ is compact and, thus, the spectra of $L^{-1}$ and of $L_h^{-1}$ are bounded subsets of the intervals $J:=[0,{\lambda }_1^{-1}]$ and $J_h:=[{\lambda }_{n_h,h}^{-1},{\lambda }_{1,h}^{-1}]\subset J$, respectively, where ${\lambda }_{1,h},{\lambda }_{n_h,h}>0$ are the smallest and the largest eigenvalue of L_h. This motivates a rational approximation r of the function $f(x):=x^{\beta }$ on J and to deduce the nonfractional Equation (3(3.1) $$Lu=-\nabla ·(\mathbf{H}\nabla u)+{\kappa }^{2}u,$$(3.1) .5) from $u_{h,m}^R=r(L_h^{-1}){\mathcal{W}}_h$.

To achieve our envisaged choice (3.8) of different polynomial degrees $m_{\ell }$ and m_r, we decompose f via $f(x)=\widehat{f}(x)x^{m_{\beta }}$, where $\widehat{f}(x):=x^{\beta -m_{\beta }}$. We approximate $\widehat{f}\approx \widehat{r}:=\frac{q_1}{q_2}$ on J_h, where $q_1(x):={\sum }_{i=0}^m{c}_i{x}^i$ and $q_2(x):={\sum }_{j=0}^{m+1}{b}_j{x}^j$ are polynomials of degree m and m + 1, respectively, and use $r(x):=\widehat{r}(x)x^{m_{\beta }}$ as an approximation for f. This construction leads (after expanding the fraction with x^m) to a rational approximation $\frac{p_r}{p_{\ell }}$ of $x^{-\beta }$, (3.9) $$\begin{array}{c}{x}^{-\beta }=f(x^{-1})\approx \widehat{r}(x^{-1})x^{-m_{\beta }}=\frac{q_1(x^{-1})}{q_2(x^{-1})x^{m_{\beta }}}\hfill \\ =\frac{\sum _{i=0}^m{c}_i{x}^{m-i}}{\sum _{j=0}^{m+1}{b}_j{x}^{m+m_{\beta }-j}},\hfill \end{array}$$(3.9) where the polynomials $p_r(x):={\sum }_{i=0}^m{c}_i{x}^{m-i}$ and $p_{\ell }(x):={\sum }_{j=0}^{m+1}{b}_j{x}^{m+m_{\beta }-j}$ are of degree m and $m+m_{\beta }$, respectively, that is, (3.8) is satisfied.

The operators $P_{\ell ,h},P_{r,h}$ in (3.5) are defined accordingly, (3.10) $$\begin{array}{c}{P}_{\ell ,h}:=p_{\ell }(L_h)=\sum _{j=0}^{m+1}{b}_j{L}_h^{m+m_{\beta }-j},\hfill \\ P_{r,h}:=p_r(L_h)=\sum _{i=0}^m{c}_i{L}_h^{m-i}.\hfill \end{array}$$(3.10)

Their continuous counterparts in (3.7) are $P_{\ell }:=p_{\ell }(L)$ and $P_r:=p_r(L)$. We note that, for (3.8) to hold, any choice $m_2\in \{0,1,\dots ,m+m_{\beta }\}$ would have been permissible for the polynomial degree of q₂, if m is the degree of q₁. The reason for setting $m_2=m+1$ is that this is the maximal choice which is universally applicable for all values of $\beta >d/4$.

In the following we refer to $u_{h,m}^R$ in (3.5) with $P_{\ell ,h},P_{r,h}$ defined by (3.10) as the rational SPDE approximation of degree m. We emphasize that this approximation relies (besides the finite element discretization) only on the rational approximation of the function $\widehat{f}$. In particular, no information about the operator L except for a lower bound of the eigenvalues is needed. In the Matérn case, we have $L={\kappa }^2-\Delta$ (with certain boundary conditions) and an obvious lower bound of the eigenvalues is therefore given by ${\kappa }^2$.

3.4 An Error Bound for the Rational Approximation

In this subsection, we justify the approach proposed in Sections 3.2 and 3.3 by providing an upper bound for the strong mean-square error $\left|\right|u-u_{h,m}^R|{|}_{L_2(\Omega ;L_2(\mathcal{D}))}$. Here u and $u_{h,m}^R$ are the solutions of (1.2) and (3.5) and the rational approximation $u_{h,m}^R$ is constructed as described in Section 3.3, assuming that $\widehat{r}={\widehat{r}}_h$ is the $L_{\infty }$-best rational approximation of $\widehat{f}(x)=x^{\beta -m_{\beta }}$ on the interval J_h for each h. This means that ${\widehat{r}}_h$ minimizes the error in the supremum norm on J_h among all rational approximations of the chosen degrees in numerator and denominator. How such approximations can be computed is discussed in Section 3.5.

The theoretical analysis presented in Appendix B results in the following theorem, showing strong convergence of the rational approximation $u_{h,m}^R$ to the exact solution u.

Theorem 3.1.

Suppose that $\beta >d/4$ and that L is given by (3.1) where H and κ satisfy the assumptions I.–II. from Section 3.1. Let u, $u_{h,m}^R$ be the solutions to (1.2) and (3.5), respectively. Then, there is a constant C > 0, independent of h, m, such that, for sufficiently small h, $$\begin{array}{c}\left|\right|u-u_{h,m}^R|{|}_{L_2(\Omega ;L_2(\mathcal{D}))}\\ \le C(h^{\text{min}\{2\beta -d/2, 2\}}+1_{\beta \notin \mathbb{N}}{h}^{\text{min}\{2(\beta -1),0-d/2\}}{e}^{-2\pi \sqrt{|\beta -m_{\beta }|m}}).\end{array}$$

Remark 3.1.

To calibrate the accuracy of the rational approximation with the finite element error, one can choose $m\in \mathbb{N}$ such that $e^{-2\pi \sqrt{|\beta -m_{\beta }|m}}\propto h^{2\max \{\beta , 1\}}$. The strong rate of mean-square convergence is then $\min \{2\beta -d/2, 2\}$.

Remark 3.2

. If the functions H and κ in (3.1) are smooth, $\mathbf{H}\in C^{\infty }{(\overline{\mathcal{D}})}^{d\times d}$ and $\kappa \in C^{\infty }(\overline{\mathcal{D}})$ (as, e.g., in the Matérn case) and if the domain $\mathcal{D}$ has a smooth boundary, the higher-order strong mean-square convergence rate $\min \{2\beta -d/2, p+1\}$ can be proven for a finite element method with continuous basis functions which are piecewise polynomial of degree at most $p\in \mathbb{N}$. Thus, for $\beta >1$, finite elements with p > 1 may be meaningful.

3.5 Computing the Coefficients of the Rational Approximation

As explained in Section 3.3, the coefficients ${\left\{c_i\right\}}_{i=0}^m$ and ${\left\{b_j\right\}}_{j=0}^{m+1}$ needed for defining the operators $P_{\ell ,h},P_{r,h}$ in (3.10) are obtained from a rational approximation $\widehat{r}={\widehat{r}}_h$ of $\widehat{f}(x)=x^{\beta -m_{\beta }}$ on J_h. For each h, this approximation can, for example, be computed with the second Remez algorithm (Remez 1934), which generates the coefficients of the $L_{\infty }$-best approximation. The error analysis for the resulting approximation $u_{h,m}^R$ in (3.5) was performed in Section 3.4. Despite the theoretical benefit of generating the $L_{\infty }$-best approximation, the Remez algorithm is often unstable in computations and, therefore, we use a different method in our simulations. However, versions of the Remez scheme were used, for example, by Harizanov et al. (2018).

A simpler and computationally more stable way of choosing the rational approximation is, for instance, the Clenshaw–Lord Chebyshev–Padé algorithm (Baker and Graves-Morris 1996). To further improve the stability of the method, we will rescale the operator L so that its eigenvalues are bounded from below by one, which for the Matérn case corresponds to reformulating the SPDE (1.3) as ${(\text{Id}-{\kappa }^{-2}\Delta )}^{\beta }(\tilde{\tau }u)=\mathcal{W}$ and using $L=\text{Id}-{\kappa }^{-2}\Delta$, where $\text{Id}$ denotes the identity on $L_2(\mathcal{D})$ and $\tilde{\tau }:={\kappa }^{2\beta }\tau$.

To avoid computing a different rational approximation $\widehat{r}$ for each finite element mesh width h, in practice we compute the approximation $\widehat{r}$ only once on the interval $J_{*}:=[\delta ,1]$, where $\delta \in (0,1)$ should ideally be chosen such that $J_h\subset J_{*}$ for all considered mesh sizes h. For the numerical experiments later, we will use $\delta ={10}^{-(5+m)/2}$ when computing rational approximations of order m, which gives acceptable results for all values of β. As an example, the coefficients computed with the Clenshaw–Lord Chebyshev–Padé algorithm on $J_{*}$ for the case of exponential covariance on ${\mathbb{R}}^2$ are shown in Table 1.

Table 1 Coefficients of the rational approximation for $\beta =3/4$ (exponential covariance on ${\mathbb{R}}^2$) for m = 1, 2, 3, normalized so that c_m = 1.

m	b0	c0	b1	c1	b2	c2	b3	c3	b4
1	1.69e–2	7.69e–2	8.06e–1	1	2.57e–1
2	8.08e–4	5.30e–3	1.98e–1	4.05e–1	1.07	1	1.41e–1
3	3.72e–5	3.27e–4	3.03e–2	8.57e–2	6.84e–1	1.00	1.28	1	9.17e–2

4 Computational Aspects of the Rational Approximation

In the nonfractional case, the sparsity of the precision matrix for the weights u in (1.1) facilitates fast computation of samples, likelihoods, and other quantities of interest for statistical inference. The purpose of this section is to show that the rational SPDE approximation proposed in Section 3 preserves these good computational properties.

The representation (3.6) shows that $u_{h,m}^R$ can be seen as a Markov random field $x_{h,m}$, transformed by the operator $P_{r,h}$. Solving this latent model as explained in Appendix A, yields an approximation of the form (1.1), where ${\mathbf{\Sigma }}_{\mathbf{u}}={\mathbf{P}}_r{\mathbf{Q}}^{-1}{\mathbf{P}}_r^{\top }$. Here ${\mathbf{P}}_{\ell },{\mathbf{P}}_r\in {\mathbb{R}}^{n_h\times n_h}$ correspond to the discrete operators $P_{\ell ,h}$ and $P_{r,h}$ in (3.10), respectively. The matrix $\mathbf{Q}:={\mathbf{P}}_{\ell }^{\top }{\mathbf{C}}^{-1}{\mathbf{P}}_{\ell }$ is sparse if the mass matrix C with respect to the finite element basis ${\left\{{\phi }_j\right\}}_{j=1}^{n_h}$ is replaced by the diagonal lumped mass matrix $\tilde{\mathbf{C}}$, see Appendix A. By defining $\mathbf{x}\sim \mathsf{N}(0,{\mathbf{Q}}^{-1})$, we have $\mathbf{u}={\mathbf{P}}_r\mathbf{x}$, which is a transformed Gaussian Markov random field (GMRF). Choosing x as a latent variable instead of u thus enables us to use all computational methods, which are available for GMRFs (see Rue and Held 2005), also for the rational SPDE approximation.

As an illustration, we consider the following hierarchical model, with a latent field u which is a rational approximation of (1.2), (4.1) $$\begin{array}{c}{y}_i=u({\mathbf{s}}_i)+{\epsilon }_i,i=1,\dots ,N,\hfill \\ P_{\ell }u=P_r\mathcal{W}\text{in} \mathcal{D},\hfill \end{array}$$(4.1) where u is observed under iid Gaussian measurement noise ${\epsilon }_i\sim \mathrm{N}(0,{\sigma }^2)$. Given that one can treat this case, one can easily adapt the method to be used for inference in combination with MCMC or INLA (Rue, Martino, and Chopin 2009) for models with more sophisticated likelihoods.

Defining the matrix A with elements $A_{ij}={\phi }_j({\mathbf{s}}_i)$ and the vector $\mathbf{y}={(y_1,\dots ,y_N)}^{\top }$ gives us the discretized model (4.2) $$\begin{array}{c}\mathbf{y}|\mathbf{x}\sim \mathsf{N}(\mathbf{A}{\mathbf{P}}_r\mathbf{x},{\sigma }^2\mathbf{I}),\hfill \\ \mathbf{x}\sim \mathsf{N}(0,{\mathbf{Q}}^{-1}).\hfill \end{array}$$(4.2)

In this way, the problem has been reduced to a standard latent GMRF model and a sparse Cholesky factorization of Q can be used for sampling x from $\mathsf{N}(0,{\mathbf{Q}}^{-1})$ as well as to evaluate its log-density $\text{log}\hspace{0.17em}{\pi }_x(\mathbf{x})$. Samples of u can then be obtained from samples of x via $\mathbf{u}={\mathbf{P}}_r\mathbf{x}$. For evaluating the log-density of u, $\text{log}\hspace{0.17em}{\pi }_u(\mathbf{u})$, the relation $\text{log}\hspace{0.17em}{\pi }_u(\mathbf{u})=\hspace{0.17em}\text{log}\hspace{0.17em}{\pi }_x({\mathbf{P}}_r^{-1}\mathbf{u})$ can be exploited. Furthermore, the posterior distribution of x is given by $\mathbf{x}|\mathbf{y}\sim \mathrm{N}({\mathbf{\mu }}_{\mathbf{x}|\mathbf{y}},{\mathbf{Q}}_{\mathbf{x}|\mathbf{y}}^{-1})$, where ${\mathbf{\mu }}_{\mathbf{x}|\mathbf{y}}={\sigma }^{-2}{\mathbf{Q}}_{\mathbf{x}|\mathbf{y}}^{-1}{\mathbf{P}}_r^{\top }{\mathbf{A}}^{\top }\mathbf{y}$ and ${\mathbf{Q}}_{\mathbf{x}|\mathbf{y}}=\mathbf{Q}+{\sigma }^{-2}{\mathbf{P}}_r^{\top }{\mathbf{A}}^{\top }\mathbf{A}{\mathbf{P}}_r$ is a sparse matrix. Thus, simulations from the distribution of $\mathbf{x}|\mathbf{y}$, and evaluations of the corresponding log-density $\text{log}\hspace{0.17em}{\pi }_{x|y}(\mathbf{x})$, can be performed efficiently via a sparse Cholesky factorization of ${\mathbf{Q}}_{\mathbf{x}|\mathbf{y}}$. Finally, the marginal data log-likelihood is proportional to $$\begin{array}{c}\hspace{0.17em}\text{log}\hspace{0.17em}\left|{\mathbf{P}}_{\ell }\right|-\frac{1}{2}\text{log}\hspace{0.17em}\left|{\mathbf{Q}}_{\mathbf{x}|\mathbf{y}}\right|-N\hspace{0.17em}\text{log}\hspace{0.17em}\sigma \\ -\frac{1}{2}\left({\mathbf{\mu }}_{\mathbf{x}|\mathbf{y}}^{\top }\mathbf{Q} {\mathbf{\mu }}_{\mathbf{x}|\mathbf{y}}+{\sigma }^{-2}{\Vert \mathbf{y}-\mathbf{A}{\mathbf{P}}_r{\mathbf{\mu }}_{\mathbf{x}|\mathbf{y}}\Vert }^2\right).\end{array}$$

We therefore conclude that all computations needed for statistical inference can be facilitated by sparse Cholesky factorizations of ${\mathbf{P}}_{\ell }$ and ${\mathbf{Q}}_{\mathbf{x}|\mathbf{y}}$.

Remark 4.1

. From the specific form of the matrices ${\mathbf{P}}_{\ell }$ and ${\mathbf{P}}_r$ addressed in Appendix A, we can infer that the number of nonzero elements in ${\mathbf{Q}}_{\mathbf{x}|\mathbf{y}}$ for a rational SPDE approximation of degree m will be the same as the number of nonzero elements in ${\mathbf{Q}}_{\mathbf{x}|\mathbf{y}}$ for the standard (nonfractional) SPDE approach with $\beta =m+m_{\beta }$. Thus, also the computational cost will be comparable for these two cases.

Remark 4.2

. The matrix ${\mathbf{Q}}_{\mathbf{x}|\mathbf{y}}$ can be ill-conditioned for m > 1 if a FEM approximation with piecewise linear basis functions is used. The numerical stability for large values of m can likely be improved by increasing the polynomial degree of the FEM basis functions, see also Remark 3.2.

5 Numerical Experiments

5.1 The Matérn Covariance on ${\mathbb{R}}^2$

As a first test, we investigate the performance of the rational SPDE approach for Gaussian Matérn fields, without including the finite element discretization in space.

The spectral density S of the solution to (1.3) on ${\mathbb{R}}^2$ is given by (2.1), whereas the spectral density for the non-discretized rational SPDE approximation $u_m^R$ in (3.7) is (5.1) $$S_R(\mathbf{k})\propto {\kappa }^{4\beta }{\left(\frac{\sum _{i=1}^m{c}_i{(1+{\kappa }^{-2}|\left|\mathbf{k}\right|{|}^2)}^{m-i}}{\sum _{j=1}^{m+1}{b}_j{(1+{\kappa }^{-2}|\left|\mathbf{k}\right|{|}^2)}^{m+m_{\beta }-j}}\right)}^2.$$(5.1)

We compute the coefficients as described in Section 3.5. To this end, we apply an implementation of the Clenshaw–Lord Chebyshev–Padé algorithm provided by the Matlab package Chebfun (Driscoll, Hale, and Trefethen 2014). By performing a partial fraction decomposition of (5.1), expanding the square, transforming to polar coordinates, and using the equality $${\int }_0^{\infty }\frac{\omega J_0(\omega h)}{({\omega }^2+a^2)({\omega }^2+b^2)} \mathrm{d}\omega =\frac{1}{(b^2-a^2)}(K_0(ah)-K_0(bh)),$$ we are able to compute the corresponding covariance function $C_R(h)$ analytically. Here, J₀ is a Bessel function of the first kind and K₀ is a modified Bessel function of the second kind. To measure the accuracy of the approximation, we compare $C_R(h)$ to the true Matérn covariance function C(h) for different values of ν, where $\kappa =\sqrt{8\nu }$ is chosen such that the practical correlation range $r=\sqrt{8\nu }/\kappa$ equals one in all cases.

To put the accuracy of the rational approximation in context, the Markov approximation by Lindgren, Rue, and Lindström (2011) and the quadrature method by Bolin, Kirchner, and Kovács (2018a) are also shown. For the quadrature method, K = 12 quadrature nodes are used, which results in an approximation with the same computational cost as a rational approximation of degree m = 11, see Appendix C. Figure 1 shows the normalized error in the L₂-norm and the error with respect to $L_{\infty }$-norm for different values of ν, both with respect to the interval $[0,2]$ of length twice the practical correlation range, that is, $${\left(\frac{{\int }_0^2(C(h)-C_a(h){)}^2 \mathrm{d}h}{{\int }_0^{2}C{(h)}^2 \mathrm{d}h}\right)}^{1/2}\hspace{1em}\text{and}\hspace{1em}\underset{h\in [0,2]}{\sup }|C(h)-C_a(h)|.$$

Fig. 1 The L₂- and $L_{\infty }$-errors of the covariance functions for different values of ν for the different approximation methods. When ν = 1, all methods are exact.

Here, C_a is the covariance function obtained by the respective approximation method.

Already for m = 3, the rational approximation performs better than both the Markov approximation and the quadrature approximation for all values of ν. It also decreases the error for the case of an exponential covariance by several orders of magnitude.

All methods are exact when ν = 1, since this is the nonfractional case. The Markov and rational methods show errors decreasing to zero as ν = 1, whereas the error of the quadrature method has a singularity at ν = 1. The performance of the quadrature method can be improved (although not the behavior near ν = 1) by increasing the number of quadrature nodes, see Appendix C. This is reasonable if the method is needed only for sampling from the model, but implementing this method for statistical applications, which require kriging or likelihood evaluations, is not feasible since the computational costs then are comparable to the standard SPDE approach with $\beta =K$.

Finally, it should be noted that the Markov method also is exact at ν = 2 ($\beta =1.5$) since the spectrum of the process then is the reciprocal of a polynomial. The rational and quadrature methods cannot exploit this fact, since these approximations are based on the corresponding differential operator instead of the spectral density. This is the prize that has to be paid to formulate a method which works not only for the stationary Matérn fields but also for nonstationary and non-Gaussian models.

5.2 Computational Cost and the Finite Element Error

From the study in the previous subsection, we infer that the rational SPDE approach performs well for Matérn fields with arbitrary smoothness. However, as for the standard SPDE approach, we need to discretize the problem to be able to use the method in practice, for example, for inference. This induces an additional error source, which means that one should balance the two errors by choosing the degree m of the rational approximation appropriately with respect to the FEM error. A calibration based on the theoretical results has been suggested in Remark 3.1. In this section, we address this issue in practice and investigate the computational cost of the rational SPDE approximation.

As a test case, we compute approximations of a Gaussian Matérn field with unit variance and practical correlation range r = 0.1 on the unit square in ${\mathbb{R}}^2$. We assume homogeneous Neumann boundary conditions for the Matérn operator ${\kappa }^2-\Delta$ in (1.3). For the discretization, we use a FEM with a nodal basis of continuous piecewise linear functions with respect to a mesh induced by a Delaunay triangulation of a regular lattice on the domain, with a total of n_h nodes. We consider three different meshes with $n_h={57}^2,{85}^2,{115}^2$, which corresponds to the mesh sizes $h\approx r/4,r/6,r/8$.

To measure the accuracy, we compute the covariances between the midpoint of the domain ${\tilde{\mathbf{s}}}_{*}$ and all other nodes in the lattice ${\left\{{\tilde{\mathbf{s}}}_j\right\}}_{j=1}^{n_h}$ for the Matérn field and the rational SPDE approximations and calculate the error similarly to the L₂-error in Section 5.1, $${\left(\frac{\sum _{j=1}^{n_h}(C(\left|\right|{\tilde{\mathbf{s}}}_{*}-{\tilde{\mathbf{s}}}_j\left|\right|)-{\Sigma }_{j,*}^{\mathbf{u}}{)}^2}{\sum _{j=1}^{n_h}C{(||{\tilde{\mathbf{s}}}_{*}-{\tilde{\mathbf{s}}}_j||)}^2}\right)}^{1/2},$$ where ${\mathbf{\Sigma }}^{\mathbf{u}}={\mathbf{P}}_r{\mathbf{P}}_{\ell }^{-1}\mathbf{C}{\mathbf{P}}_{\ell }^{-\top }{\mathbf{P}}_r^{\top }$ is the covariance matrix of u, see Appendix A. As a consequence of imposing boundary conditions, the error of the covariance is larger close to the boundary of the domain. However, we compare this error to the error of the nonfractional SPDE approach, which has the same boundary effects. As measures of the computational cost, we consider the time it takes to sample u and to evaluate $\text{log}\hspace{0.17em}\left|{\mathbf{Q}}_{\mathbf{x}|\mathbf{y}}\right|$ for the model (4.2) with σ = 1, when y is a vector of noisy observations of the latent field at 1000 locations, drawn at random in the domain (a similar computation time is needed to evaluate ${\mathbf{\mu }}_{\mathbf{x}|\mathbf{y}}$).

The results for rational SPDE approximations of different degrees for the case $\beta =3/4$ (exponential covariance) are shown in Table 2. Furthermore, we perform the same experiment when the standard (nonfractional) SPDE approach is used for $\beta =2,3,4$. As previously mentioned in Remark 4.1, the computational cost of the rational SPDE approximation of degree m should be comparable to the standard SPDE approach with $\beta =m+1$. Table 2 validates this claim. One can also note that the errors of the rational SPDE approximations are similar to those of the standard SPDE approach, and that the reduction in error when increasing from m = 2 to m = 3 is small for all cases, indicating that the error induced by the rational approximation is small compared to the FEM error, even for a low degree m. This is also the reason for why, in particular in the pre-asymptotic region, one can in practice choose the degree m smaller than the value suggested in Remark 3.1, which gives $m\approx 6,7,8$ for $\beta =3/4$ and the three considered finite element meshes.

Table 2 Covariance errors ($\times 100$) and computing times in seconds ($\times 100$) for sampling from the rational SPDE approximation u (with $\beta =3/4$) and, in parentheses, for evaluating $\text{log}\hspace{0.17em}\left|{\mathbf{Q}}_{\mathbf{x}|\mathbf{y}}\right|$.

		Rational SPDE approximation			Standard SPDE approach
n		m = 1	m = 2	m = 3	β = 2	β = 3	β = 4
57^2	Error	1.849	1.339	1.415	2.259	2.173	2.147
	Time	1.5 (3.2)	1.8 (5.2)	2.7 (8.7)	1.7 (2.6)	1.7 (3.9)	2.2 (6.3)
85^2	Error	1.720	0.757	0.807	0.953	0.928	0.921
	Time	3.1 (8.4)	5.0 (14)	7.6 (25)	3.0 (8.2)	5.8 (13)	7.9 (22)
115^2	Error	1.559	0.526	0.501	0.509	0.498	0.494
	Time	7.6 (22)	11 (34)	18 (57)	6.3 (18)	11 (35)	18 (53)

NOTE: These values are also given for the standard SPDE approach with $\beta =2,3,4$.

6 Likelihood-Based Inference of Matérn Parameters

The computationally efficient evaluation of the likelihood of the rational SPDE approximation facilitates likelihood-based inference for all parameters of the Matérn model, including ν which until now had to be fixed when using the SPDE approach. In this section, we first discuss the identifiability of the model parameters and then investigate the accuracy of this approach within the scope of a simulation study.

6.1 Parameter Identifiability

A common reason for fixing the smoothness in Gaussian Matérn models is the result by Zhang (2004) which shows that all three Matérn parameters cannot be estimated consistently under infill asymptotics. More precisely, for a fixed smoothness parameter ν, one cannot estimate both the variance of the field, ${\varphi }^2$, and the scale parameter, κ, consistently. However, the quantity ${\varphi }^2{\kappa }^{2\nu }$ can be estimated consistently. The derivation of this result relies on the equivalence of Gaussian measures corresponding to Matérn fields (Zhang 2004, Theorem 2). The following theorem provides the analogous result for the Gaussian measures induced by the class of random fields specified via (1.3) on a bounded domain.

Theorem 6.1

. Let $\mathcal{D}\subset {\mathbb{R}}^d,d\in \{1,2,3\}$, be bounded, open, and connected. For $i\in \{1,2\}$, let ${\beta }_i>d/4$, ${\kappa }_i,{\tau }_i>0$, and consider the Gaussian measure ${\mu }_i:=\mathsf{N}(m_i,{\mathcal{Q}}_i^{-1})$ on $L_2(\mathcal{D})$ with mean $m_i:=0$ and precision operator ${\mathcal{Q}}_i:={\tau }_i^2{L}_i^{2{\beta }_i}$, where, for $i\in \{1,2\}$, the operators $L_i:={\kappa }_i^2-\Delta$ are augmented with the same homogeneous Neumann or Dirichlet boundary conditions. Then, μ₁ and μ₂ are equivalent if and only if ${\beta }_1={\beta }_2$ and ${\tau }_1={\tau }_2$.

The proof can be found in Appendix D. Note that, for $\mathcal{D}:={\mathbb{R}}^d$, the parameter τ is related to the variance of the Gaussian random field via ${\varphi }^2=\Gamma (\nu ){({\tau }^2\Gamma (2\beta ){(4\pi )}^{d/2}{\kappa }^{2\nu })}^{-1}$. Thus, ${\tau }^{-2}\propto {\varphi }^2{\kappa }^{2\nu }$, which means that Theorem 6.1 is in accordance with the result by Zhang (2004). Since the Gaussian measures induced by the operators $L_1=\tau {({\kappa }_1+\Delta )}^{\beta }$ and $L_2=\tau {({\kappa }_2+\Delta )}^{\beta }$ are equivalent, we will not be able to consistently estimate κ under infill asymptotics. Yet, Theorem 6.1 suggests that it is possible to estimate τ and β consistently. In fact, with Theorem 6.1 available, it is straightforward to show that τ can be estimated consistently for a fixed ν by exploiting the same arguments as in the proof of (Zhang 2004, Theorem 3). However, it is beyond the scope of this article to show that both ν and τ can be estimated consistently which would also extend the results by Zhang (2004).

6.2 Simulation Study

To numerically investigate the accuracy of likelihood-based parameter estimation using the rational SPDE approach, we again assume homogeneous Neumann boundary conditions for the Matérn operator in (1.3) and consider the standard latent model (4.1) from Section 4. We take the unit square as the domain of interest, set ${\sigma }^2=0.1,\nu =0.5$ and choose κ and τ so that the latent field has variance ${\varphi }^2=1$ and practical correlation range r = 0.2. For the FEM, we take a mesh based on a regular lattice on the domain, extended by twice the correlation range in each direction to reduce boundary effects, yielding a mesh with approximately 3500 nodes.

As a first test case, we use simulated data from the discretized model. We simulate 50 replicates of the latent field, each with corresponding noisy observations at 1000 measurement locations drawn at random in the domain. This results in a total of 50,000 observations, which we use to estimate the parameters of the model. We draw initial values for the parameters at random and then numerically optimize the likelihood of the model with the function fminunc in Matlab. This procedure is repeated 100 times, each time with a new simulated dataset.

As a second test case, we repeat the simulation study, but this time we simulate the data from a Gaussian Matérn field with an exponential covariance function instead of from the discretized model. For the estimation, we compute the rational SPDE approximation for the same finite element mesh as in the first test case. To investigate the effect of the mesh resolution on the parameter estimates, we also estimate the parameters using a uniformly refined mesh with twice as many nodes. The average computation time for evaluating the likelihood is approximately 0.16 sec for the coarse mesh and 0.4 sec for the fine mesh. This computation time is affine with respect to the number of replicates, and with only one replicate it is 0.09 sec for the coarse mesh and 0.2 sec for the fine mesh.

The results of the parameter estimation can be seen in Table 3, where the true parameter values are shown together with the mean and standard deviations of the 100 estimates for each case. Notably, we are able to estimate all parameters accurately in the first case. For the second case, the finite element discretization seems to induce a small bias, especially for the nugget estimate (${\sigma }^2$) that depends on the resolution of the mesh. The bias in the nugget estimate is not surprising since the increased nugget compensates for the FEM error. The bias could be decreased by choosing the mesh more carefully, also taking the measurement locations into account. In practice, however, this bias will not be of great importance, since the optimal nugget for the discretized model should be used.

Table 3 Results of the parameter estimation.

		Rational samples	Matérn samples
	Truth	Estimate	Coarse mesh	Fine mesh
κ	10	10.026 (0.5661)	10.966 (1.8060)	10.864 (0.4414)
${\varphi }^2$	1.0	1.0014 (0.0228)	1.1089 (0.6155)	0.9743 (0.0210)
${\sigma }^2$	0.1	0.1001 (0.0009)	0.3016 (0.0036)	0.2320 (0.0044)
ν	0.5	0.5011 (0.0168)	0.5554 (0.0991)	0.5462 (0.0138)

NOTE: For each parameter estimate, the mean of 100 different estimates is shown, with the corresponding standard deviation in parentheses.

It should be noted that there are several other methods for decreasing the computational cost of likelihood-based inference for stationary Matérn models. The major advantage of the rational SPDE approach is that it is directly applicable to more complicated nonstationary models, which we will use in the next section when analyzing real data.

7 Application

In this section, we illustrate for the example of a climate reanalysis dataset how the rational SPDE approach can be used for spatial modeling.

Climate reanalysis data are generated by combining a climate model with observations to obtain a description of the recent climate. We use reanalysis data generated with the Experimental Climate Prediction Center Regional Spectral Model (ECPC-RSM) which was originally prepared for the North American Regional Climate Change Assessment Program (NARCCAP) by means of NCEP/DOE Reanalysis (Mearns et al. 2009, 2014). As variable we consider average summer precipitation over the conterminous United States for a 26-year period from 1979 to 2004. The average value for each grid cell and year is computed as the average of the corresponding daily values for the days in June, July, and August. To obtain data which can be modeled by a Gaussian distribution, we follow Genton and Kleiber (2015) and transform the data by taking the cube root. We then subtract the mean over the 26 years from each grid cell so that we can assume that the data has zero mean and focus on the correlation structure of the residuals. The resulting residuals for the year 1979 are shown in Figure 2.

Fig. 2 Average summer precipitation residuals (in cm) for 1979 and the FEM mesh.

The 4112 observed residuals for each year are modeled as independent realizations of a zero-mean Gaussian random field with a nugget effect. That is, the measurement Y_ij at spatial location ${\mathbf{s}}_i$ for year j is modeled as $Y_{ij}=u_j({\mathbf{s}}_i)+{\epsilon }_{ij}$, where ${\epsilon }_{ij}\sim \mathsf{N}(0,{\sigma }^2)$ are independent, and ${\{u_j(\mathbf{s})\}}_j$ are independent realizations of a zero-mean Gaussian random field $u(\mathbf{s})$. The analysis of Genton and Kleiber (2015) revealed that an exponential covariance model is suitable for a subset of this dataset. Because of this, a natural first choice is to use a stationary Matérn model (1.3), either with $\beta =0.75$ (exponential covariance) or with a general β which we estimate from the data. However, since we have data for a larger spatial region than Genton and Kleiber (2015), one would suspect that a nonstationary model for $u(\mathbf{s})$ might be needed. The standard nonstationary model for the SPDE approach, as first suggested by Lindgren, Rue, and Lindström (2011) and used in many applications since then, is (7.1) $${(\kappa {(\mathbf{s})}^2-\Delta )}^{\beta } (\tau (\mathbf{s})u(\mathbf{s}))=\mathcal{W}(\mathbf{s}),\hspace{1em}\mathbf{s}\in \mathcal{D},$$(7.1) where β = 1 is fixed. Until now, it has not been possible to use the model (7.1) with fractional smoothness. Therefore, our main question is now: What is more important for this data—the fractional smoothness β or the nonstationary parameters? We thus consider four different SPDE models for $u(\mathbf{s})$. Two of them are nonfractional models, where β = 1 is fixed, and for the other two (fractional) models, we estimate the fractional order β jointly with the other parameters from the data. For both cases, we consider stationary Matérn and nonstationary models, where the latter are formulated via (7.1) with $$\text{log}\hspace{0.17em}\kappa (\mathbf{s})={\kappa }_0+{\kappa }_{\mathrm{a}}{\psi }_{\mathrm{a}}(\mathbf{s})+\sum _{i,j=1}^2\sum _{k,\ell =1}^2{\kappa }_{ij}^{k\ell } {\psi }_i^k({\tilde{s}}_1) {\psi }_j^{\ell }({\tilde{s}}_2),$$ and the same model is used for $\tau (\mathbf{s})$. Here, ${\psi }_j^1(\tilde{s}):=\hspace{0.17em}\text{sin}\hspace{0.17em}(j\pi \tilde{s})$, ${\psi }_j^2(\tilde{s}):=\hspace{0.17em}\text{cos}\hspace{0.17em}(j\pi \tilde{s}),{\psi }_{\mathrm{a}}(\mathbf{s})$ is the altitude at location s, and $\tilde{\mathbf{s}}=({\tilde{s}}_1,{\tilde{s}}_2)$ denotes the spatial coordinate after rescaling so that the observational domain is mapped to the unit square. Thus, $\text{log}\hspace{0.17em}\kappa (\mathbf{s})$ and $\text{log}\hspace{0.17em}\tau (\mathbf{s})$ are modeled by the altitude covariate and 16 additional Fourier basis functions to capture large-scale trends in the parameters. The altitude covariate and the eight Fourier basis functions $\left\{{\psi }_1^k({\tilde{s}}_1){\psi }_j^{\ell }({\tilde{s}}_2):j,k,\ell =1,2\right\}$ are shown in Figure 3.

Fig. 3 Nine basis functions modeling the parameters for the nonstationary models.

We discretize each model with respect to the finite element mesh shown in Figure 2, assuming homogeneous Neumann boundary conditions. The mesh has 5021 nodes and was computed using R-INLA (Lindgren and Rue 2015). For the fractional models, we set m = 1 in the rational approximation and, for each model, the model parameters are estimated by numerical optimization of the log-likelihood as described in Section 4.

The log-likelihood values for the four models can be seen in Table 4. The parameter estimates for the stationary nonfractional ($\beta =\nu =1$) model are $\kappa =0.67,\tau =5.44$, and $\sigma =0.014$, which implies a standard deviation $\varphi =0.077$ and a practical range $\rho =4.21$. The estimates for the fractional model are $\kappa =0.20,\tau =10.58,\sigma =0.012$, and $\beta =0.72$, corresponding to $\varphi =0.081,\rho =9.21$, and a smoothness parameter $\nu =0.44$. We note that the fractional model has a longer correlation range. This is likely to be caused by the nonfractional model underestimating the range ρ to compensate for the wrong local behavior of the covariance function induced by the smoothness parameter ν = 1.

Table 4 Model-dependent results for (i) the log-likelihood, (ii) the pseudo cross-validation scores (RMSE, CRPS, LS, each $\times 100$) averaged over ten replicates, and (iii) the computational time for one evaluation of the likelihood averaged over 100 computations.

	Log-likelihood	RMSE	CRPS	LS	Time
Stationary β = 1	219773	4.206	2.295	177.2	0.125
Stationary fractional	220255	4.167	2.274	178.2	0.412
Nonstationary β = 1	225969	4.194	2.266	182.1	0.121
Nonstationary fractional	226095	4.170	2.254	182.4	0.416

Figure 4 shows the estimated marginal standard deviation $\varphi (\mathbf{s})$ for the two nonstationary models (computed using the estimates of the parameters for $\kappa (\mathbf{s})$ and $\tau (\mathbf{s})$) and 0.7 contours of the correlation function for selected locations in the domain. The estimate of β for the nonstationary fractional model is 0.723. Also for the nonstationary models, we observe a slightly longer practical correlation range $\rho (\mathbf{s})$ for the fractional model.

Fig. 4 Estimated marginal standard deviations (top row) and contours of 0.7 correlation of the correlation function for selected locations (bottom row), for the fractional (left column) and β = 1 (right column) models.

To investigate the predictive accuracy of the models, a pseudo cross-validation study is performed. We choose 10% of the spatial observation locations at random, and use the corresponding observations for each year to predict the values at the remaining locations. The accuracy of the four models is measured by the root mean square error (RMSE), the average continuous ranked probability score (CRPS), and the average log-score (LS). This procedure is repeated ten times, where in each iteration new locations are chosen at random to base the predictions on. The average scores for the ten iterations are shown in Table 4. Recall that low RMSE and CRPS values, respectively, high LS values correspond to good scores.

We observe that the predictive performance of the nonstationary nonfractional (β = 1) model is similar to the stationary fractional model in terms of CRPS, and actually worse in terms of RMSE. This clearly indicates that the data should be analyzed by a fractional model. Although the nonstationary fractional model has a better performance in terms of CRPS and LS than the stationary fractional model, the difference is quite small given that the nonstationary model has 38 parameters, compared to 4 for the stationary model. Thus, the fractional smoothness seems to be the most important aspect for this data. The fact that the rational SPDE approach allows us to make these comparisons and to verify the smoothness parameter, for stationary and nonstationary models, is one of its most important features.

8 Discussion

We have introduced the rational SPDE approach providing a new type of computationally efficient approximations for a class of Gaussian random fields. These are based on an extension of the SPDE approach by Lindgren, Rue, and Lindström (2011) to models with general second-order differential operators of arbitrary order $\beta >d/4$. For these approximations, explicit rates of strong convergence have been derived and we have shown how to calibrate the degree of the rational approximation with the mesh size of the FEM to achieve these rates. The results can also be combined with the results in Bolin et al. (2018b) to obtain explicit rates of weak convergence (convergence of functionals of the random field).

Our approach can, for example, be used to approximate stationary Matérn fields with general smoothness, and it is also directly applicable to more complicated nonstationary models, where the covariance function may be unknown. A general fractional order of the differential operator opens up for new applications of the SPDE approach, such as to Gaussian fields with exponential covariances on ${\mathbb{R}}^2$. For the Matérn model and its extensions, it furthermore facilitates likelihood-based (or Bayesian) inference of all model parameters. The specific structure of the approximation then in turn enables a combination with INLA or MCMC in situations where the Gaussian model is a part of a more complicated non-Gaussian hierarchical model.

We have illustrated the rational SPDE approach for stationary and nonstationary Matérn models. A topic for future research is to apply the approach to other random field models in statistics which are difficult to approximate by GMRFs, such as to models with long-range dependence (Lilly et al. 2017) based on the fractional Brownian motion. Another topic for future research is to modify the fractional SPDE approach by replacing the FEM basis by a multiresolution basis and to compare this approach to other multiresolution approaches such as (Katzfuss 2017). Finally, it is also of interest to extend the method to non-Gaussian versions of the SPDE-based Matérn models (Wallin and Bolin 2015), since the Markov approximation considered by Wallin and Bolin (2015) is only computable under the restrictive requirement $\beta \in \mathbb{N}$.

Supplementary Materials

Appendices: The four appendices of the manuscript. (supplementary.pdf)

R code: Code for replicating the results of the application considered in Section 7. (code.zip)

Acknowledgments

The authors thank Stig Larsson, Finn Lindgren, the editor, and the anonymous reviewer for valuable comments on the article, which led to an improvement of the results and the presentation.

Additional information

Funding

This work has been supported in part by the Swedish Research Council under grant no. 2016-04187 and the Knut and Alice Wallenberg Foundation (KAW 20012.0067).