GAUSSIAN RANDOM FIELDS: WITH AND WITHOUT COVARIANCES

N. H. BINGHAM and Tasmin L. SYMONS

To Mikhailo Iosifovich Yadrenko, on his 90th birthday

Abstract
We begin with isotropic Gaussian random fields, and show how the Bochner-
Godement theorem gives a natural way to describe their covariance structure.
We continue with a study of Matérn processes on Euclidean space, spheres,
manifolds and graphs, using Bessel potentials and stochastic partial
differential equations (SPDEs). We then turn from this continuous setting to
approximating discrete settings, Gaussian Markov random fields (GMRFs), and
the computational advantages they bring in handling large data sets, by
exploiting the sparseness properties of the relevant precision (concentration)
matrices.

Key words Gaussian random field, covariance function, Bessel potential,
stochastic partial differential equation, Matérn process, Gaussian Markov
random field, precision matrix, sparseness, numerical linear algebra

0\. Introduction

In this survey, we study Gaussian random fields – Gaussian processes
parametrized by some space $\Omega$, which will usually be Euclidean space
${\mathbb{R}}^{d}$ of dimension $d$ or the sphere ${\mathbb{S}}^{d}$ of
dimension $d$ as a Riemannian manifold, embedded in ${\mathbb{R}}^{d+1}$; we
normalise to radius 1 for convenience. Gaussian processes are specified by
their mean function (which we shall take to be zero for simplicity), and
covariance function. So specifying Gaussian processes is the same as
specifying covariance functions, equivalently, positive-definite functions.
All this is in a continuous setting, where we can use our most powerful
weapon, calculus.
Data, however, and the computation needed to handle it, are discrete. The
tension between these two inspired the title of the 2012 paper [SimLR1]: ‘In
order to make spatial statistics computationally feasible, we need to forget
about the covariance function’ (of course, the authors do not dismiss
covariances for theory, only for computational statistics with large data
sets). Regarding Euclidean space and spheres, compare their later comment
[SimLP2]: ‘the great tragedy of spatial statistics is that the Earth turned
out not to be flat’.
We begin in §1 with a study of covariance functions in symmetric spaces (this
class, aptly named, includes both Euclidean space and spheres). Our main tool
is the powerful and under-utilised Bochner-Godement theorem. We continue in §2
with Bessel potentials (again, powerful and under-utilised in this area),
stochastic partial differential equations (SPDEs), and Matérn processes,
ubiquitous in the field.
Quoting the above authors again: ‘In contrast to traditional statistical
modelling, practical problems in spatial statistics are, by and large,
computational in nature’. We turn in §3 to Gaussian Markov random fields
(GMRFs), which arise [LinRL1] via graphs obtained by triangulation of the
space $\Omega$. Dependence occurs only between neighbouring vertices of the
graph. This means that the precision (or concentration) matrix (the inverse of
the covariance matrix) is sparse: most of its elements are zero. The
computational burden is then carried by sparse numerical linear algebra.
For further background on random fields, see [Mal], [MarP], [Yad].

1\. Covariances; the Bochner-Godement theorem

1.1. The Bochner and Bochner-Schoenberg theorems
In Euclidean space, the positive definite functions (normalised to be 1 at the
origin) are the characteristic functions, by Bochner’s theorem of 1933. These
are also the covariance functions [Kal1, 50]: for a process $\xi$ parametrised
by $T$, ${\rho}_{ij}=cov(\xi(\tau_{i}),\xi(\tau_{j}))$,

$\sum_{ij}a_{i}a_{j}\rho_{ij}=\sum
a_{i}a_{j}cov(\xi(\tau_{i}),\xi(\tau_{j}))=var(\sum_{i}a_{i}\xi(\tau_{i}))\geq
0,$

for all $a_{i}\in\mathbb{R},\ \tau_{i}\in T$. This result was extended from
processes on Euclidean space to isotropic processes on spheres by Bochner and
Schoenberg in 1940-42 (see [BinS4] for references). Here one needs the
ultraspherical (or Gegenbauer) polynomials $P_{n}^{\lambda}(u)$ of index
$\lambda$ ($u\in[-1,1]$, $\lambda=\frac{1}{2}(d-1)$); these are classical
orthogonal polynomials on the interval $[-1,1]$; see e.g. Szegő [Sze, §4.7].
The Bochner-Schoenberg theorem gives the general isotropic covariance on
${\mathbb{S}}^{d}$ to within scale by a convex combination (or mixture) of
Gegengauer polynomials:

$c\sum_{0}^{\infty}a_{n}P_{n}^{\lambda}(u),\ \ a_{n}\geq 0,\ \
\sum_{0}^{\infty}a_{n}=1,\ \ u=d({\bf x},{\bf y}),\ \ {\bf
x,y}\in{\mathbb{S}}^{d},\ \ \lambda=\frac{1}{2}(d-1)$ $None$

(here and below, isotropic means that the covariance between the values at two
points depends only on their geodesic distance $u=d({\bf x},{\bf y})$ as
here).

1.2. Symmetric spaces
One looks for a general framework to include both these classical results, and
this is given by symmetric spaces. This field, which stems from Elie Cartan in
the 1920s, belongs to differential and Riemannian geometry and Lie theory
rather than to probability. Below we summarise briefly what we need; for
background and detail, we refer to the standard works by Helgason [Hel1,2,3,4]
and Wolf [Wol1, Ch. 11]. Readers without a geometric background may prefer to
think of what are for us the prime examples: spheres, lines and half-lines,
and their products. For the relevant harmonic analysis (extending the Fourier
transform to the spherical transform), see §1.3, §1.4 and §1.8 below.
A symmetric space is a Riemannian manifold $M$ with curvature tensor invariant
under parallel translation. Equivalently, a geodesic symmetry exists: this
fixes some base point $o$ (the origin in ${\mathbb{R}}^{d}$, the North Pole in
${\mathbb{S}}^{d}$), and reverses the direction of the geodesics through $o$.
This gives an involutive isomorphism, mapping a point into its ‘mirror image’
on the geodesic through it and $o$. Then $M$ is a homogeneous space $G/K$: the
coset space of $G$, a closed subgroup of the isometry group of $M$ (containing
the transvections, [Wol2, 11.1A]) and $K$, the isotropy subgroup of $G$ fixing
$o$.

1.3. The spherical transform and spherical dual
To proceed, we need a version of the Fourier transform on symmetric spaces.
This is provided by the spherical transform, involving spherical measures and
spherical functions. A spherical measure $m$ is a $K$-bi-invariant
multiplicative linear functional on $C_{c}(K\backslash G/K)$ (continuous
functions of compact support on the double-coset space). The spherical
functions (originally, zonal spherical functions) are the continuous functions
$\omega:G\to\mathbb{C}$ such that the measure

