The SPDE approach for Gaussian random fields with general smoothness

A popular approach for modeling and inference in spatial statistics is to represent Gaussian random fields as solutions to stochastic partial differential equations (SPDEs) $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 covariance models such as the exponential covariance on $\mathbb{R}^2$. We demonstrate how this restriction can be avoided by combining a finite element discretization in space with a rational approximation of the function $x^{-\beta}$ to approximate the solution $u$. For the resulting approximation, an explicit rate of strong convergence is derived and we show that the method has the same computational benefits as in the restricted case $2\beta\in\mathbb{N}$ when used for statistical inference and prediction. Several numerical experiments are performed to illustrate the accuracy of the method, and to show how it can be used for likelihood-based inference for all model parameters including $\beta$.