$m_{\omega}(f):=\int_{G}f(g)\omega(g^{-1})d\mu_{G}(g)$

is spherical. The map

$f\mapsto\hat{f}(\omega):=m_{\omega}(f)=\int_{G}f(g)\omega(g^{-1})d\mu_{G}(g)$

is the spherical transform. The positive definite spherical functions $\phi$
on $(G,K)$ are in bijection with the irreducible unitary representations $\pi$
of $G$ with a $K$-fixed unit vector $u$ by

$\phi(g)=\langle u,\pi(g)u\rangle.$

These form the spherical dual, $\Lambda$.
We confine ourselves to the rank-one case [Hel1, V.6]. These fall into three
classes, the compact and non-compact (dual to each other) and the Euclidean
(self-dual) [Hel1, V.2]. These are the two-point homogeneous spaces, the
spaces of constant curvature $\kappa$: spheres (curvature $\kappa>0$),
Euclidean ($\kappa=0$) and hyperbolic spaces ($\kappa<0$) respectively [Hel1,
IX.5, X.3], [Wol2].
In the compact rank-one case, the $\pi$ and $\phi$ above, and $\Lambda$, are
in bijection with the Cartan-Weyl dominant weights [Wol1, §6.3], and with a
subset of $\mathbb{R}$, which we may again write as $\Lambda$, by the Cartan-
Helgason theorem [Hel3, V.1.1, 534-538], [Wol1, 11.4B].

1.4. The Bochner-Godement theorem
That Bochner’s theorem may be extended beyond ${\mathbb{R}}^{d}$ and
${\mathbb{S}}^{d}$ to symmetric spaces goes back to Godement in 1957, giving
the Bochner-Godement theorem. In its modern formulation, this very useful and
powerful result is as follows.

Theorem BG (Bochner-Godement theorem). The general positive definite function
$\psi$ on a symmetric space $M$ is given (to within scale $c$) by a mixture of
positive-definite spherical functions ${\phi}({\lambda})$ over the spherical
dual $\Lambda$ by a probability measure $\mu$:

$\psi(m)=c\int_{\Lambda}{\phi}_{\lambda}(m)\ d\mu(\lambda)\quad(m\in
M):\qquad\psi=c\int_{\Lambda}{\phi}_{\lambda}\ d\mu(\lambda).$ $None$

For background and details, see e.g. Wolf [Wol1, Th. 9.3.4], van Dijk [Dij,
Prop. 8.3.3] (there as the ‘Bochner-Godement-Schwartz theorem’). An early use
in a probability context (Gaussian processes on compact symmetric spaces) is
Askey and Bingham [AskB] in 1976.
Bochner’s theorem is the special case $G={\mathbb{R}}^{d}$, $K=\\{0\\}$ (here
the spherical functions are the characters, the complex exponentials
$e^{ix.}=[t\mapsto e^{ixt}]$). The Bochner-Schoenberg theorem is the case

${\mathbb{S}}^{d}=SO(d+1)/SO(d)\leftrightarrow(SO(d+1),SO(d)).$

Here the spherical functions are the Gegenbauer polynomials above.

1.5. Products of symmetric spaces; spatio-temporal and geotemporal processes
By the Schur product theorem (1911; [Sch], [HorJ, §7.5]), products of positive
definite functions are positive definite: the class ${\cal{P}}(M)$ of positive
definite functions is closed under pointwise products (which correspond to
tensor products of representations). For product spaces, if
$f_{i}\in{\cal{P}}(M_{i})$, then $f_{1}f_{2}\in{\cal{P}}(M_{1}\times M_{2})$.
The spherical dual of the product is the product of the spherical duals:

$M=M_{1}\times
M_{2},\qquad(G,K)=(G_{1},K_{1})\times(G_{2},K_{2}),\qquad\Lambda={\Lambda}_{1}\times{\Lambda}_{2}.$

The products of interest here are those in which the first factor is spatial
(Euclidean or spherical) and the second is time. We refer to the first case as
spatio-temporal and the second as geotemporal. Both are of great importance,
in both theory and applications. We recall Whittle’s advice, that one should
think of a spatial process as the equilibnrium process of a spatio-temporal
one [Whi3, Preface].

1.6. Gelfand pairs
For $G$ a locally compact group (Lie group will suffice here) with $K$ a
compact subgroup, $(G,K)$ is called a Gelfand pair if the convolution algebra
of $K$-bi-invariant continuous measures on $G$ with compact support is
commutative. For symmetric spaces $M=(G,K)$, $(G,K)$ is a Gelfand pair ([Hel1,
X, Th. 2.9, 4.1], [Hel3, IV, Th. 3,1]). We use the language of symmetric
spaces rather than Gelfand pairs, but Gelfand pairs (which can be discrete or
continuous) are very interesting; see van Dijk [Dij] for a monograph
treatment. Their relevance here is shown by there being a form of the Bochner-
Godement theorem for them ([Dij, Prop. 8.3.3] as noted above).

1.7. Sphere cross line; the Berg-Porcu theorem
From 1.5 above, one can produce geotemporal covariances by taking an isotropic
covariance on the sphere (from the Bochner-Schoenberg theorem), a stationary
covariance on the line (by Bochner’s theorem) and taking their product. This,
though a valid covariance on sphere cross line, is separable (one might as
well deal with the two separately). The interesting and useful covariances
here are inseparable. The question was raised by Mijatović and the authors
[BinMS, §4.4] in 2016 of finding the general (isotropic, stationary)
covariance on sphere cross line, and answered by Berg and Porcu [BerP] in
2017.
We state the Berg-Porcu theorem below, and illustrate the power and usefulness
of the Bochner-Godement theorem by using it to give a very short and simple
proof (indeed, two proofs).

Theorem BP (Berg-Porcu Theorem). The class of isotropic stationary covariances
on sphere cross line coincides (to within scale) with the class of mixtures of
products of Gegenbauer polynomials $P_{n}^{\lambda}(u)$ and characteristic
functions ${\phi}_{n}(t)$ on the line: with $u$, $\lambda$ as before, they are
given by

$c\sum_{0}^{\infty}a_{n}P_{n}^{\lambda}(u){\phi}_{n}(t),\quad c>0,\quad
a_{n}\geq 0,\quad\sum_{0}^{\infty}a_{n}=1.$ $None$

Proof. In the notation above, with the first factor the sphere
${\mathbb{S}}^{d}$ the spherical dual ${\Lambda}_{1}$ is the set of Gegenbauer
polynomials $P_{n}^{\lambda}$; with the second factor the line,
${\Lambda}_{2}$ is the set of characters, which can also be identified with
the line:

$t\leftrightarrow e^{it.}=(x\mapsto e^{ixt}).$

We can now proceed using the language of either measure theory or probability
theory.
1\. Measure theory. We use disintegration of measures in $(BG)$ (Fubini’s
theorem extended beyond product measures: see e.g. [Kal1, Th. 6.4], Bogachev
[Bog, §10.6]). Here,
${\Lambda}_{1}={\mathbb{N}}_{0}:={\mathbb{N}}\cup\\{0\\}$,
${\Lambda}_{2}=\mathbb{R}$. Integrate the probability measure $\mu$ on
$\Lambda$ over the second (time) variable for fixed $n$, then sum over $n$.
This disintegrates $\mu$ into a sequence of probability measures ${\mu}_{n}$
on the line, and a probability measure $(a_{n})$ on ${\mathbb{N}}_{0}$. Now
integrating $e^{itx}$ over ${\mu}_{n}(dx)$ gives its characteristic function
${\phi}_{n}(t)$, the second factor in $(BP)$. The remaining integration is a
sum over $n$ of the product of this and the first factor $P_{n}^{\lambda}(u)$
with weight $a_{n}$, giving the result.
2\. Probability theory. Take $\lambda=({\lambda}_{1},{\lambda}_{2})$ in $(BG)$
as a random variable with law $\mu$. Condition on its first coordinate, and
use the Conditional Mean Formula [Wil2, 390] (a special case of the tower
property [Wil1, 9.7i] or chain rule [Kal1, 105]). $\square$

Products of spheres. There is a similar result due to Guella, Menegatto and
Peron in 2016 [GueMP]. Its proof is immediate from the Bochner-Godement
theorem. See [BinMS] for details, [BinS3,4] for background.

1.8. Transforms: commutative and non-commutative settings
The prototypical transform for us is the Fourier transform. There are two
major extensions to commutative settings: Pontryagin duality on locally
compact groups (1934), and the Gelfand transform on commutative Banach
algebras in the 1940s [Rud, Ch. 11]. The spherical transform above derives
from Gelfand in 1950 [Gel] (see its review by Godement [God]). The
commutativity in a Gelfand pair captures what is needed to extend to a
genuinely non-commutative setting.

2\. Bessel potentials, SPDEs and Matérn processes

2.1. Riesz potentials
The two most lasting contributions of Marcel Riesz (1886-1969, younger brother
of Frederick Riez) are Riesz means (typical means) in summability theory, and
his work of 1937-38 on fractional potentials [Rie]. As potential theory is
intimately linked to the Laplacian, it is natural that one can treat Riesz
potentials in terms of fractional Laplacians (an important type of
pseudodifferential operators, [Tay]). See the now classic treatment in Stein
[Ste, V]. Now the Laplacian $\Delta$ has non-positive eigenvalues (as
$D^{2}e^{i\omega.}=-{\omega}^{2}e^{i\omega.}$), so it is preferable to deal
with $-\Delta$, which has non-negative eigenvalues, before taking fractional
powers. The operators

$I_{\alpha}:=(-\Delta)^{-\alpha/2}$

are called Riesz potentials, and in $d$ dimensions ${\mathbb{R}}^{d}$ are
well-behaved for $\alpha\in(0,d)$. Then for suitable functions $f$,

$\displaystyle(-\Delta)^{-\alpha/2}(f)$ $\displaystyle=$
$\displaystyle\frac{\Gamma(\frac{d}{2}-\frac{\alpha}{2})}{{\pi}^{d/2}2^{\alpha}\Gamma(\alpha/2)}\int_{{\mathbb{R}}^{d}}|x-y|^{-d+\alpha}f(y)dy\quad(0<\alpha<d)$
$\displaystyle=$ $\displaystyle\int_{{\mathbb{R}}^{d}}k(x,y)f(y)dy,$

say, where $k(.,.)$ is called the kernel of the operator $I_{-\alpha/2}$ [Ste,
V.1].

2.2. Bessel potentials
The restriction $\alpha<d$ is often inconvenient [Ste V.3]. We can avoid it by
working instead with $(I-\Delta)$, which has positive eigenvalues. For
suitable $f$, the Fourier transforms of $-\Delta f$ and $(I-\Delta)f$ are
$4{\pi}^{2}|x|^{2}\hat{f}(x)$ and $(1+4{\pi}^{2}|x|^{2})\hat{f}(x)$. The
corresponding fractional powers

${\cal{J}}_{\alpha}:=(I-\Delta)^{-\alpha/2}$

are called Bessel potentials (the calligraphic ${\cal{J}}_{\alpha}$ is to
avoid confusion with $J_{\alpha}$, the Bessel function of the first kind). The
kernel of ${\cal{J}}_{\alpha}$ is the function $G_{\alpha}$ with Fourier
transform

$(1+4{\pi}^{2}|x|^{2})^{-\alpha/2}.$

The name Bessel potential derives from the fact that the kernel $G_{\alpha}$
involves a Bessel function of the third kind with imaginary argument, more
briefly a Macdonald function [Wat, 3.7(6)], $K_{\nu}$ say. With all three
kinds of Bessel function, $J_{\nu},I_{\nu},K_{\nu}$, the standard notation for
the order is $\nu$ [Wat].
We now have [Den], [AroS]

$G_{\alpha}(x):=\frac{K_{(d-\alpha)/2}(|x|)\
|x|^{(\alpha-d)/2}}{2^{(d+\alpha-2)/2}{\pi}^{d/2}\ \Gamma(\alpha/2)}.$

The functions in the two last displays are essentially Fourier transforms of
each other; see below.
The dimension $d$ is fixed; the parameter $\alpha$ governs smoothness, as we
shall see. It is convenient, following Simpson. Lindgren and Rue [SimLR1], to
use

$\eta:=(\alpha-d)/2$ $None$

instead as smoothness parameter. We can then write the kernel $G_{\alpha}$ as

$G_{\alpha}(x)=|x|^{\eta}K_{\eta}(|x|)/c_{\eta},\quad
c_{\eta}:=2^{(d+\alpha-2)/2}{\pi}^{d/2}\ \Gamma(\alpha/2).$

Note that $\eta$ may be negative (not a problem, as $K_{-\eta}=K_{\eta}$; see
below).
Macdonald functions.
As Watson remarks, the importance of the Macdonald function is as a Bessel
function with exponential decay at infinity [Wat, 3.7 and 7.23]. The integral
representation

$K_{\nu}(x)=\frac{1}{2}\int_{0}^{\infty}u^{\nu}\ \exp\ \\{x(u+1/u)\\}\ du/u$

is often useful [Erd, II, 7.12(23)], [Jør].
The behaviour of $K_{\nu}$ near the origin is important, as we shall see. We
note here that $K_{0}$ has a logarithmic singularity at the origin [Wat,
(3.71)(14)]: $K_{0}(x)\sim\log(1/x)$ ($x\to 0+$). For $\nu\neq 0$ [Wat,
3.7(2),(6)]:

$K_{\nu}(z):=\frac{1}{2}\pi\
\frac{I_{-{\nu}}(z)-I_{\nu}(z)}{\sin\pi\nu},\qquad
I_{\nu}(z)=\sum_{0}^{\infty}\frac{(\frac{1}{2}z)^{\nu+2m}}{m!\Gamma(\nu+m+1)}.$

(using the obvious ‘L’Hospital convention’ for $\nu=0$), so

$K_{\nu}=K_{-\nu}.$

For small $z$ $K_{\nu}(z)\sim c_{\nu}z^{\nu}$ for constant $c_{\nu}\neq 0$.

2.3. Macdonald-Student pairs
Stein gives the following useful integral formula for $G_{\alpha}$ [Ste,
V.3(26)]:

$G_{\alpha}(x)=\frac{1}{(4\pi)^{\alpha/2}\Gamma(\alpha/2)}\int_{0}^{\infty}e^{-\pi|x|^{2}/u}e^{-u/4\pi}u^{(-d+\alpha)/2}du/u$
$None$

(as he points out, this is just a rephrasing of the definition of the Gamma
function). Thus $G_{\alpha}$ is positive, even, integrable, and decreasing on
$(0,\infty)$. So it is (to within a scale factor) a probability density, let
us call it the Macdonald density. So the Macdonald and Student-$t$ densities
form a Fourier pair (‘self-reciprocal pair’ in Feller’s terminology [Fel,
503]): each is the other’s characteristic function to within scale (as one is
always integrable [Fel, XV.3 Th. 3]). For this Macdonald-Student Fourier pair
(Macdonald-Student pair for brevity), see e.g. Guttorp and Gneiting [GutG].
They give a historical study of of the family below, studied in 1960 by the
Swedish statistician Bertil Matérn (1921-2010) [Mat], and now named after him.

2.4. The Matérn family of correlation functions; Matérn processes
Guttorp and Gneiting [GutG, (1)] give the correlation function in the Matérn
family between values of a spatial random function at locations
$s\in{\mathbb{R}}^{d}$ apart as

${\rho}_{\nu}(s):=\frac{2^{1-\nu}}{\Gamma(\nu)}(\kappa|s|)^{\nu}K_{\nu}(\kappa|s|)\propto\int_{{\mathbb{R}}^{d}}\frac{\exp\\{is^{T}x\\}}{({\kappa}^{2}+|x|^{2})^{(2\nu+d)/2}}\
dx$ $None$

($\kappa>0$ is the scale parameter, $\nu>0$ is their notation for the
smoothness parameter; the variance is scaled so that ${\rho}_{\nu}(0)=1$).
Bochner’s theorem of 1933 shows that characteristic functions of probability
distributions are exactly the continuous positive definite functions taking
the value 1 at the origin. Khintchine’s theorem of 1934 shows that for
stationary processes, the correlation function is the Fourier transform of a
probability distribution, the spectral measure. More is true: by the Cramér
representation of 1942, the process itself is the Fourier transform of a
stochastic process with orthogonal increments (for details and references, see
e.g. Doob [Doo, X, XI], Cramér and Leadbetter [CraL, §7.5]). Thus the Matérn
family of correlations are those of stationary processes, the Matérn
processes, with spectral densities of Student type. Note that here the
correlation function does not depend on the dimension $d$, while the spectral
density does.
Matérn processes have long been widely used, and have become ‘the workhorses
of spatial statistics’. Their theoretical study is thus justified, not only by
its intrinsic interest, but also by the decades of practical experience of
practitioners. This is well summarised by M. L. Stein, whose advice in his
book on spatial data and kriging is ‘use the Matérn model’ [Stei, 14]. He
‘showed that the behaviour of the covariance function near the origin has
fundamental implications on predictive distributions, particularly predictive
uncertainty. The key feature of the Matérn is the inclusion of a smoothness
parameter that directly controls correlation at small distances’, and ‘The
smoothness parameter is aptly named as it implies levels of mean-square
differentiability of the random process, with large $\nu$ yielding very smooth
processes that are many times differentiable, and small $\nu$ yielding rough
processes’. There is also a link with Hausdorff dimension [GenK, §3].

2.5. Bessel potentials and the Matérn family
The kernel $G_{\alpha}$ of ${\cal J}_{\alpha}$ involves a Macdonald function
$K_{\eta}$. Crucial for us is the stochastic partial differential equation
(SPDE)

$({\kappa}^{2}-\Delta)^{\alpha/2}\ X(s)=\sigma\
W(s)\qquad(s\in{\mathbb{R}}^{d}),$ $None$

where $W$ is white noise. We note here that white noise involves generalised
functions (Schwartz distributions) rather than ordinary ones. For background
on white noise analysis, we refer to the standard work by Hida et al.
[HidKPS], and for background on generalized random fields to Gelfand and
Vilenkin [GelV].
Following Simpson, Lindgren and Rue [SimLR1], we write the Matérn covariance
as

$c_{\eta}(s_{i},s_{j})=C_{\eta}(\|s_{i}-s_{j}\|):=\frac{{\sigma}^{2}}{\Gamma(\eta+d/2)(4\pi)^{d/2}{\kappa}^{2\eta}2^{\eta-1}}\
(\kappa\|s_{i}-s_{j}\|)^{\eta}K_{\eta}(\kappa\|s_{i}-s_{j}\|^{2})$

(this differs slightly from the notation of [GutG]: $\kappa$ is a scale
parameter as there; $\sigma$ appears as this is a covariance, not a
correlation; we use $\eta$ for the smoothness parameter here rather than the
traditional $\nu$ to show the link with the dimension $d$ and the Bessel
parameter $\alpha$ in (eta) above).
It is clear from the effect of the Bessel potential being to introduce a
multiplier into the Fourier transform that a convolution is involved, and
clear from the stochastic context that this convolution will involve a
stochastic process. So it is to be expected that the deterministic factor is
the kernel $c_{\eta}$ above and the stochastic factor is white noise.

2.6. Convolution (kernel) representations; Whittle’s theorem
The link between Matérn processes and SPDEs goes back to Whittle [Whi1] [Whi2]
in 1954 and 1963. This was the starting-point of the 2011 study [LinRL], and
its 2012 sequel [SimLR1]. So we shall call the result below Whittle’s theorem;
we follow [SimLR1] in the proof. The result is as in Simpson’s contribution to
the discussion of [LinRL, 65-66].

Theorem W (Whittle’s theorem). The stationary solutions to $(SPDE)$ above are
the convolutions

$X(s)=\int_{{\mathbb{R}}^{d}}c_{\eta}(s,t)\ dW(t)\qquad(s\in{\mathbb{R}}^{d})$

with $c_{\eta}$ the Matérn covariance above and $W$ white noise. That is,

$\int_{{\mathbb{R}}^{d}}ds\ \psi(s)X(s)=\int_{{\mathbb{R}}^{d}}ds\
\psi(s)\int_{{\mathbb{R}}^{d}}c_{\eta}(s,t)\ dW(t)$

for all test functions $\psi$.

Proof. Following Walsh [Wal], we define a solution of $(SPDE)$ to be any
random field $X(s)$ satisfying

$\int X(s)({\kappa}^{2}-\Delta)^{\alpha/2}\phi(s)ds=\sigma\int\phi(s)dW(s)$

for every test function $\phi$. Choose $\phi$ to solve

$({\kappa}^{2}-\Delta)^{\alpha/2}\phi(s)ds=\psi(s):$

then

$\phi(s)=({\kappa}^{2}-\Delta)^{-\alpha/2}\psi(s)$

is smooth [Ste]. Substituting in the above gives

$\int X(s)\psi(s)ds=\sigma\int({\kappa}^{2}-\Delta)^{-\alpha/2}\psi(s)dW(s).$

But as the kernel of the Bessel potential is the Matérn covariance $c_{\eta}$,

$({\kappa}^{2}-\Delta)^{-\alpha/2}\psi(s)=\int c_{\eta}(s,t)\psi(t)dt.$

Substituting and using Fubini’s theorem gives the result. $\square$

As above, Whittle’s result [Whi1], [Whi2] was the starting-point for [LinRL].
The authors proceeded, as in [GutG], to use the above results on Fourier
transforms to identify the finite-dimensional distributions of the two sides
of $(SPDE)$ by showing that they integrate the same way over each finite set
of test functions.
Note. 1. One important advantage of the approach above is that the parameter
$\alpha>0$ is not subject to the restriction to half-integer values in
[LinRL], [SimLR1]. Such situations are typical of this area, where the
relevant analysis and/or special-function theory involves a continuous
parameter, which has geometric significance as a dimension for integer (or
half-integer) values. See e.g. [Bin1,2,3], [BinS1,2].
2\. The representation of the process $X$ in Whittle’s Theorem as a
convolution of the function (kernel) $C_{\eta}$ and the process $W$ is a
prototype of a wider class, the convolution processes (process convolutions,
kernel convolutions). These were introduced by Higdon [Hig1,2] in 1998, in a
climate-modelling setting, and have been widely used since; see e.g. Rodrigues
and Diggle [RodD].
3\. Kernel convolutions have the great advantage that one can generalise by
replacing the Matérn kernel $c_{\eta}$ by any $L_{2}$-kernel $k$:

$X(s)=\int_{{\mathbb{R}}^{d}}k(s,t)\ dW(t)$

still gives a Gaussian random field. This gives one great flexibility when
modelling. One can also replace white noise by any independently scattered
random measure (measures of disjoint sets are independent: [Kal2]), and hence
one can construct non-Gaussian random fields [AbeP].
4\. For more on smoothness of paths (regularity) and SPDEs, see [LanS],
[BinS1], [BroKLO].

2.7. Matern processes on spheres, manifolds and graphs
Much of the theory above, motivated by the needs of spatial statistics, has
involved Fourier analysis in Euclidean space, which is self-dual under the
Fourier transform. While this is the right framework for a locality or region,
it is not on the larger scale (‘maps and globes’). We live on Planet Earth,
which is (approximately) a sphere, which we may scale to be the unit sphere,
$\mathbb{S}:={\mathbb{S}}^{2}$ as in §1.
As on any compact Riemannian manifold $M$, there exists a complete orthonormal
basis $\\{{\phi}_{n}\\}$ of $L_{2}(\mathbb{S})$ of eigenfunctions of
$-\Delta=-{\Delta}_{\mathbb{S}}$, with eigenvalues ${\lambda}_{j}$ with

$0={\lambda}_{0}<{\lambda}_{1}\leq\cdots\leq{\lambda}_{j}\cdots\uparrow+\infty$

One can expand a suitable function on $M$ in an eigenexpansion, the Sturm-
Liouville expansion (see e.g. Chavel [Cha1, VI.1], where this is done for the
heat kernel):

$-\Delta f=\sum_{0}^{\infty}{\lambda}_{n}\langle
f,{\phi}_{n}\rangle{\phi}_{n}.$

For suitable $\Phi$, one can extend this to

$\Phi(-\Delta)f=\sum_{0}^{\infty}\Phi({\lambda}_{n})\langle
f,{\phi}_{n}\rangle{\phi}_{n},$

an example of the functional calculus (or symbolic calculus, of F. Riesz and
N. Dunford); see e.g. [Rud, 290-292].
To exploit the geometry, one expands in a basis of interest; for the sphere
here, one uses the spherical harmonics ([AndAR Ch.9], [SteW, IV.2]) as in
[BinS1 §3].
Using $\cal{F}$ for the Fourier transform, the last display in the proof above
may be summarised as [LinRL, B.3.2]

$\\{{\cal{F}}((\kappa^{2}-\Delta)^{\alpha/2})\phi\\}({\bf
k})=(\kappa^{2}+\|{\bf k}\|^{2})^{\alpha/2}({\cal{F}}\phi)({\bf k})\qquad({\bf
k}\in{\mathbb{R}}^{d}).$

The analogous statement for the sphere is [LinRL, B.3.2]

$\\{{\cal{F}}((\kappa^{2}-\Delta)^{\alpha/2})\phi\\}(k)=(\kappa^{2}+{\lambda}_{k}^{2})^{\alpha/2}({\cal{F}}\phi)(k)\qquad(k=0,1,2,\cdots).$

For such fractional calculus on spheres, see [BinS1, §4 Remark 2] and the
references there (to Askey and Wainger, Bavinck etc.)
For Matérn Gaussian processes on general Riemannian manifolds, see Borovitskiy
et al. [BorTMD]. For graphs, one proceeds similarly using the graph Laplacian
(see e.g. [Chu]). For Matérn processes on graphs, see Borovitskiy et al.
[BorATMDD, §3].

2.8. Extreme values
The area of spatial (or spatio-temporal, or geotemporal) extremes is vast and
topical. For reasons of space, we can only refer here to [AsaDE], [CooNN],
[HaaRo], [ShaW] and the first author’s survey [BinO] with Ostaszewski.

3\. Gaussian Markov random fields; sparse numerical linear algebra

3.1. Precision matrices and sparseness
If the space $\Omega$ is discretised by triangulation, the vertices and edges
of the triangles form a graph, $G=(V,E)$, $V=\\{v_{1},\cdots,v_{n}\\}$ say. We
may thus use the powerful graphical models of statistics; see e.g. Lauritzen
[Lau].
The values $x_{i}$ of our Gaussian process $X$ at the vertices $v_{i}$ are
multivariate Gaussian, with covariance matrix $\Sigma=({\sigma}_{ij})$, say,
and precision (concentration, inverse covariance) matrix
$K:={\Sigma}^{-1}=(k_{ij})$ (‘K for Konzentration’). Two Gaussian variables
$x_{i}$ and $x_{j}$ are independent if and only if their covariance
${\sigma}_{ij}=0$, i.e. they are uncorrelated (for the simple proof, see e.g.
[Kal1, Lemma 13.1]).
That $x_{i}$, $x_{j}$ are conditionally independent given all the others if
and only if $k_{ij}=0$ is much more recent. Now taking $x_{1}$ as the random
$2$-vector and $x_{2}$ as the $(p-2)$-vector we condition on, the Gaussian
Regression Formula (GRF; [Lau, Prop. C5], [BinK, 5.2]) shows us that the
conditional density $f_{1|2}$ of $x_{1}|x_{2}$ has (conditional) covariance
matrix $K_{11}^{-1}$. So one has conditional independence if and only if
$K_{11}^{-1}$ is diagonal, that is (the matrices being $2\times 2$), $K_{11}$
is diagonal, that is, $k_{12}=0$.
This important result derives from Dempster [Dem2, p.159], [Dem1], so as it
needs a name let us call it Dempster’s theorem. It owes its great importance
in modern statistics to the role it plays in graphical models [Lau, Prop.
5.2]. Sparseness properties of concentration matrices are highly revealing
about structure, as in [HasTW], as well as being numerically very convenient
(we note in passing the great importance of numerical linear algebra, here and
elsewhere; see e.g. [GolV]). Conditional independence statements of this sort
are important in Markov properties on undirected graphs (Hammersley-Clifford
theorem; [Lau, 3.2.1]), multivariate normal models [Lau, 5.1.3], covariance
selection models ([Lau, 5.2], following [Dem2]), etc. They come into their own
with Gaussian Markov random fields; see [RueH].
Building on [RueH], [LinRL] proceed by triangulating the space $\Omega$. For
this one can use Delaunay triangulation, giving the dual of the Voronoi
diagram (see e.g. [PrepS, Ch. 5,6]). Here one maximises the minimum angle in a
triangle, avoiding ‘long thin triangles’ which are numerically (and
geometrically) awkward.
Although the flexibility given by the smoothness parameter $\nu$ in the Matérn
model to reflect smoothness of the process being modelled is valuable, ‘$\nu$
usually fixed, since it is poorly identified in typical applications. A more
natural interpretation of the scale parameter $\kappa$ is as a range
parameter’ [LinRL, §2.1]. From the nature of the data, one may be able to
select a range $\rho$ beyond which one can neglect dependence between two
sites (data points $x_{i},x_{j}$ at vertices $v_{i},v_{j}$). An empirical rule
used here is $\rho=\sqrt{8\nu}/\kappa$.
We have stressed the importance of the smoothness parameter $\nu$ (or $\eta$
in 2.2, 2.5), yet as we have just seen it is difficult to estimate from data.
This situation is not surprising, and is reminiscent of that in density
estimation [Sil]. A discrete distribution can be approximated arbitrarily
closely by an arbitrarily smooth one (in any reasonable metric), and vice
versa. So the choice of what smoothness to assume is to be made by the
statistician, with the specifics of the particular situation in mind, and the
resulting flexibility is an important advantage.

3.2. Discretization: From Gaussian process to Gaussian Markov random field
Having chosen the range parameter $\rho$, one can now assume that the data
points distant more than $\rho$ apart are conditionally independent given the
rest, that is, by Dempster’s theorem, that the corresponding entries in the
precision matrix $K$ are zero. As this will be the case for most entries, the
$n\times n$ matrix $K$ is sparse. This in contrast to the covariance matrix
$\Sigma=K^{-1}$, which will be dense (most of its entries non-zero). So matrix
operations on $K$ will have a computational cost of $O(n^{3/2})$, while those
on $\Sigma$ have cost $O(n^{3})$. While the ‘arms race’ between increasing
computing power and increasing size of data sets is ongoing, such a dramatic
difference in computational efficiency is and will remain decisive. This
explains the deliberately dramatic titles of the papers [SimLR1,2] – which in
turn suggested our own.
In the above, only those edges $e_{ij}$ for which the corresponding
$\kappa_{ij}\neq 0$ need to be retained (thus the graph will itself be
‘sparse’, reflecting the inter-point distances). As above, the conditional
independence corresponds to the Markov property. So the continuous Gaussian
process has been approximated by a Gaussian (as all distributions are still
multivariate normal) Markov random field (GMRF); for background on GMRFs, see
Rue and Held [RueH].
Although $\Sigma$ and its inverse $K$ are both symmetric and positive
definite, have a Cholesky decomposition etc., they are diametrically opposite
computationally as only one is sparse. Numerical linear algebra for sparse
matrices is a subject in its own right; see e.g. Davis et al. [DavRS] and its
many references. We are at liberty to permute the order of the vertices
$v_{i}$, and this can lead to worthwhile improvements in computational
efficiency, etc. Statistical learning with sparsity is also a field in its own
right, in which the lasso (least absolute shrinkage and selection operator)
plays a central role. For a monograph treatment, see [HasTW].
Discretization of Riemannian manifolds is widely practised; see e.g. Chavel
[Cha2, §4.4]. There are connections with rough isometries.
There is much more to say here on the numerics side, particularly on the
finite-element method (FEM) and integrated nested Laplace approximation (INLA)
[RueMC], but for reasons of space we must refer to the sources cited above,
particularly the RSS discussion papers [LinRL], [RueMC]. There is also much to
be said about related areas in which sparsity is the key, such as wavelets
[Dau] and compressed sensing [FouR]. We mention some of the many application
areas. The first is digitization and reproduction of images: see e.g. [FouR
Fig. 1] for a photograph of one of the authors’ children, and a near-perfect
reconstruction of it using only the 1% of the wavelet coefficients largest in
magnitude. The second is automatic fingerprint (and more recently, iris)
identification systems; see e.g. [MeyC, Fig. 39-42], [RajWMC, Fig. 18-20].
Related applications include automatic number plate recognition for vehicles,
bank card recognition for cashless transactions, etc.

Dedication and acknowledgement
It is a pleasure to dedicate this paper to Mikhailo Iosifovich Yadrenko on his
90th birthday, in acknowledgement of his long and productive scientific career
and his many contributions to random fields.
We thank the editors for their kind invitation to contribute to this special
volume, and the referees for their helpful comments.

References
[AbeP] S. Aberg and K. Podgórski, A class of non-Gaussian second-order random
fields. Extremes 14 (2011), 187-222.
[AndAR] G. E. Andrews, R. Askey and R. Roy, Special functions. Encycl. Math.
Appl. 71, Cambridge University Press, 1999.
[AroS] N. Aronszajn and K. T. Smith, Theory of Bessel potentials I. Ann. Inst.
Fourier 11 (1961), 385-475.
[AsaDE] P. Asadi, A. C. Davison and S. Engelke, Extremes on river networks.
Ann. Appl. Stat. 9 (2015), 2023-2050.
[AskB] R. Askey and N. H. Bingham, Gaussian processes on compact symmetric
spaces. Z. Wahrscheinlichkeitstheorie verw. Geb. 37 (1976), 127-143.
[BerP] C. Berg and E. Porcu, From Schoenberg coefficients to Schoenberg
functions. Constructive Approximation 45 (2017), 217-241.
[Bin1] N. H. Bingham, Random walk on spheres. Z. Wahrscheinlichkeitstheorie
verw. Geb. 22 (1972), 169-192.
[Bin2] N. H. Bingham, Integral representations for ultraspherical polynomials.
J. London Math. Soc. (2) 6 (1972), 1-11.
[Bin3] N. H. Bingham, Positive definite functions on spheres. Proc. Cambridge
Phil. Soc. 73 (1973), 145-156.
[BinK] N. H. Bingham and W. J. Krzanowski, The penetration of matrix and
linear algebra into multivariate analysis and statistics. Submitted.
[BinMS] N. H. Bingham, A. Mijatović and Tasmin L. Symons, Brownian manifolds,
negative type and geo-temporal covariances. Herbert Heyer Festschrift (ed. D.
Applebaum and H. H. Kuo). Communications in Stochastic Analysis 10 no.4
(2016), 421-432.
[BinO] N. H. Bingham and A. J. Ostaszewski, Extremes and regular variation. R.
A. Doney Festschrift (ed. L. Chaumont and G. Peskir); arXiv:2001.05420.
[BinS1] N. H. Bingham and Tasmin L. Symons, Gaussian random fields on sphere
and sphere cross line. Stochastic Proc. Appl. (Larry Shepp Memorial Issue);
arXiv:1812.02103; https://doi.org/10.1016/j.spa.2019.08.007.
[BinS2] N. H. Bingham and Tasmin L. Symons, Tasmin L. Symons, Integral
representations for ultraspherical polynomials II, arXiv:2101.11809.
[BinS3] N. H. Bingham and Tasmin L. Symons, Tasmin L. Symons, Mathematics of
Planet Earth, I: Geotremporal covariances. arXiv:1706.02972.
[BinS4] N. H. Bingham and Tasmin L. Symons, Tasmin L. Symons, Mathematics of
Planet Earth, II: The Bochner-Godement theorem for symmetric spaces.
arXiv:1707.05204.
[Bog] V. I. Bogachev, Measure theory, Volumes 1,2. Springer, 2007.
[BorATMDD] V. Borovitskiy, I. Azangulov, A. Terenin, P. Mostowsky, M. P.
Diesenroth and N. Durrande, Matérn Gaussian processes on graphs.
arXiv:2010.15538v3.
[BorTMD] V. Borovitskiy, A. Terenin, P. Mostowsky and M. P. Diesenroth, Matérn
Gaussian processes on Riemannian manifolds. arXiv:2006.10160v2.
[Cha1] I. Chavel, Eigenvalues in Riemannian geometry. Academic Press, 1984.
[Cha2] I. Chavel, Riemannian geometry: A modern introduction. Cambridge Tracts
Math. 108, Cambridge University Press, 1993.
[Chu] F. R. K. Chung, Spectral graph theory. CBMS Reg. Conf. Ser. Math. 92,
Amer. Math. Soc., 1997.
[CooNN] D. Cooley, D. Nychka and P. Naveau, Bayesian spatial modelling of
extreme precipitation return levels. J. Amer. Stat. Assoc. 102 (2007),
824-860.
[CraL] H. Cramér and M. R. Leadbetter, Stationary and related stochastic
processes. Sample function properties and their applications. Wiley, 1967.
[Dav] T. A. Davis, S. Rajamanikam and W. M. Sid-Lakhdar, A survey of direct
methods for sparse linear systems. Acta Numerica (2016), 383-566.
[Dem1] A. P. Dempster, Elements of continuous multivariate analysis. Addison-
Wesley, 1969.
[Dem2] A. P. Dempster, Covariance selection. Biometrics 28 (1972), 157-175.
[Den] J. Deny, Review of [AroS]. Mathematical Reviews 26#1485 (MR0143935).
[Dij] G. van Dijk, Introduction to harmonic analysis and generalized Gelfand
pairs. De Gruyter Studies in Math. 36, Walter de Gruyter, 2009.
[Doo] J. L. Doob, Stochastic processes. Wiley, 1953.
[Erd] A. Erdélyi (ed.), Higher transcendental functions, Vol. I-III, McGraw-
Hill, 1953.
[Fel] W. Feller, An introduction to probability theory and its applications,
Volume II, end ed. Wiley, 1971 (1st ed. 1966).
[FouR] S. Foucart and H. Rauhut, A mathematical introduction to compressive
sensing. Springer, 2013.
[Gel] I. M. Gelfand, Spherical functions in symmetric Riemann spaces
(Russian). Doklady Akac. Nauk SSSR 70 (1950), 5-8.
[GelV] I. M. Gelfand and N. Ya. Vilenkin, Generalized functions. Vol. 4.
Applications of harmonic analysis. Academic Press, 1964 (reprinted, AMS
Chelsea, 2016).
[GenK] M. G. Genton and W. Kleiber, Cross-covariance functions for
multivariate geostatistics. Stat. Sci. 30 (2015), 147-163.
[God] R. Godement, Review of [Gel]. Mathematical Reviews 11, 498b (MR0033832).
[GolV] G. H. Golub and C. F. Van Loan, Matrix computations, 3rd ed. The Johns
Hopkins University Press, 1996.
[GueMP] J. C. Guella, V. A. Menegatto and A. P. Peron, An extension of a
theorem of Schoenberg to products of spheres. Banach J. Math. Analysis 10
(2016), 671-685.
[GutG] P. Guttorp and T. Gneiting, Studies in the history of probability and
statistics XLIX: On the Matérn correlation family. Biometrika 93 (2006),
989-995.
[HaaRo] L. de Haan and J. de Ronde, Sea and wind: multivariate extremes at
work. Extremes 1 (1998), 7-45.
[HasTW] T. Hastie, R. Tibshirani and M. Wainwright, Statistical learning with
sparsity: the lasso and its generalizations. Monogr. Stat. Appl. Prob. 143,
Chapman & Hall/CRC 2015.
[Hel1] S. Helgason, Differential geometry and symmetric spaces. Academic
Press, 1962\.
[Hel2] S. Helgason, Differential geometry, Lie groups and symmetric spaces.
Academic Press, 1978.
[Hel3] S. Helgason, Groups and geometric analysis: Integral geometry,
invariant differential operators, and spherical functions. Academic Press,
1978 (2nd ed., Amer. Math. Soc., 2001).
[Hel4] S. Helgason, Geometric analysis on symmetric spaces. Math. Surveys
Monog. 39, Amer. Math. Soc., 1994 (2nd ed. 2008).
[HidKPS] T. Hida, H.-H. Kuo, J. Potthoff and L. Streit, White noise: An
infinite-dimensional calculus. Kluwer, 1993.
[Hig1] D. Higdon, A process-convolution approach to modelling temperature in
the North Atlantic Ocean. Environmental and Ecological Statistics 5 (1998),
173-190.
[Hig2] D. Higdon, Space and space-time modelling using process convolutions.
Quantitative methods for current environmental issues (ed. C. W. Anderson, V.
Barnett, P. C. Chatwin and A. H. El-Shaarwi), 37-56. Springer, 2002.
[HorJ] R. A. Horn and C. R. Johnson, Matrix analysis. Cambridge University
Press, 1985.
[Jør] B. Jørgensen, Statistical properties of the generalized inverse Gaussian
distribution. Lecture Notes in Stat. 9, Springer, 1982.
[Kal1] O. Kallenberg, Foundations of modern probability, 2nd ed. Springer,
2002 (1st ed. 1997).
[Kal2] O. Kallenberg, Random measures, theory and applications. Springer,
2017.
[Lau] S. L. Lauritzen, Graphical models. Oxford University Press, 1996.
[LinRL] F. Lindgren, H. Rue and T. Lindström, An explicit link between
Gaussian fields and Gaussian Markov random fields. The SPDE approach (with
discussion). J. Roy. Stat. Soc. B 73 (2011), 423-498.
[Mal] A. Malyarenko, Invariant random fields on spaces with a group action.
Springer, 2013.
[MarP] D. Marinucci and G. Peccati, Random fields on the sphere:
Representation, limit theorems and cosmological applications. London Math.
Soc. Lecture Notes Series 389, Cambridge University Press, 2011.
[Mat] B. Matérn, Spatial variation, 2nd ed. Springer, 1986 (1st ed. 1960).
[MeyC] F. G. Meyer and R. R. Coifman, Brushlets: A tool for directional image
analysis and image compression. App. Comput. Harm. Anal. 4 (1997), 147-187.
[PrepS] F. P. Preparata and M. I. Shamos, Computational geometry: An
introduction. Springer, 1985.
[RajWMC] N. M. Rajpoot, R. G. Wilson, F. G. Meyer and R. R. Coifman, Adaptive
wavelet packet basis selection for zerotree image coding. IEEE Transations on
Image Processing 12 (2003), 1460-1471.
[Rie] M. Riesz, Intégrales de Riemann-Liouville et potentiels. Acta Szeged
Sect. Math. 9 (1937-38), 1-42; reprinted, Paper 42 p.482-523, in Collected
papers (ed. L. Garding and L. Hörmander). Springer, 1988.
[RodD] A. Rodrigues and P. J. Diggle, A class of convolution-based models for
spatio-temporal processes with non-separable covariance structure. Scand. J.
Stat. 37 (2010), 553-567.
[Rud] W. Rudin, Functional analysis, 2nd ed. McGraw-Hill, 1991 (1st ed. 1973).
[RueH] H. Rue and L. Held, Markov random fields: Theory and applications.
Monographs in Statistics and Applied Probability 104, Chapman & Hall/CRC,
2005.
[RueMC] H. Rue, S. Martino and N. Chopin, Approximate Bayesian inference for
latent Gaussian models by using integrated nested Laplace approximation (with
discussion). J. Roy. Stat. Soc. B 71 (2009), 319-392.
[Sch] I. Schur, Bemerkungen zur Theorie der beschränkten Bilinearformen mit
unendlich vielen Veränderlichen. J. Reine Angew. Math. 140 (1911), 1-28.
[ShaW] P. Sharkey and H. C. Winter, A Bayesian spatial hierarchical model for
extreme precipitation in Great Britain. Environmetrics 30 (2019), e2529.
[Sil] B. W. Silverman, Density estimation for statistics and data analysis.
Chapman & Hall, 1986.
[SimLR1] D. Simpson, F. Lindgren and H. Rue, In order to make spatial
statistics computationally feasible, we need to forget about the covariance
function. Environmetrics 23(1) (Special Issue on Spatio-Temporal Stochastic
Modelling) (2012), 65-74.
[SimLR2] D. Simpson, F. Lindgren and H. Rue, Beyond the valley of the
covariance function. Stat. Sci. 30 (2015), 164-166.
[Ste] E. M. Stein, Singular integrals and differentiability properties of
functions. Princeton University Press, 1970.
[SteW] E. M. Stein and G. Weiss, Introduction to Fourier analysis on Euclidean
spaces. Princeton University Press, 1971.
[Stei] M. L. Stein, Interpolation of spatial data: Some theory for kriging.
Springer, 1999.
[Sze] G. Szegő, Orthogonal polynomials. AMS Colloquium Publications XXIII,
Amer. Math. Soc., 1959.
[Tay] M. E. Taylor, Pseudodifferential operators. Princeton University Press,
1981.
[Wal] J. B. Walsh, An introduction to stochastic partial differential
equations. Ecole d’Eté de Saint-Flour XIV, 265-439; Lecture Notes in Math.
1180, Springer, 1984.
[Wat] G. N. Watson, A treatise on the theory of Bessel functions, 2nd ed.
Cambridge University Press. 1944 (1st ed. 1922).
[Whi1] P. Whittle, On stationary processes in the plane. Biometrika 41 (1954),
434-449.
[Whi2] P. Whittle, Stochastic processes in several dimensions. Bull. Internat.
Stat. Inst. 40 (1963), 974-994.
[Whi3] P. Whittle, Systems in stochastic equilibrium. Wiley, 1986.
[Wil1] D. Williams, Probability with martingales. Cambridge University Press,
1991.
[Wil2] D. Williams, Weighing the odds: A course in probability and statistics.
Cambridge University Press, 2001.
[Wol1] Wolf, J. A.: Spaces of constant curvature. Amer. Math. Soc., 1967 (6th
ed. 2011).
[Wol2] Wolf, J. A.: Harmonic analysis on commutative spaces. Amer. Math. Soc.,
2007.
[Yad] M. I. Yadrenko, Spectral theory of random fields. Optimization Software,
1983.

Mathematics Department, Imperial College, London SW7 2AZ, UK;
<EMAIL_ADDRESS>

Telethon Kids Institute, 15 Hospital Avenue, Perth, WA 6009, Australia;
<EMAIL_ADDRESS>