# Sample Path Regularity of Gaussian Processes
from the Covariance Kernel

Nathaël Da Costa1<EMAIL_ADDRESS>
Marvin Pförtner1<EMAIL_ADDRESS>
Lancelot Da Costa2,3,4<EMAIL_ADDRESS>
Philipp Hennig1<EMAIL_ADDRESS>

1Tübingen AI Center, University of Tübingen, Tübingen, Germany
2VERSES AI Research Lab, Los Angeles, USA
3Department of Mathematics, Imperial College, London, UK
4Wellcome Centre for Human Neuroimaging, UCL, London, UK

###### Abstract

Gaussian processes (GPs) are the most common formalism for defining
probability distributions over spaces of functions. While applications of GPs
are myriad, a comprehensive understanding of GP sample paths, i.e. the
function spaces over which they define a probability measure, is lacking. In
practice, GPs are not constructed through a probability measure, but instead
through a mean function and a covariance kernel. In this paper we provide
necessary and sufficient conditions on the covariance kernel for the sample
paths of the corresponding GP to attain a given regularity. We use the
framework of Hölder regularity as it grants particularly straightforward
conditions, which simplify further in the cases of stationary and isotropic
GPs. We then demonstrate that our results allow for novel and unusually tight
characterisations of the sample path regularities of the GPs commonly used in
machine learning applications, such as the Matérn GPs.

Keywords: Gaussian processes, Gaussian random fields, differentiability,
Hölder regularity, Sobolev regularity, stationary kernels, isotropic kernels,
Matérn kernels

##  1 Introduction

Gaussian processes (GPs) provide a formalism to assign probability
distributions over spaces of functions. That distribution is principally
characterised and controlled by the process’ covariance function, which is a
positive definite kernel. Such kernels are also associated with reproducing
kernel Hilbert spaces (RKHSs). However, it is relatively widely known that the
_sample paths_ of the associated GP are not generally elements of the RKHS,
but form a “larger” space of typically less regular functions (Kanagawa et
al., 2018, Section 4). This sample path space is _much harder to characterise_
than the RKHS.

GP regression is widely applied in statistics and machine learning for
inference from physical observations. In such settings, practitioners mostly
concern themselves with the posterior mean function (which _is_ an element of
the RKHS) and the marginal variance (which is _not_ in the RKHS, but inherits
its regularity) while largely ignoring the regularity of the sample paths and
other properties of the support of the prior probability measure. But recently
relevant use cases for GPs in computational tasks urgently require a more
careful, and ideally tight analysis of the sample path regularity, which we
provide in this work. For instance, Pförtner et al. (2022) recently showed
that a large class of classic solution methods for linear partial differential
equations (PDEs) can be interpreted as GP inference, i.e. as base instances of
probabilistic numerical methods. To infer the solution to a PDE using a
probabilistic numerical method, we want to construct a GP prior for it. It
should be tailored to the problem as tightly as possible, for the following
reasons:

  1. 1.

We need to ensure that the PDE’s differential operator is well-defined on all
samples of the GP. Hence, the sample paths must be regular enough. Here,
regularity typically refers to the existence of a number of strong or weak
partial derivatives, which is encapsulated in the frameworks of Hölder and
Sobolev regularity respectively.

  2. 2.

But we also want the GP posterior to provide a useful uncertainty estimate
over the solution of the PDE. Hence, we do not want to needlessly impose
_additional_ regularity constraints on the sample paths. The sample paths
should be as irregular as possible, to avoid overconfidence. Such cautious
models may also be advisable on numerical grounds, to avoid instabilities such
as Gibbs or Runge phenomena.

The above is just one example of a setting in which one would like to
characterise the regularity of GP sample paths as tightly as possible. Further
examples include generalised coordinates, a method that underlies the analysis
and filtering of stochastic differential equations driven by noise admitting a
high number of derivatives (Friston et al., 2023; Heins et al., 2023). These
noise signals are usually modelled with GPs, and hence characterising GP
sample path regularity is crucial for understanding in which situations these
methods can be applied.

Figure 1: Left: surface plot of a kernel defined as the tensor product of two
Wendland kernels (with $d=1$, $n=0$ on the $x_{1}$-axis and $d=1$, $n=1$ on
the $x_{2}$-axis). Right: surface plot of a sample path from the corresponding
centered GP. The sample path displays different regularity along each axis
(non-differentiable along the $x_{1}$-axis; once differentiable along the
$x_{2}$-axis), reflecting the rotational asymmetry in the regularity of the
kernel.

###  1.1 Summary of Contributions

A consequence of our main result, Theorem 6, may be written as follows:

###### Corollary 1

Let $k\colon\mathbb{R}^{d}\times\mathbb{R}^{d}\to\mathbb{R}$ be a symmetric
positive definite kernel. If either

  * •

all partial derivatives of the form
$\frac{\partial^{2n}k}{\partial\bm{x}^{\bm{\alpha}}\partial\bm{y}^{\bm{\alpha}}}$
for multi-indices $\bm{\alpha}\in\mathbb{N}_{0}^{d}$ with $|\bm{\alpha}|=n$
exist and are (locally) Lipschitz,

  * •

for stationary $k(\bm{x},\bm{y})=k_{\delta}(\bm{x}-\bm{y})$,
$\frac{\partial^{2n}k_{\delta}}{\partial\bm{x}^{\bm{\alpha}}}$ for
$|\bm{\alpha}|=2n$ exist and are (locally) Lipschitz,

  * •

for isotropic $k(\bm{x},\bm{y})=k_{r}(\|\bm{x}-\bm{y}\|)$, $k_{r}^{(2n)}$
exists and is (locally) Lipschitz,

then the sample paths of $f\sim\mathcal{GP}(0,k)$ are $n$ times continuously
differentiable.

Note that, by the mean value theorem, the existence of one additional
continuous derivative of $k$ is sufficient for local Lipschitz continuity of
the lower derivatives.

The actual Theorem 6 is sharper than this. It provides necessary and
sufficient conditions on the regularity of the kernel for the sample paths to
attain a given Hölder regularity. Just as the above corollary, the theorem
encompasses statements for stationary and isotropic GPs.

We apply Theorem 6 in particular to the Matérn GPs to obtain Proposition 9, a
consequence of which is the following:

###### Corollary 2

The sample paths of a centered Matérn GP with smoothness parameter
$\nu\not\in\mathbb{N}_{0}$ are $\lfloor\nu\rfloor$ times continuously
differentiable and no more. In particular, when $\nu=n+\nicefrac{{1}}{{2}}$
for some $n\in\mathbb{N}_{0}$, the sample paths are $n$ times continuously
differentiable and no more.

We then obtain similar results for other GPs, such as the Wendland
(Proposition 11), squared exponential and rational quadratic GPs (Remark 13).

In Section 4.6 we describe how to apply Theorem 6 to covariance kernels
constructed through algebraic transformations of other kernels, including
conic combinations (Proposition 16), products (Proposition 17), tensor
products (Proposition 20) and coordinate transformations (Proposition 22). We
note in Section 4.7 how the results in this work may be generalised to
manifold GPs.

Theorem 24 summarises results about GP sample path Sobolev regularity.
Comparing these with Theorem 6, we argue that one should not expect the
existence of a greater number of weak (Sobolev) derivatives than of strong
derivatives for GP sample paths. This provides further evidence that Theorem 6
contains all that is needed in practice.

We then plot a number of stationary kernel functions as well as sample paths
from the corresponding centered GPs, as in Figure 1.

All regularity properties considered in this work are _local_ , in the sense
that they capture the infinitesimal behaviour of functions around each point
in their domain. Note that our proofs can be extended and our assumptions
strengthened to consider global regularity properties, such as global Hölder
continuity.

###  1.2 Related Work

Some of the earliest work on sample path regularity can be attributed to
Kolmogorov and Chentsov (1956). Later, Fernique (1975) proved the first
necessary and sufficient conditions for sample path continuity of stationary
GPs. Since then, many authors have extended these continuity results. Potthoff
(2009) provides an approachable overview of how such results can be obtained.
It also investigates conditions for uniform continuity, which is a _global_
regularity property. Azmoodeh et al. (2014), which we rely on in this work,
proves a converse to the Kolmogorov continuity theorem for GPs, providing
necessary and sufficient conditions for Hölder continuity.

Differentiability has also been investigated by many authors, including
Scheuerer (2010b) for general random fields, Adler and Taylor (2007), Potthoff
(2010) and Henderson (2024) for Sobolev regularity. The sample path regularity
of Matérn GPs was studied in Scheuerer (2010a) (see specifically Scheuerer
(2010a, Examples 5.3.19, 5.5.6 & 5.5.12)) and of tensor product of one
dimensional Matérn GPs in Wang et al. (2021). Pierce (1970, Proposition I.3)
combines differentiability results with the Kolmogorov continuity theorem, and
can be seen as the global version of the forward implication of our Theorem 6
(1).

Compared to the existing literature, the present work provides _general_ ,
_versatile_ , _easily applicable_ and _unusually tight_ results for the study
of GP sample path continuity and differentiability.

##  2 Preliminaries

###### Definition 3

A Gaussian process on a set $O$ is a map $f\colon O\times\Omega\to\mathbb{R}$,
where $\Omega$ is a probability space, such that for all $x\in O$,
$f(x,\cdot)\colon O\to\mathbb{R}$ is measurable, and such that for each
$X:=(x_{1},\dots,x_{N})\in O^{N}$, the map
$f(X,\cdot)\colon\Omega\to\mathbb{R}^{N}$ given by
$f(X,\omega)=(f(x_{1},\omega),\dots,f(x_{N},\omega))$ is a multivariate
Gaussian random variable.

The maps $f(\cdot,\omega)\colon O\to\mathbb{R}$ for $\omega\in\Omega$ are the
sample paths of the GP.

We would like to study the regularity of the sample paths of a GP. However,
Definition 3 is not convenient to work with in practice: one rarely constructs
the probability space $\Omega$ of a GP. Instead, one characterises a GP by its
_mean_ and its _covariance kernel_.

###### Definition 4

The mean $\mu$ of a GP $f\colon O\times\Omega\to\mathbb{R}$ is the map

$\mu\colon O\to\mathbb{R},\;x\mapsto\mathbb{E}(f(x,\cdot)).$

The covariance kernel $k$ is the map

$k\colon O\times
O\to\mathbb{R},\;(x,y)\mapsto\operatorname{cov}(f(x,\cdot),f(y,\cdot)).$

The GP $f$ is said to be centered if $\mu=0$.

Conversely, given maps $\mu\colon O\to\mathbb{R}$ and $k\colon O\times
O\to\mathbb{R}$ with $k$ symmetric positive definite – meaning that for any
finite set of points $\\{x_{1},\dots,x_{N}\\}\subset\Omega$, the matrix
$(k(x_{i},x_{j}))_{i,j}$ is symmetric positive semi-definite – one can
construct a probability space $\Omega$, and a GP $f\colon
O\times\Omega\to\mathbb{R}$ with mean $\mu$ and covariance kernel $k$ (Klenke,
2014, Theorem 14.36). We then write $f\sim\mathcal{GP}(\mu,k)$.

At this point it is important to note that $f\sim\mathcal{GP}(\mu,k)$ does not
uniquely specify $f$ or the probability space $\Omega$. However it uniquely
specifies the finite dimensional distributions of $f$ – since multivariate
Gaussians are entirely characterised by their first two moments. In machine
learning applications all that will ever be observed are evaluations of $f$ at
finitely many points. Therefore we would like to study sample path regularity
of $f$ _up to modification_. By this we mean that $f\sim\mathcal{GP}(\mu,k)$
will be said to have samples in $\mathcal{F}(O)$, where $\mathcal{F}(O)$ is
some space of functions $O\to\mathbb{R}$, if there exists a construction of
the GP $f$, say $f\colon O\times\Omega\to\mathbb{R}$, such that
$f(\cdot,\omega)\in\mathcal{F}(O)$ for all $\omega\in\Omega$.

Finally, note that $f\sim\mathcal{GP}(\mu,k)$ having samples in
$\mathcal{F}(O)$ is equivalent to $\tilde{f}+\mu$ having samples in
$\mathcal{F}(O)$, where $\tilde{f}\sim\mathcal{GP}(0,k)$. For this it suffices
that $\mu\in\mathcal{F}(O)$ and $\tilde{f}$ has samples in $\mathcal{F}(O)$.
Therefore, in what follows, we will consider only centered GPs.

Thus we are studying sample path regularity of GPs _from the covariance
kernel_ , and our goal is to link the regularity of the kernel to the
regularity of the GP sample paths.

###  2.1 Setup

We study GPs defined on open subsets of Euclidean spaces:
$O\subset\mathbb{R}^{d}$ is open, for some $d\in\mathbb{N}$, and
$f\sim\mathcal{GP}(0,k)$ is a centered GP on $O$. An important special case is
when $k$ is stationary, i.e. $O=\mathbb{R}^{d}$ and there is an even function
$k_{\delta}\colon\mathbb{R}^{d}\to\mathbb{R}$ such that

$k(\bm{x},\bm{y})=k_{\delta}(\bm{x}-\bm{y})$ (1)

for all $\bm{x},\bm{y}\in\mathbb{R}^{d}$. A further special case of $k$
stationary is when $k$ is isotropic, i.e. $O=\mathbb{R}^{d}$ and there is an
even function $k_{r}\colon\mathbb{R}\to\mathbb{R}$ such that

$k(\bm{x},\bm{y})=k_{r}(\|\bm{x}-\bm{y}\|)$ (2)

for all $\bm{x},\bm{y}\in\mathbb{R}^{d}$, where $\|\cdot\|$ denotes the
standard Euclidean norm on $\mathbb{R}^{d}$.

##  3 Hölder Regularity

The main type of sample path regularity we would like to investigate is
continuity and the order of continuous differentiability. These are
encapsulated in the framework of _Hölder regularity_.

We choose $O$ to be an open set in $\mathbb{R}^{d}$ in order to capture
differentiability of the sample paths. However this implies that Hölder spaces
on $O$ are inconvenient for our purposes, as they do not solely characterise
local regularity properties of functions on $O$. Indeed, the Hölder spaces not
only constrain the local behaviour of functions, but also their behaviour at
infinity, or near the boundary of $O$. To retain solely local constraints we
employ local Hölder spaces.

######  Definition 5 (Local Hölder and almost-Hölder spaces)

Let $n\in\mathbb{N}_{0}$ and $\gamma\in[0,1]$.

  1. (1)

The local Hölder space $C^{n,\gamma}_{loc}(O)$ is the space of functions $f$
on $O$ for which $\partial^{\bm{\alpha}}f$ exists for all multi-indices
$\bm{\alpha}=(\alpha_{1},\dots,\alpha_{d})\in\mathbb{N}_{0}^{d}$ with
$|\bm{\alpha}|:=\alpha_{1}+\dots+\alpha_{d}\leq n$, and such that the highest
order partial derivatives satisfy a Hölder condition of the form: for all
compact subsets $K\subset O$ there is a constant $C_{K}>0$ such that

$|\partial^{\bm{\alpha}}f(\bm{x})-\partial^{\bm{\alpha}}f(\bm{y})|\leq
C_{K}\|\bm{x}-\bm{y}\|^{\gamma}$

for all $\bm{x},\bm{y}\in K$ and $|\bm{\alpha}|=n$.

  2. (2)

The local almost-Hölder space $C^{(n+\gamma)^{-}}_{loc}(O)$ is defined as
$\bigcap_{n^{\prime}+\gamma^{\prime}<n+\gamma}C^{n^{\prime},\gamma^{\prime}}_{loc}(O)$.

###### Remark 6

For $n,n^{\prime}\in\mathbb{N}_{0}$ and $\gamma,\gamma^{\prime}\in[0,1]$,
$n^{\prime}+\gamma^{\prime}<n+\gamma$ implies
$C^{n^{\prime},\gamma^{\prime}}_{loc}(O)\supsetneq
C^{(n+\gamma)^{-}}_{loc}(O)\supsetneq C^{n,\gamma}_{loc}(O)$. Moreover
$C^{n}(O)=C^{n,0}_{loc}(O)$.

The key point in the definition of $C^{(n+\gamma)^{-}}_{loc}(O)$ is the strict
inequality $n^{\prime}+\gamma^{\prime}<n+\gamma$, since allowing for equality
we have $\bigcap_{n^{\prime}+\gamma^{\prime}\leq
n+\gamma}C^{n^{\prime},\gamma^{\prime}}_{loc}(O)=C^{n,\gamma}_{loc}(O)$.

As will become clear, the local almost-Hölder spaces are natural for
characterising GP sample path regularity.

To describe the regularity of the covariance kernel $k$, we will require
derivatives of the form $\partial^{\bm{\alpha},\bm{\beta}}k$. Here
$\partial^{\bm{\alpha},\bm{\beta}}$ stands for the application of
$\partial^{\bm{\alpha}}$ with respect to the first variable, followed by the
application of $\partial^{\bm{\beta}}$ with respect to the second variable. In
fact the order in which they are applied will not matter, by continuity of
these partial derivatives. We moreover write $C^{n\otimes n}(O\times O)$ for
the space of functions $k$ on $O\times O$ for which
$\partial^{\bm{\alpha},\bm{\beta}}k$ exists and is continuous for all
$|\bm{\alpha}|,|\bm{\beta}|\leq n$.

Also recall the “big O” notation: the functions
$f,g\colon\mathbb{R}^{d}\setminus\\{\bm{0}\\}\to\mathbb{R}$ satisfy
$f(\bm{h})=\mathcal{O}(g(\bm{h}))$ as $\bm{h}\to\bm{0}$ if and only if there
is $C>0$ such that
$\limsup_{\bm{h}\to\bm{0}}\left|f(\bm{h})/g(\bm{h})\right|\leq C$. Moreover,
for a family of functions $(f_{\bm{x}})_{\bm{x}\in O}$, we say
$f_{\bm{x}}(\bm{h})=\mathcal{O}(g(\bm{h}))$ as $\bm{h}\to\bm{0}$ locally
uniformly in $\bm{x}\in O$ if, for every compact subset $K\subset O$, $C>0$
may be chosen independently of $\bm{x}\in K$.

We are now in a position to state the main result of this paper. [Sample path
Hölder regularity]theoremthmholder Let $n\in\mathbb{N}_{0}$ and
$\gamma\in(0,1]$. The process $f\sim\mathcal{GP}(0,k)$ has samples in
$C^{(n+\gamma)^{-}}_{loc}(O)$ if and only if,

  1. (1)

for general $k$,

     * •

$k\in C^{n\otimes n}(O\times O)$,

     * •

$|\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x}+\bm{h},\bm{x}+\bm{h})-\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x}+\bm{h},\bm{x})-\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x},\bm{x}+\bm{h})+\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x},\bm{x})|=\mathcal{O}(\|\bm{h}\|^{2\epsilon})$
as $\bm{h}\to\bm{0}$, locally uniformly in $\bm{x}\in O$, for all
$\epsilon\in(0,\gamma)$ and $|\bm{\alpha}|=|\bm{\beta}|=n$.

  2. (2)

for stationary $k(\bm{x},\bm{y})=k_{\delta}(\bm{x}-\bm{y})$,

     * •

$k_{\delta}\in C^{2n}(\mathbb{R}^{d})$,

     * •

$|\partial^{\bm{\alpha}}k_{\delta}(\bm{h})-\partial^{\bm{\alpha}}k_{\delta}(\bm{0})|=\mathcal{O}(\|\bm{h}\|^{2\epsilon})$
as $\bm{h}\to\bm{0}$ for all $\epsilon\in(0,\gamma)$ and $|\bm{\alpha}|=2n$.

  3. (3)

for isotropic $k(\bm{x},\bm{y})=k_{r}(\|\bm{x}-\bm{y}\|)$,

     * •

$k_{r}\in C^{2n}(\mathbb{R})$,

     * •

$|k_{r}^{(2n)}(h)-k_{r}^{(2n)}(0)|=\mathcal{O}(|h|^{2\epsilon})$ as $h\to 0$
for all $\epsilon\in(0,\gamma)$.

In each case, differentiating sample-wise we have
$\partial^{\bm{\alpha}}f\sim\mathcal{GP}(0,\partial^{\bm{\alpha},\bm{\alpha}}k)$
for all $|\bm{\alpha}|\leq n$. The proof can be found in Appendix A.

###### Remark 7

The positive definiteness of $k$ allows us to extrapolate regularity: $k$,
$k_{\delta}$ or $k_{r}$ have no more derivatives around the diagonal
$\\{(\bm{x},\bm{x})\in O\times O:\bm{x}\in O\\}$, $\bm{0}$ or $0$ respectively
as they do on the rest of the domain (see Lemma A.3). Therefore, when applying
Theorem 6, it suffices to check the existence of the respective derivatives
around the diagonal, $\bm{0}$ or $0$ respectively.

Another way to see that information about $k$ in a neighbourhood of the
diagonal must be sufficient to characterise any _local_ sample path regularity
property, as it must suffice to consider the covariances at arbitrarily close
by points.

Theorem 6 combines results about sample Hölder continuity of GPs (Azmoodeh et
al., 2014) to results about sample differentiability of GPs (Potthoff, 2010),
applied inductively on the partial derivatives. It provides necessary and
sufficient conditions for $f$ to have almost-Hölder continuous samples of a
certain degree. Let us note however that it does not quite give us necessary
conditions for $f$ to have samples in $C^{n}(O)$, as the latter is not an
almost-Hölder space. To achieve this, we can adapt the proof of Theorem 6
combining necessary and sufficient conditions for sample continuity of GPs
with the results about sample differentiability from Potthoff (2010) (see
Scheuerer (2010a, Theorem 5.3.16) for a result of this flavour).
Characterising necessary and sufficient conditions on the kernel for the
sample continuity of the GP is much more involved than almost-Hölder
continuity, and therefore we do not expand our results to this setting (see
Fernique (1975) for the stationary case, Talagrand (1987) for the general
case). Moreover we demonstrate in the examples below that Theorem 6 is all we
need.

##  4 Examples

In this section we demonstrate how Theorem 6 can be applied in practice. We
investigate the sample path regularity of the most widely used GP families,
recovering known results as well as proving novel ones.

The reverse implications in Theorem 6 allow us to obtain sharp sample path
regularity characterisations in the following sense: $f$ is said to have
samples in $C^{(n+\gamma)^{-}}_{loc}(O)$ _and no more_ , if $f$ is not sample
$C^{(n^{\prime}+\gamma^{\prime})^{-}}_{loc}(O)$ for any
$n^{\prime}+\gamma^{\prime}>n+\gamma$.

###  4.1 Wiener Kernel

The Wiener process is a centered GP on $O=\mathbb{R}_{>0}$ with covariance
kernel

$k(x,y)=\min(x,y)$

for $x,y\in\mathbb{R}_{>0}$. $\min(\cdot,\cdot)$ is Lipschitz but non-
differentiable in both its arguments. So by Theorem 6 (1) the Wiener process
has samples in $C^{\nicefrac{{1}}{{2}}^{-}}_{loc}(\mathbb{R}_{>0})$ and no
more.

###### Remark 8

By the fundamental theorem of calculus, the $n$-times integrated Wiener
process (Schober et al., 2014, Section B) has samples in
$C^{(n+\nicefrac{{1}}{{2}})^{-}}_{loc}(\mathbb{R}_{>0})$ and no more.

###  4.2 Matérn Kernels

The Matérn kernels are isotropic kernels on $\mathbb{R}^{d}$ given by

$k(\bm{x},\bm{y})=k_{r}(\|\bm{x}-\bm{y}\|)=\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\sqrt{2\nu}\|\bm{x}-\bm{y}\|\right)^{\nu}K_{\nu}\left(\sqrt{2\nu}\|\bm{x}-\bm{y}\|\right)$
(3)

for $\bm{x},\bm{y}\in\mathbb{R}^{d}$, where $K_{\nu}$ is the modified Bessel
function of the second kind and $\nu>0$ is the smoothness parameter. The cases
where $\nu=n+\nicefrac{{1}}{{2}}$ for some $n\in\mathbb{N}_{0}$ are
particularly interesting because the expression in Equation (3) simplifies to
a product of a polynomial and an exponential of the radial distance (Rasmussen
and Williams, 2005, Equation 4.16).

###### Proposition 9

A centered Matérn GP with smoothness parameter $\nu>0$ has samples in
$C^{\nu^{-}}_{loc}(O)$ and no more. In particular, if
$\nu=n+\nicefrac{{1}}{{2}}$ for some $n\in\mathbb{N}_{0}$, the GP has samples
in $C^{(n+\nicefrac{{1}}{{2}})^{-}}_{loc}(\mathbb{R}^{d})$ and no more.

This follows from Theorem 6 (3) and the proof can be found in Appendix B.1.

###### Remark 10

Proposition 9 also covers the (multivariate) Ornstein-Uhlenbeck process, the
centered GP with covariance kernel given by
$k(\bm{x},\bm{y})=\exp(-\|\bm{x}-\bm{y}\|)$ which corresponds to the Matérn GP
with $\nu=\nicefrac{{1}}{{2}}$. So the (multivariate) Ornstein-Uhlenbeck
process has samples in $C^{\nicefrac{{1}}{{2}}^{-}}_{loc}(\mathbb{R}^{d})$ and
no more.

###  4.3 Wendland Kernels

The Wendland kernels are isotropic kernels which are compactly supported and
piecewise polynomial in the radial distance. On $\mathbb{R}^{d}$, they are
defined as follows (Wendland, 2004, Definition 9.11):

$k(\bm{x},\bm{y})=k_{r}(\|\bm{x}-\bm{y}\|)=\mathcal{I}^{n}\phi_{\lfloor\nicefrac{{d}}{{2}}\rfloor+n+1}(\|\bm{x}-\bm{y}\|)$

for $\bm{x},\bm{y}\in\mathbb{R}^{d}$, where
$\phi_{j}(\rho):=\max((1-\rho)^{j},0)$ and
$\mathcal{I}\phi(\rho):=\int_{\rho}^{\infty}t\phi(t)\,dt/\int_{0}^{\infty}t\phi(t)\,dt$
for $\rho\geq 0$. $n\in\mathbb{N}_{0}$ is a parameter controlling the degree
of the polynomial (precisely, $k_{r}(x)$ has degree
$\lfloor\nicefrac{{d}}{{2}}\rfloor+3n+1$ in $\rho$ around 0).

###### Proposition 11

A centered Wendland GP with degree parameter $n$ has samples in
$C^{(n+\nicefrac{{1}}{{2}})^{-}}_{loc}(O)$ and no more.

This follows from Theorem 6 (3) and the proof can be found in Appendix B.2.

###  4.4 Smooth Kernels

From Theorem 6 immediately follows the following important corollary:

###### Corollary 12

The process $f\sim\mathcal{GP}(0,k)$ has samples in $C^{\infty}(O)$ if and
only if,

  1. (1)

for general $k$, $k\in C^{\infty}(O\times O)$.

  2. (2)

for stationary $k(\bm{x},\bm{y})=k_{\delta}(\bm{x}-\bm{y})$, $k_{\delta}\in
C^{\infty}(\mathbb{R}^{d})$.

  3. (3)

for isotropic $k(\bm{x},\bm{y})=k_{r}(\|\bm{x}-\bm{y}\|)$, $k_{r}\in
C^{\infty}(\mathbb{R})$.

###### Remark 13

The centered GPs with squared exponential
($k(\bm{x},\bm{y})=\exp(-\|\bm{x}-\bm{y}\|^{2})$), rational quadratic
($k(\bm{x},\bm{y})=(1+\|\bm{x}-\bm{y}\|^{2})^{-a}$ for some $a>0$) and
periodic ($k(x,y)=\exp(-\sin(\pi(x-y))^{2})$ for $x,y\in\mathbb{R}$)
covariance kernels are therefore sample $C^{\infty}(\mathbb{R}^{d})$.

###  4.5 Feature Kernels

For any map $\bm{\phi}\colon O\to\mathbb{R}^{m}$ is associated a feature
kernel

$k(\bm{x},\bm{y})=\bm{\phi}(\bm{x})^{\top}\bm{\phi}(\bm{y})$

for $\bm{x},\bm{y}\in O$. From Theorem 6 (1) we obtain the following result:

###### Proposition 14

For $n\in\mathbb{N}_{0}$ and $\gamma\in(0,1]$, if $\phi_{i}\in
C^{(n+\gamma)^{-}}_{loc}(O)$ for all $1\leq i\leq m$ then
$f\sim\mathcal{GP}(0,k)$ has samples in $C^{(n+\gamma)^{-}}_{loc}(O)$.

The proof can be found in Appendix B.3.

###### Remark 15

As a special case, linear kernels
($k(\bm{x},\bm{y})=\langle\bm{x},\bm{y}\rangle$) and polynomial kernels
($k(\bm{x},\bm{y})=(1+\langle\bm{x},\bm{y}\rangle)^{m}$ for some
$m\in\mathbb{N}$) are feature kernels, with the feature maps being the
coordinate maps or polynomials of the coordinate maps respectively. Therefore
the associated centered GPs have samples in $C^{\infty}(\mathbb{R}^{d})$.

###  4.6 Kernel Algebra

The versatility of the GP framework in applications is largely due to the fact
that covariance kernels can be recursively combined using algebraic operations
like pointwise sums or products which allows for constructing highly flexible
and complex prior models in a systematic and interpretable fashion. It is
therefore important to understand how the sample path regularity of a GP
depends on the recursive structure of the covariance kernel.

In the following the $k_{i}$ are always assumed to be symmetric positive
definite kernels. Note that in this section we work solely with general
kernels, leveraging Theorem 6 (1). It would then be straightforward to
similarly leverage Theorem 6 (2) and Theorem 6 (3) for stationary and
isotropic kernels respectively.

####  4.6.1 Conic Combinations

Conic combinations of kernels are pointwise linear combinations of the kernel
functions with positive weights, i.e.

$k(\bm{x},\bm{y})=\sum_{i=1}^{m}a_{i}k_{i}(\bm{x},\bm{y})$

for $\bm{x},\bm{y}\in O$, where $a_{i}>0$ for all $1\leq i\leq m$. Sums of
kernels and positive scaling of kernels are special cases of conic
combinations with $a_{i}=1$ and $m=1$, respectively.

###### Proposition 16

For $n\in\mathbb{N}_{0}$ and $\gamma\in(0,1]$, if the $k_{i}$ satisfy the
condition in Theorem 6 (1) for all $1\leq i\leq m$ then
$f\sim\mathcal{GP}(0,k)$ has samples in $C^{(n+\gamma)^{-}}_{loc}(O)$.

This follows from the fact that for $f\sim\mathcal{GP}(0,k)$, we can express
the sample paths as
$f(\cdot,\omega)=\sum_{i=1}^{m}\sqrt{a_{i}}f_{i}(\cdot,\omega)$ where
$f_{i}\sim\mathcal{GP}(0,k_{i})$, and that $C^{(n+\gamma)^{-}}_{loc}(O)$ is a
vector space.

####  4.6.2 Products

Define the product kernel

$k(\bm{x},\bm{y})=\prod_{i=1}^{m}k_{i}(\bm{x},\bm{y})$

for $\bm{x},\bm{y}\in O$. $k$ is also positive definite (Berg et al., 1984,
Chapter 3 Theorem 1.12).

###### Proposition 17

For $n\in\mathbb{N}_{0}$ and $\gamma\in(0,1]$, if the $k_{i}$ satisfy the
condition in Theorem 6 (1) for all $1\leq i\leq m$ then
$f\sim\mathcal{GP}(0,k)$ has samples in $C^{(n+\gamma)^{-}}_{loc}(O)$.

The proof can be found in Appendix B.4.

###### Remark 18

Combining Proposition 17 with Proposition 14 we get that if $\tilde{k}$
satisfies the condition in Theorem 6 (1) for some $n\in\mathbb{N}_{0}$ and
$\gamma\in(0,1]$, and $\phi\in C^{(n+\gamma)^{-}}_{loc}(O)$, then defining
$k(\bm{x},\bm{y})=\phi(\bm{x})\tilde{k}(\bm{x},\bm{y})\phi(\bm{y})$ for all
$\bm{x},\bm{y}\in O$, $f\sim\mathcal{GP}(0,k)$ has samples in
$C^{(n+\gamma)^{-}}_{loc}(O)$. We could equivalently prove this by noting that
we can express the sample paths as
$f(\cdot,\omega)=\phi(\cdot)\tilde{f}(\cdot,\omega)$ where
$\tilde{f}\sim\mathcal{GP}(0,\tilde{k})$.

####  4.6.3 Tensor Products

Let $O_{i}\subset\mathbb{R}^{d_{i}}$ be open sets for $1\leq i\leq m$ and
$O_{1}\times\dots\times
O_{m}=:O\subset\mathbb{R}^{d}=\mathbb{R}^{d_{1}}\times\dots\times\mathbb{R}^{d_{m}}$.
For $k_{i}\colon O_{i}\times O_{i}\to\mathbb{R}$ for $1\leq i\leq m$, the
tensor product kernel $k\colon O\times O\to\mathbb{R}$, is given by

$k(\bm{x},\bm{y})=(k_{1}\otimes\dotsb\otimes
k_{m})(\bm{x},\bm{y}):=\prod_{i=1}^{m}k_{i}(\bm{x}_{i},\bm{y}_{i})$

for $\bm{x},\bm{y}\in\mathbb{R}^{d}$, where
$\bm{x}_{i},\bm{y}_{i}\in\mathbb{R}^{d_{i}}$ are the orthogonal projections of
$\bm{x}$, $\bm{y}$ respectively onto $\mathbb{R}^{d_{i}}$.

GPs with tensor product covariance kernels are interesting as they allow
_spherically asymmetric_ sample path regularity. We need the following
definition to formalise this.

###### Definition 19

Let $n_{1},\dots,n_{m}\in\mathbb{N}_{0}$ and
$\gamma_{1},\dots,\gamma_{m}\in[0,1]$.

  1. (1)

$C^{(n_{1},\gamma_{1})\otimes\dotsb\otimes(n_{m},\gamma_{m})}_{loc}(O_{1}\times\dotsb\times
O_{m})$ is the space of functions $f$ on $O$ for which
$\partial^{\bm{\alpha}}f$ exists for all multi-indices
$\bm{\alpha}\in\mathbb{N}_{0}^{d}$ with $|\bm{\alpha}_{i}|\leq n_{i}$, where
$\bm{\alpha}_{i}\in\mathbb{N}_{0}^{d_{i}}$ is the projection of $\bm{\alpha}$
onto $\mathbb{N}_{0}^{d_{i}}$, and such that the highest order partial
derivatives satisfy a Hölder condition of the form: for all compact subsets
$K\subset O$ there is a constant $C_{K}>0$ such that

$|\partial^{\bm{\alpha}}f(\bm{x})-\partial^{\bm{\alpha}}f(\bm{y})|\leq
C_{K}\|\bm{x}_{i}-\bm{y}_{i}\|^{\gamma_{i}}$

for all $1\leq i\leq d$, $\bm{x},\bm{y}\in K$ with $\bm{x}_{j}=\bm{y}_{j}$ for
$j\neq i$, and $|\bm{\alpha}_{i}|=n_{i}$.

  2. (2)

$C^{(n_{1}+\gamma_{1})^{-}\otimes\dotsb\otimes(n_{m}+\gamma_{m})^{-}}_{loc}(O_{1}\times\dotsb\times
O_{m})$
$:=\bigcap_{n_{i}^{\prime}+\gamma_{i}^{\prime}<n_{i}+\gamma_{i}}C^{(n_{1}^{\prime},\gamma_{1}^{\prime})\otimes\dotsb\otimes(n_{m}^{\prime},\gamma_{m}^{\prime})}_{loc}(O_{1}\times\dotsb\times
O_{m})$.

The following proposition generalises Theorem 6.

###### Proposition 20

For $n_{1},\dots,n_{m}\in\mathbb{N}_{0}$ and
$\gamma_{1},\dots,\gamma_{m}\in(0,1]$, $f\sim\mathcal{GP}(0,k)$ has samples in
$C^{(n_{1}+\gamma_{1})^{-}\otimes\dots\otimes(n_{m}+\gamma_{m})^{-}}_{loc}(O_{1}\times\dotsb\times
O_{m})$ if and only if the $k_{i}$ satisfy the condition in Theorem 6 (1) with
$n=n_{i}$ and $\gamma=\gamma_{i}$ for all $1\leq i\leq m$.

The proof of Proposition 20 follows the same steps as the one of Theorem 6 in
Appendix A.

###### Remark 21

We could generalise Theorem 6 to allow for even more flexibility in the
asymmetrical derivatives in Theorem 6. This can be done by defining Hölder
spaces extending the spaces $C^{A}(O)$, where $A\subset\mathbb{N}_{0}^{d}$ is
a downward closed set of multi-indices which specifies which partial
derivatives exist (Pförtner et al., 2022, Definition B.9).

####  4.6.4 Coordinate Transformations

For a map $\bm{\varphi}\colon O\to\mathbb{R}^{m}$ and a symmetric positive
definite kernel $\tilde{k}$ defined on
$\operatorname{Im}(\bm{\varphi})\times\operatorname{Im}(\bm{\varphi})$, we
have a kernel

$k(\bm{x},\bm{y})=\tilde{k}(\bm{\varphi}(\bm{x}),\bm{\varphi}(\bm{y})).$

###### Proposition 22

For $n\in\mathbb{N}_{0}$ and $\gamma,\delta\in(0,1]$, if $\tilde{k}$ satisfies
the condition in Theorem 6 (1) and $\varphi_{i}\in
C^{(n+\delta)^{-}}_{loc}(O)$ for all $1\leq i\leq m$ then
$f\sim\mathcal{GP}(0,k)$ has samples in $C^{(n+\gamma\delta)^{-}}_{loc}(O)$.

This follows from the fact that for $\tilde{f}\sim GP(0,\tilde{k})$ we have
$\tilde{f}\circ\bm{\varphi}\sim\mathcal{GP}(0,k)$, and that the composition of
a $C^{n,\zeta}_{loc}(O)$ with a $C^{n,\epsilon}_{loc}(O)$ function is
$C^{n,\epsilon\zeta}_{loc}(O)$.

###  4.7 Manifold Kernels

We can generalise the results in the present work to manifold domains; instead
of considering GPs on open subsets $O\subset\mathbb{R}^{d}$, we consider GPs
on manifolds $M$. There exists two ways of constructing positive definite
kernels, and hence GPs, on manifolds $M$.

_Extrinsic kernels_ are defined by viewing $M\subset\mathbb{R}^{d}$, taking a
positive definite kernel $k$ on $\mathbb{R}^{d}$, and restricting it to $M$.
This restriction is analogous to a coordinate transformation (Section 4.6.4),
since $k|_{M\times M}(\bm{x},\bm{y})=k(\bm{\iota}(\bm{x}),\bm{\iota}(\bm{y}))$
for $\bm{x},\bm{y}\in M$, where $\bm{\iota}:M\to\mathbb{R}^{d}$ is the
inclusion map. Therefore, if the manifold is smooth, $k|_{M\times M}$ has the
same sample regularity as $k$.

_Intrinsic kernels_ are constructed directly on the manifold. Examples of
these include the intrinsic Matérn and heat kernels (Borovitskiy et al.,
2020), and the hyperbolic secant kernel (Da Costa et al., 2023). In this case,
since the sample path regularity results in this work are all local, we may
study sample path regularity in each coordinate patch $O\subset M$, treating
them as open subsets of $\mathbb{R}^{d}$.

##  5 Sobolev Regularity

In this section we investigate how the regularity of the kernel affects the
weak differentiability of the sample paths. Specifically, we consider the
important $L^{2}$-Sobolev regularity. The more general $L^{p}$-Sobolev
regularity is studied in Henderson (2024). To characterise only the local
regularity of the sample paths, as for Hölder spaces, we define the local pre-
Sobolev spaces.

######  Definition 23 (Local pre-Sobolev spaces)

Let $n\in\mathbb{N}_{0}$. The local pre-Sobolev space
$\mathcal{H}^{n}_{loc}(O)$ is the space of functions111A (local) Sobolev space
is an $L^{2}$ quotient of a (local) pre-Sobolev space. Since we are interested
in _function_ spaces, pre-Sobolev spaces are the right framework for us. $f$
on $O$ for which for every compact $K\subset O$ the $L^{2}$ weak derivative
$\partial^{\bm{\alpha}}_{w}f$ exists on $K$ for all multi-index
$\bm{\alpha}\in\mathbb{N}_{0}^{d}$ with $|\bm{\alpha}|\leq n$, i.e. for all
such $\bm{\alpha}$ there is a function $\partial^{\bm{\alpha}}_{w}f\in
L^{2}(K)$ such that

$\int_{K}\partial^{\bm{\alpha}}_{w}f(\bm{x})\varphi(\bm{x})\,d\bm{x}=(-1)^{|\bm{\alpha}|}\int_{K}f(\bm{x})\partial^{\bm{\alpha}}\varphi(\bm{x})\,d\bm{x}$

for all $\varphi\in C^{\infty}(K)$.

######  Theorem 24 (Sample path Sobolev regularity)

Let $n\in\mathbb{N}$. The process $f\sim\mathcal{GP}(0,k)$ has samples in
$\mathcal{H}^{n}_{loc}(O)$ if,

  1. (1)

for general $k$, $k\in C^{n\otimes n}(O\times O)$.

  2. (2)

for stationary $k(\bm{x},\bm{y})=k_{\delta}(\bm{x}-\bm{y})$,
$\partial^{\bm{\alpha}}k_{\delta}$ exists at $\bm{0}$ for all
$|\bm{\alpha}|\leq 2n$.

  3. (3)

for isotropic $k(\bm{x},\bm{y})=k_{r}(\|\bm{x}-\bm{y}\|)$, $k_{r}^{(j)}$
exists at $0$ for all $j\leq 2n$.

If we assume that $k$ is continuous, or more generally that $f$ has a
measurable modification, then the converse holds in the scenarios (2) and (3).

Proof  (1) follows from the general theorem for second order measurable random
fields (Scheuerer, 2010b, Theorem 1). (2) and (3) follow from its corollary
(Scheuerer, 2010b, Corollary 1). The converse in (2) follows from (Scheuerer,
2010b, Proposition 1). Finally, if (3) holds then, by the propagation of
regularity from Gneiting (1999), we have $k_{r}\in C^{2n}(\mathbb{R})$. So by
Lemma A.1 (3) $\Rightarrow$ (2) (to be precise an adapted version without the
Hölder conditions) we deduce (2), and hence we obtain the converse statement
in this case too.
As observed in Scheuerer (2010b), Sobolev differentiability is a natural
regularity notion for sample paths of general random fields as it can be
deduced from continuous mean square differentiability. In the case of GPs
however, if, in addition to the condition in Theorem 24 (1), we assume Hölder
control on the highest order partial derivatives of the kernel at the diagonal
as in Theorem 6 (1), we deduce strong differentiability of the sample paths.
The Hölder condition in Theorem 6 (1) being a weak one, one can usually not
obtain a greater number of weak derivatives using Theorem 24 than of strong
derivatives using Theorem 6. The reverse implications in Theorem 24 (2) and
Theorem 24 (3) show that, in the stationary and isotropic cases, this is not a
limitation of the theorem, but an inherent property of stationary and
isotropic GPs. For example, the proof of Proposition 9 Appendix B.1 shows
that, by the converse to Theorem 24 (3), even in the edge case
$\nu\in\mathbb{N}$, the samples of the Matérn GPs have as many Sobolev
derivatives as of strong derivatives.

Note however that one way to construct non-stationary GPs with more weak
derivatives than strong derivatives is through feature kernels, with features
that admit more weak derivatives than strong derivatives. This also shows that
it is not possible to obtain a converse for Theorem 24 (1) as we have for
Theorem 24 (2) and Theorem 24 (3).

Also note that deducing sample path Hölder regularity from applying Theorem 24
and Sobolev embedding theorems introduces a superfluous dependence on
dimension, and yields weaker results than the dimension independent Theorem 6.

We therefore expect Theorem 24 to have less practical value for GPs as its
general version (Scheuerer, 2010b, Theorem 1) for second order measurable
random fields. For GPs we encourage the use of Theorem 6.

Acknowledgments and Disclosure of Funding

NDC, MP and PH gratefully acknowledge financial support by the European
Research Council through ERC StG Action 757275 / PANAMA; the DFG Cluster of
Excellence “Machine Learning - New Perspectives for Science”, EXC 2064/1,
project number 390727645; the German Federal Ministry of Education and
Research (BMBF) through the Tübingen AI Center (FKZ: 01IS18039A); and funds
from the Ministry of Science, Research and Arts of the State of Baden-
Württemberg. NDC acknowledges the support of the Fonds National de la
Recherche, Luxembourg. The authors thank the International Max Planck Research
School for Intelligent Systems (IMPRS-IS) for supporting NDC and MP.

##  Appendix A Proof of Theorem 6

We start by showing the equivalences of the various kernel conditions in
Theorem 6, namely (1) $\Leftrightarrow$ (2) in the stationary case, and (1)
$\Leftrightarrow$ (2) $\Leftrightarrow$ (3) in the isotropic case.

###### Lemma A.1

Let $n\in\mathbb{N}_{0}$ and $\epsilon\in(0,1]$. Then the condition

  1. (1)

     * •

$k\in C^{n\otimes n}(O\times O)$,

     * •

$|\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x}+\bm{h},\bm{x}+\bm{h})-\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x}+\bm{h},\bm{x})-\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x},\bm{x}+\bm{h})+\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x},\bm{x})|=\mathcal{O}(\|\bm{h}\|^{2\epsilon})$
as $\bm{h}\to\bm{0}$, locally uniformly in $\bm{x}\in O$, for all
$|\bm{\alpha}|=|\bm{\beta}|=n$.

is equivalent to the following:

  1. (2)

for stationary $k(\bm{x},\bm{y})=k_{\delta}(\bm{x},\bm{y})$,

     * •

$k_{\delta}\in C^{2n}(\mathbb{R}^{d})$,

     * •

$|\partial^{\bm{\alpha}}k_{\delta}(\bm{h})-\partial^{\bm{\alpha}}k_{\delta}(\bm{0})|=\mathcal{O}(\|\bm{h}\|^{2\epsilon})$
as $\bm{h}\to\bm{0}$ for all $|\bm{\alpha}|=2n$.

  2. (3)

for isotropic $k(\bm{x},\bm{y})=k_{r}(\|\bm{x}-\bm{y}\|)$,

     * •

$k_{r}\in C^{2n}(\mathbb{R})$,

     * •

$|k_{r}^{(2n)}(h)-k_{r}^{(2n)}(0)|=\mathcal{O}(|h|^{2\epsilon})$ as $h\to 0$.

Proof of Lemma A.1 (1) $\Leftrightarrow$ (2) Let $k$ be stationary, i.e.
$k(\bm{x},\bm{y})=k_{\delta}(\bm{x}-\bm{y})$.

(2) $\Rightarrow$ (1): For $|\bm{\alpha}|,|\bm{\beta}|\leq n$,

$\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x},\bm{y})=(-1)^{|\bm{\beta}|}\partial^{\bm{\alpha}+\bm{\beta}}k_{\delta}(\bm{x}-\bm{y})$

for all $\bm{x},\bm{y}\in\mathbb{R}^{d}$, by the chain rule. So the existence
and the continuity of $\partial^{\bm{\alpha}+\bm{\beta}}k_{\delta}$ implies
the existence and the continuity of $\partial^{\bm{\alpha},\bm{\beta}}k$.
Moreover when $|\bm{\alpha}|=|\bm{\beta}|=n$,

$\displaystyle|\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x}+\bm{h},\bm{x}+\bm{h})-\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x}+\bm{h},\bm{x})-\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x},\bm{x}+\bm{h})+\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x},\bm{x})|$
(4)
$\displaystyle=2|\partial^{\bm{\alpha}+\bm{\beta}}k_{\delta}(\bm{0})-\partial^{\bm{\alpha}+\bm{\beta}}k_{\delta}(\bm{h})|$

for all $\bm{x},\bm{h}\in\mathbb{R}^{d}$, so the Hölder conditions on the
highest order partial derivatives of $k$ at the diagonal and of $k_{\delta}$
at $\bm{0}$ correspond.

(1) $\Rightarrow$ (2): For $|\bm{\alpha}|\leq 2n$, say
$\bm{\alpha}=\bm{\beta}+\bm{\gamma}$ with $|\bm{\beta}|,|\bm{\gamma}|\leq n$,
$\partial^{\bm{\alpha}}k_{\delta}(\bm{x})$ exists and is equal to
$(-1)^{|\bm{\gamma}|}\partial^{\bm{\beta},\bm{\gamma}}k(\bm{x},\bm{0})$, for
all $\bm{x}\in\mathbb{R}^{d}$. Moreover the Hölder conditions on the partial
derivatives of $k$ and $k_{\delta}$ correspond by Equation (4).
To prove the second equivalence in Lemma A.1, we need the following result:

###### Lemma A.2

Let $m\in\mathbb{N}_{0}$, $g\in C^{m}(\mathbb{R}\setminus\\{0\\})$, and
$\nu\in\mathbb{R}$ such that $g^{(j)}(h)=\mathcal{O}(|h|^{\nu-j})$ as $h\to
0$, for all $0\leq j\leq m$. Define
$f:=g\circ\|\cdot\|\colon\mathbb{R}^{d}\setminus\\{\bm{0}\\}\to\mathbb{R}$.
Then $f\in C^{m}(\mathbb{R}^{d}\setminus\\{\bm{0}\\})$ and
$\partial^{\bm{\alpha}}f(\bm{h})=\mathcal{O}(\|\bm{h}\|^{\nu-j})$ as
$\bm{h}\to\bm{0}$, for $0\leq j\leq m$ and $|\bm{\alpha}|=j$.

Proof  By induction on $m$. $m=0$ is clear; the growth/decay of $f$ at
$\bm{0}$ is the same as that of $g$ at $0$. For $m>0$,
$\bm{x}\in\mathbb{R}^{d}\setminus\\{\bm{0}\\}$, and $1\leq i\leq d$, note that

$\partial^{\bm{e}_{i}}f(\bm{x})=g^{\prime}(\|\bm{x}\|)\frac{x_{i}}{\|\bm{x}\|}=\tilde{g}(\|\bm{x}\|)x_{i}$
(5)

where $\bm{e}_{i}$ is the $i$th unit vector and
$\tilde{g}(x):=\frac{g^{\prime}(x)}{x}$ for $x\in\mathbb{R}\setminus\\{0\\}$.
Clearly $\tilde{g}\in C^{m-1}(\mathbb{R}\setminus\\{0\\})$ and

$\tilde{g}^{(j)}(h)=\sum_{l=0}^{j}\binom{j}{l}g^{(l+1)}(h)(-1)^{j-l}(j-l)!h^{-1-j+l}=\mathcal{O}(|h|^{(\nu-2)-j})$

as $h\to 0$, for all $0\leq j\leq m-1$. Let
$\tilde{f}=\tilde{g}\circ\|\cdot\|$. By the induction hypothesis,
$\partial^{\bm{\alpha}}\tilde{f}(\bm{h})=\mathcal{O}(\|\bm{h}\|^{(\nu-2)-j})$
as $\bm{h}\to\bm{0}$, for $|\bm{\alpha}|=j$ and $0\leq j\leq m-1$. Thus, by
(5),

$\partial^{\bm{\alpha}}\partial^{\bm{e}_{i}}f(\bm{h})=\left\\{\begin{aligned}
&h_{i}\partial^{\bm{\alpha}}\tilde{f}(\bm{h})&\ &\text{if }\alpha_{i}=0\\\
&h_{i}\partial^{\bm{\alpha}}\tilde{f}(\bm{h})+\partial^{\bm{\alpha}-\bm{e}_{i}}\tilde{f}(\bm{h})&&\text{if
}\alpha_{i}>0\end{aligned}\right\\}=\mathcal{O}(\|\bm{h}\|^{\nu-(j+1)})$

as $\bm{h}\to\bm{0}$. Consequently,
$\partial^{\bm{\alpha}}f(\bm{h})=\mathcal{O}(\|\bm{h}\|^{\nu-j})$ as
$\bm{h}\to\bm{0}$ for $|\bm{\alpha}|=j$ and $0\leq j\leq m$. This concludes
the induction step.

Proof of Lemma A.1 (1) $\Leftrightarrow$ (3) Let $k$ be isotropic, i.e.
$k(\bm{x},\bm{y})=k_{r}(\|\bm{x}-\bm{y}\|)$. By Lemma A.1 (1)
$\Leftrightarrow$ (2), it suffices to show that the condition on $k_{\delta}$
in (2) is equivalent to the condition on $k_{r}$ in (3).

(2) $\Rightarrow$ (3): This follows from the fact that $x\mapsto
k_{\delta}(x\bm{e}_{1})$ is exactly $k_{r}$, and analogously
$x\mapsto\partial^{j\bm{e}_{1}}k_{\delta}(x\bm{e}_{1})$ is equal to
$k_{r}^{(j)}$ for $0\leq j\leq 2n$.

(3) $\Rightarrow$ (2): $k_{r}\colon\mathbb{R}\to\mathbb{R}$ is even, so the
odd derivatives of $k_{r}$ which exist at 0 vanish there. Let

$g(x):=k_{r}(x)-\sum_{j=0}^{n}\frac{x^{2j}}{j!}k_{r}^{(2j)}(0),$

for all $x\in\mathbb{R}$. The Lagrange form of the remainder in an order
$2n-j-1$ Taylor expansion of $g^{(j)}$ then reveals that

$g^{(j)}(h)=\frac{h^{2n-j}}{(2n-j)!}(k_{r}^{(2n)}(\xi_{j}(h))-k_{r}^{(2n)}(0))=\mathcal{O}(|h|^{2n-j+2\epsilon})$

as $h\to 0$, for all $0\leq j\leq 2n$ and where $\xi_{j}(h)\in(0,h)$ (or
$(h,0)$ if $h<0$). So $g|_{\mathbb{R}\setminus\\{0\\}}$ satisfies the
conditions of Lemma A.2 with $m=2n$ and $\nu=2n+2\epsilon$. Let
$f:=g\circ\|\cdot\|\colon\mathbb{R}^{d}\to\mathbb{R}$. Then, by Lemma A.2,
$f|_{\mathbb{R}^{d}\setminus\\{\bm{0}\\}}\in
C^{2n}(\mathbb{R}^{d}\setminus\\{\bm{0}\\})$, and
$\partial^{\bm{\alpha}}f(\bm{h})=\mathcal{O}(\|\bm{h}\|^{2n+2\epsilon-j})$ as
$\bm{h}\to\bm{0}$ for $|\bm{\alpha}|=j$ and $0\leq j\leq 2n$. In particular,
$\partial^{\bm{\alpha}}f(\bm{h})\to 0$ as $\bm{h}\to\bm{0}$, and by the mean
value theorem this is sufficient to deduce that
$\partial^{\bm{\alpha}}f(\bm{0})$ exists and is 0. Now by noting that
$\|\bm{x}\|^{2j}=(x_{1}^{2}+\dots+x_{d}^{2})^{j}$ is smooth in $\bm{x}$ for
all $j$, we have that

$k_{\delta}(\bm{x})=f(\bm{x})+\sum_{j=0}^{n}\frac{\|\bm{x}\|^{2j}}{j!}k_{r}^{(2j)}(0)$

is in $C^{2n}(\mathbb{R})$. Moreover, for $\bm{h}\in\mathbb{R}^{d}$ and
$|\bm{\alpha}|=2n$ there is a constant $C_{\bm{\alpha}}\in\mathbb{R}$ such
that
$\partial^{\bm{\alpha}}k_{\delta}(\bm{h})=\partial^{\bm{\alpha}}f(\bm{h})+C_{\bm{\alpha}}$,
implying that

$|\partial^{\bm{\alpha}}k_{\delta}(\bm{h})-\partial^{\bm{\alpha}}k_{\delta}(\bm{0})|=|\partial^{\bm{\alpha}}f(\bm{h})-\underbrace{\partial^{\bm{\alpha}}f(\bm{0})}_{=0}|=\mathcal{O}(\|\bm{h}\|^{2n+2\epsilon-2n})=\mathcal{O}(\|\bm{h}\|^{2\epsilon})$

as $\bm{h}\to\bm{0}$. Thus $k_{\delta}$ satisfies the Hölder condition in (2).
Note that Lemma A.1 is purely a result about stationary and isotropic
functions; the proof does make use of the positive definiteness of $k$. In
fact, the positive definiteness of $k$ gives us additional information about
its regularity, which can be seen through Lemma A.3 below.

For a GP $f$, the mean-square partial derivative
$\partial^{\bm{e}_{i}}_{ms}f$, when it exists, is a stochastic process defined
on the same probability space such that $\partial^{\bm{e}_{i}}_{ms}f(\bm{x})$
is the $L^{2}$ limit of $\frac{f(\bm{x}+h\bm{e}_{i})-f(\bm{x})}{h}$ as $h\to
0$. The existence of a mean-square partial derivative is neither sufficient
nor necessary for the existence of sample derivatives. However such
derivatives are relevant for our work as they are directly related to the
derivatives of the covariance kernel:

###### Lemma A.3

Let $n\in\mathbb{N}_{0}$. For $f\sim\mathcal{GP}(0,k)$, the following
conditions are equivalent:

  1. (1)

$k\in C^{n\otimes n}(O\times O)$,

  2. (2)

$\partial^{\bm{\alpha},\bm{\alpha}}k$ exists in a neighbourhood the diagonal
and is continuous there for all $|\bm{\alpha}|\leq n$,

  3. (3)

the mean-square partial derivative $\partial^{\bm{\alpha}}_{ms}f$ exists and
is mean-square continuous for all $|\bm{\alpha}|\leq n$.

Moreover, if either condition holds then we have
$\partial^{\bm{\alpha}}_{ms}f\sim\mathcal{GP}(0,\partial^{\bm{\alpha},\bm{\alpha}}k)$
for all $|\bm{\alpha}|\leq n$.

Proof  We prove this for $n=1$; for general $n$ it then suffices to apply the
same argument inductively on the partial derivatives.

(1) $\Rightarrow$ (2): This implication is clear.

(2) $\Rightarrow$ (3): Let $1\leq i\leq d$. Note that

$\displaystyle\mathbb{E}\left[\left(\frac{f(\bm{x}+h\bm{e}_{i})-f(\bm{x})}{h}\right)\left(\frac{f(\bm{y}+h^{\prime}\bm{e}_{i})-f(\bm{y})}{h^{\prime}}\right)\right]$
(6)
$\displaystyle=\frac{\mathbb{E}[f(\bm{x}+h\bm{e}_{i})f(\bm{y}+h^{\prime}\bm{e}_{i})]-\mathbb{E}[f(\bm{x}+h\bm{e}_{i})f(\bm{y})]-\mathbb{E}[f(\bm{x})f(\bm{y}+h^{\prime}\bm{e}_{i})]+\mathbb{E}[f(\bm{x})f(\bm{y})]}{hh^{\prime}}$
$\displaystyle=\frac{k(\bm{x}+h\bm{e}_{i},\bm{y}+h^{\prime}\bm{e}_{i})-k(\bm{x}+h\bm{e}_{i},\bm{y})-k(\bm{x},\bm{y}+h^{\prime}\bm{e}_{i})+k(\bm{x},\bm{y})}{hh^{\prime}}$
$\displaystyle\to\partial^{\bm{e}_{i},\bm{e}_{i}}k(\bm{x},\bm{y})$

as $h,h^{\prime}\to 0$, for $\bm{x},\bm{y}\in O$ close enough, since
$\partial^{\bm{e}_{i},\bm{e}_{i}}k$ is assumed to be continuous in a
neighbourhood of the diagonal. Thus we have

$\displaystyle\mathbb{E}\left[\left(\frac{f(\bm{x}+h\bm{e}_{i})-f(\bm{x})}{h}-\frac{f(\bm{x}+h^{\prime}\bm{e}_{i})-f(\bm{x})}{h^{\prime}}\right)^{2}\right]$
(7)
$\displaystyle=\mathbb{E}\left[\left(\frac{f(\bm{x}+h\bm{e}_{i})-f(\bm{x})}{h}\right)^{2}\right]-2\mathbb{E}\left[\left(\frac{f(\bm{x}+h\bm{e}_{i})-f(\bm{x})}{h}\right)\left(\frac{f(\bm{x}+h^{\prime}\bm{e}_{i})-f(\bm{x})}{h^{\prime}}\right)\right]$
$\displaystyle\quad+\mathbb{E}\left[\left(\frac{f(\bm{x}+h^{\prime}\bm{e}_{i})-f(\bm{x})}{h^{\prime}}\right)^{2}\right]$
$\displaystyle\to\partial^{\bm{e}_{i},\bm{e}_{i}}k(\bm{x},\bm{x})-2\partial^{\bm{e}_{i},\bm{e}_{i}}k(\bm{x},\bm{x})+\partial^{\bm{e}_{i},\bm{e}_{i}}k(\bm{x},\bm{x})=0$

where the limit is taken as $h,h^{\prime}\to 0$. Equation (7) shows that
$\lim_{h\to 0}\frac{f(\bm{x}+h\bm{e}_{i})-f(\bm{x})}{h}$ exists in $L^{2}$,
i.e. that $\partial_{ms}^{\bm{e}_{i}}f(\bm{x})$ exists. Moreover Equation (6)
then implies that
$\mathbb{E}\left[\partial^{\bm{e}_{i}}_{ms}f(\bm{x})\partial^{\bm{e}_{i}}_{ms}f(\bm{y})\right]=\partial^{\bm{e}_{i},\bm{e}_{i}}k(\bm{x},\bm{y})$
for $\bm{x},\bm{y}\in O$ close enough. Hence

$\displaystyle\mathbb{E}\left[(\partial^{\bm{e}_{i}}_{ms}f(\bm{x}+\bm{h})-\partial^{\bm{e}_{i}}_{ms}f(\bm{x}))^{2}\right]$
$\displaystyle=\mathbb{E}\left[\partial^{\bm{e}_{i}}_{ms}f(\bm{x}+\bm{h})^{2}\right]-2\mathbb{E}\left[\partial^{\bm{e}_{i}}_{ms}f(\bm{x}+\bm{h})\partial^{\bm{e}_{i}}_{ms}f(\bm{x})\right]+\mathbb{E}\left[\partial^{\bm{e}_{i}}_{ms}f(\bm{x})^{2}\right]$
$\displaystyle=\partial^{\bm{e}_{i},\bm{e}_{i}}k(\bm{x}+\bm{h},\bm{x}+\bm{h})-2\partial^{\bm{e}_{i},\bm{e}_{i}}k(\bm{x}+\bm{h},\bm{x})+\partial^{\bm{e}_{i},\bm{e}_{i}}k(\bm{x},\bm{x})$
$\displaystyle\to 0$

as $\bm{h}\to\bm{0}$, so $\partial^{\bm{e}_{i}}_{ms}f$ is mean-square
continuous.

(3) $\Rightarrow$ (1): For $\bm{x},\bm{y}\in O$ and $1\leq i\leq d$ we have

$\frac{k(\bm{x}+h\bm{e}_{i},\bm{y})-k(\bm{x},\bm{y})}{h}=\mathbb{E}\left[\frac{f(\bm{x}+h\bm{e}_{i})-f(\bm{x})}{h}f(\bm{y})\right]\to\mathbb{E}\left[\partial^{\bm{e}_{i}}_{ms}f(\bm{x})f(\bm{y})\right]$

as $h\to 0$. Hence $\partial^{\bm{e}_{i},\bm{0}}k(\bm{x},\bm{y})$ exists.
Furthermore, for $1\leq j\leq d$,

$\displaystyle\frac{\partial^{\bm{e}_{i},\bm{0}}k(\bm{x},\bm{y}+h\bm{e}_{j})-\partial^{\bm{e}_{i},\bm{0}}k(\bm{x},\bm{y})}{h}$
$\displaystyle=\mathbb{E}\left[\partial^{\bm{e}_{i}}_{ms}f(\bm{x})\frac{f(\bm{y}+h\bm{e}_{j})-f(\bm{y})}{h}\right]$
(8)
$\displaystyle\to\mathbb{E}\left[\partial^{\bm{e}_{i}}_{ms}f(\bm{x})\partial^{\bm{e}_{j}}_{ms}f(\bm{y})\right]$

as $h\to 0$. Hence $\partial^{\bm{e}_{i},\bm{e}_{j}}k$ exists, and is
continuous by Equation (8), since $\partial^{\bm{e}_{i}}_{ms}f$ and
$\partial^{\bm{e}_{j}}_{ms}f$ are mean-square continuous.

Now it remains to show that in all cases
$\partial^{\bm{e}_{i}}_{ms}f\sim\mathcal{GP}(0,\partial^{\bm{e}_{i},\bm{e}_{i}}k)$.
Let $\\{\bm{x}_{1},\dots,\bm{x}_{N}\\}\subset O$ a finite set of points. Then
we need to show

$(\partial^{\bm{e}_{i}}_{ms}f(\bm{x}_{1}),\dots,\partial^{\bm{e}_{i}}_{ms}f(\bm{x}_{N}))\sim\mathcal{N}\left(\bm{0},\left(\partial^{\bm{e}_{i},\bm{e}_{i}}k(\bm{x}_{p},\bm{x}_{q})\right)_{p,q=1}^{N}\right).$
(9)

$(\partial^{\bm{e}_{i}}_{ms}f(\bm{x}_{1}),\dots,\partial^{\bm{e}_{i}}_{ms}f(\bm{x}_{N}))$
is the limit in $L^{2}(\Omega;\mathbb{R}^{N})$ as $h\to 0$ of

$\displaystyle\left(\frac{f(\bm{x}_{p}+h\bm{e}_{i})-f(\bm{x}_{p})}{h}\right)_{p=1}^{N}$
(10) $\displaystyle\sim\mathcal{N}\bigg{(}\bm{0},\bigg{(}$
$\displaystyle\frac{k(\bm{x}_{p}+h\bm{e}_{i},\bm{x}_{q}+h\bm{e}_{i})-k(\bm{x}_{p}+h\bm{e}_{i},\bm{x}_{q})-k(\bm{x}_{p},\bm{x}_{q}+h\bm{e}_{i})+k(\bm{x}_{p},\bm{x}_{q})}{h^{2}}\bigg{)}_{p,q=1}^{N}\bigg{)}.$

Now convergence in $L^{2}$ implies convergence in distribution, and we see
that the multivariate normal distribution in Equation (10) converges to the
multivariate normal distribution in Equation (9) as $h\to 0$. So this shows
$\partial^{\bm{e}_{i}}_{ms}f\sim\mathcal{GP}(0,\partial^{\bm{e}_{i},\bm{e}_{i}}k)$.
We can now prove our main result, Theorem 6, which we restate here for
convenience: *

Proof  We only prove the general case (1). Cases (2) and (3) then follow by
taking intersections over $\epsilon\in(0,\gamma)$ in Lemma A.1.

$\Longleftarrow$: For $n=0$ the result follows from applying the Kolmogorov
continuity theorem to GPs, see Azmoodeh et al. (2014, Theorem 1), or Potthoff
(2009, Page 347) for the result in arbitrary dimension $d$.

For $n=1$, the existence and continuity of $\partial^{\bm{e}_{i},\bm{e}_{i}}k$
implies by Lemma A.3 (1) $\Rightarrow$ (3) that $f$ is mean-square
differentiable in direction $i$. Moreover this mean-square partial derivative
satisfies
$\partial^{\bm{e}_{i}}_{ms}f\sim\mathcal{GP}(0,\partial^{\bm{e}_{i},\bm{e}_{i}}k)$.
Now applying the case $n=0$ to $\partial^{\bm{e}_{i},\bm{e}_{i}}k$ we deduce
that $\partial^{\bm{e}_{i}}_{ms}f$ has samples in222The mean-square partial
derivative $\partial^{\bm{e}_{i}}_{ms}f(\bm{x})$ at a point $\bm{x}\in O$ is
only well-defined almost surely on $\Omega$. But we still talk of sample paths
for the process $\partial^{\bm{e}_{i}}_{ms}f$, since we mean this up to
modification, as described in Section 2. $C^{\gamma^{-}}_{loc}(O)$. In
particular it has continuous samples, and this implies by (Potthoff, 2010,
Theorem 3.2) that $f$ has differentiable samples in direction $i$, with
$\partial^{\bm{e}_{i}}f(\bm{x})=\partial^{\bm{e}_{i}}_{ms}f(\bm{x})$ almost
surely, for all $\bm{x}\in O$. Thus
$\partial^{\bm{e}_{i}}f\sim\mathcal{GP}(0,\partial^{\bm{e}_{i},\bm{e}_{i}}k)$
and $\partial^{\bm{e}_{i}}f$ has samples in $C^{\gamma^{-}}_{loc}(O)$. $1\leq
i\leq d$ was arbitrary, so $f$ has samples in $C^{(1+\gamma)^{-}}_{loc}(O)$.

For $n>1$, we apply the same argument inductively on the partial derivatives.

$\Longrightarrow$: For $n=0$, Azmoodeh et al. (2014, Theorem 1) gives the
converse to Kolmogorov’s theorem for GPs. To be precise, this theorem is only
stated for dimension $d=1$, but we can see that the proof generalises to
arbitrary $d$.

Now suppose $n=1$. Pick $1\leq i\leq d$. For any finite set of points
$\\{\bm{x}_{1},\dots,\bm{x}_{N}\\}\subset O$,
$(\partial^{\bm{e}_{i}}f(\bm{x}_{1}),\dots,\partial^{\bm{e}_{i}}f(\bm{x}_{N}))$
is the almost sure limit of centered multivariate Gaussians with distribution

$\mathcal{N}\bigg{(}\bm{0},\bigg{(}\frac{k(\bm{x}_{p}+h\bm{e}_{i},\bm{x}_{q}+h\bm{e}_{i})-k(\bm{x}_{p}+h\bm{e}_{i},\bm{x}_{q})-k(\bm{x}_{p},\bm{x}_{q}+h\bm{e}_{i})+k(\bm{x}_{p},\bm{x}_{q})}{h^{2}}\bigg{)}_{p,q=1}^{N}\bigg{)}$

(see Equation (10)), so is itself a centered multivariate Gaussian
distribution. Hence $\partial^{\bm{e}_{i}}f$ is a GP.

Let $\bm{x},\bm{x}+h\bm{e}_{i}\in O$. By the mean value theorem, for each
$\omega\in\Omega$ there is $\xi_{\bm{x},\omega}(h)$ between $\bm{x}$ and
$\bm{x}+h\bm{e}_{i}$ such that

$\frac{f(\bm{x}+h\bm{e}_{i},\omega)-f(\bm{x},\omega)}{h}=\partial^{\bm{e}_{i}}f(\xi_{\bm{x},\omega}(h),\omega).$

Thus

$\left|\frac{f(\bm{x}+h\bm{e}_{i},\omega)-f(\bm{x},\omega)}{h}-\partial^{\bm{e}_{i}}f(\bm{x},\omega)\right|^{2}=|\partial^{\bm{e}_{i}}f(\xi_{\bm{x},\omega}(h),\omega)-\partial^{\bm{e}_{i}}f(\bm{x},\omega)|^{2}\leq
C_{\omega}^{2}|h|^{\gamma}$

since $\gamma\in(0,2\gamma)$, for some constant $C_{\omega}$ depending on
$\omega$ but not on $h$, assuming $h$ is small enough. Thus we have

$\mathbb{E}\left[\left(\frac{f(\bm{x}+h\bm{e}_{i})-f(\bm{x})}{h}-\partial^{\bm{e}_{i}}f(\bm{x})\right)^{2}\right]\leq\underbrace{\mathbb{E}\left[C_{\omega}^{2}\right]}_{(*)}\cdot\underbrace{|h|^{\gamma}}_{(**)}$

for $h$ small enough. $(*)<\infty$ by Azmoodeh et al. (2014, Theorem 1), which
we can apply since we showed that $\partial^{\bm{e}_{i}}f$ is a GP. Also
$(**)\to 0$ as $h\to 0$. Hence

$\frac{f(\bm{x}+h\bm{e}_{i})-f(\bm{x})}{h}\xrightarrow{L^{2}}\partial^{\bm{e}_{i}}f(\bm{x})$

as $h\to 0$, i.e. $f$ is mean-square differentiable in direction $i$. Hence by
Lemma A.3 (3) $\Rightarrow$ (1), $\partial^{\bm{e}_{i},\bm{e}_{j}}k$ exists
and is continuous for all $1\leq i,j\leq d$, and
$\partial^{\bm{e}_{i}}f\sim\mathcal{GP}(0,\partial^{\bm{e}_{i},\bm{e}_{i}}k)$.
Now

$\displaystyle\partial^{\bm{e}_{i},\bm{e}_{j}}k(\bm{x}+\bm{h},\bm{x}+\bm{h})-\partial^{\bm{e}_{i},\bm{e}_{j}}k(\bm{x}+\bm{h},\bm{x})-\partial^{\bm{e}_{i},\bm{e}_{j}}k(\bm{x},\bm{x}+\bm{h})+\partial^{\bm{e}_{i},\bm{e}_{j}}k(\bm{x},\bm{x})$
(11)
$\displaystyle=\mathbb{E}\left[(\partial^{\bm{e}_{i}}f(\bm{x}+\bm{h})-\partial^{\bm{e}_{i}}f(\bm{x}))(\partial^{\bm{e}_{j}}f(\bm{x}+\bm{h})-\partial^{\bm{e}_{j}}f(\bm{x}))\right]$
$\displaystyle\leq\mathbb{E}\left[(\partial^{\bm{e}_{i}}f(\bm{x}+\bm{h})-\partial^{\bm{e}_{i}}f(\bm{x}))^{2}\right]^{\nicefrac{{1}}{{2}}}\mathbb{E}\left[(\partial^{\bm{e}_{j}}f(\bm{x}+\bm{h})-\partial^{\bm{e}_{j}}f(\bm{x}))^{2}\right]^{\nicefrac{{1}}{{2}}}$
$\displaystyle\leq\underbrace{\mathbb{E}\left[A_{\omega}^{2}\right]^{\nicefrac{{1}}{{2}}}}_{<\infty}\|\bm{h}\|^{\epsilon}\cdot\underbrace{\mathbb{E}\left[B_{\omega}^{2}\right]^{\nicefrac{{1}}{{2}}}}_{<\infty}\|\bm{h}\|^{\epsilon}=\mathcal{O}(\|\bm{h}\|^{2\epsilon})$

as $\bm{h}\to\bm{0}$, locally uniformly in $\bm{x}\in O$, for all
$\bm{\epsilon}\in(0,\gamma)$, where $A_{\omega}$, $B_{\omega}$ are some
constants depending on $\omega$, by Azmoodeh et al. (2014, Theorem 1) applied
to $\partial^{\bm{e}_{i}}f$ and $\partial^{\bm{e}_{j}}f$, for $1\leq i,j\leq
d$.

Finally, for $n>1$ we apply the same argument inductively on the partial
derivatives.

##  Appendix B Proofs in Section 4

###  B.1 Proof of Proposition 9

Proof  For $\rho>0$ we can write for $\nu\not\in\mathbb{N}$

$K_{\nu}(\rho)=\frac{\pi}{2}\frac{I_{-\nu}(\rho)-I_{\nu}(\rho)}{\sin(\nu\pi)}$

where

$I_{\pm\nu}(\rho)=\left(\frac{\rho}{2}\right)^{\pm\nu}\sum_{j=0}^{\infty}\frac{1}{j!\,\Gamma(j\pm\nu+1)}\left(\frac{\rho}{2}\right)^{2j}$

are modified Bessel functions of the first kind (Abramowitz and Stegun, 1965,
Equations 9.6.2 & 9.6.10). So we can write the Matérn kernel (3) as

$k_{r}(x)=C_{\nu}\bigg{(}\underbrace{\sum_{j=0}^{\infty}\frac{1}{j!\,\Gamma(j-\nu+1)}\left(\frac{x}{2}\right)^{2j}}_{(*)}-\frac{|x|^{2\nu}}{2^{\nu}}\underbrace{\sum_{j=0}^{\infty}\frac{1}{j!\,\Gamma(j+\nu+1)}\left(\frac{x}{2}\right)^{2j}}_{(**)}\bigg{)}$

for $x\in\mathbb{R}$, where $C_{\nu}>0$ is some constant. The power series
$(*)$ and $(**)$ are smooth in $x$, and $|x|^{2\nu}$ has precisely $\lceil
2\nu-1\rceil$ continuous derivatives with $(2\nu-2\lceil\nu-1\rceil)$-Hölder
continuous $2\lceil\nu-1\rceil$th derivative at 0. So, for
$\nu\not\in\mathbb{N}$, the result follows by Theorem 6 (3).

For $\nu=n\in\mathbb{N}$, we instead use the following formula for the
modified Bessel function of the second kind (Abramowitz and Stegun, 1965,
Equations 9.6.11): for $\rho>0$

$\displaystyle K_{n}(\rho)$
$\displaystyle=\frac{1}{2}\left(\frac{\rho}{2}\right)^{-n}\sum_{j=0}^{n-1}\frac{(n-j-1)!}{j!}\left(-\frac{\rho^{2}}{4}\right)^{j}+(-1)^{n+1}\log\left(\frac{\rho}{2}\right)I_{n}(\rho)$
$\displaystyle+(-1)^{n}\frac{1}{2}\left(\frac{\rho}{2}\right)^{n}\sum_{j=0}^{\infty}\left(\psi(j+1)+\psi(n+j+1)\right)\frac{\left(\frac{\rho^{2}}{4}\right)^{j}}{j!(n+j)!}$

where $\psi$ is the digamma function. So we can write the Matérn kernel (3) as

$\displaystyle k_{r}(x)$
$\displaystyle=A_{n}\underbrace{\sum_{j=0}^{n-1}\frac{(n-j-1)!}{j!}(-1)^{j}\left(-\frac{x^{2}}{4}\right)^{j}}_{(I)}+B_{n}\underbrace{x^{2n}\log\left(\frac{x}{2}\right)}_{(II)}\underbrace{\sum_{j=0}^{\infty}\frac{1}{j!(n+j)!}\left(\frac{x}{2}\right)^{2j}}_{(III)}$
$\displaystyle+C_{n}\underbrace{x^{2n}\sum_{j=0}^{\infty}\left(\psi(j+1)+\psi(n+j+1)\right)\frac{\left(\frac{x}{2}\right)^{2j}}{j!(n+j)!}}_{(IV)}$

for some constants $A_{n},B_{n},C_{n}\neq 0$. $(I)$ and $(IV)$ are smooth in
$x$. $(III)$ is smooth in $x$ and is non-zero at $x=0$. $(II)$ has precisely
$2n-1$ continuous derivatives, with almost-$2$-Hölder continuous $(2n-2)$th
derivative at 0 (precisely, a Hölder decay of $\mathcal{O}(h^{2}\log h)$). So
the result follows by Theorem 6 (3).

###  B.2 Proof of Proposition 11

Proof  It is shown in (Wendland, 2004, Theorem 9.12) that the Wendland kernels
may be written as

$k_{r}(x)=\sum_{j=0}^{\lfloor\nicefrac{{d}}{{2}}\rfloor+3n+1}d^{(\lfloor\nicefrac{{d}}{{2}}\rfloor+n+1)}_{j,n}|x|^{j}$

for $x\in\mathbb{R}$, where
$d^{(\lfloor\nicefrac{{d}}{{2}}\rfloor+n+1)}_{j,n}\in\mathbb{R}$ are
coefficients. Furthermore the odd degree coefficients satisfy
$d^{(\lfloor\nicefrac{{d}}{{2}}\rfloor+n+1)}_{2j+1,n}=0$ if and only if $0\leq
j\leq 2n-1$. Now for all $j\in\mathbb{N}_{0}$, $|x|^{2j}$ is smooth and
$|x|^{2j+1}$ is precisely $2n$ times continuously differentiable with
Lipschitz $2n$th derivative. So the result follows by Theorem 6 (3).

###  B.3 Proof of Proposition 14

Proof  This follows from Theorem 6 (1), by noting that

$\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x},\bm{y})=\partial^{\bm{\alpha}}\bm{\phi}(\bm{x})^{T}\partial^{\bm{\beta}}\bm{\phi}(\bm{y})$

for all $\bm{x},\bm{y}\in O$ and $\bm{\alpha},\bm{\beta}\in\mathbb{N}_{0}^{d}$
with $|\bm{\alpha}|,|\bm{\beta}|\leq n$, and moreover when
$|\bm{\alpha}|=|\bm{\beta}|=n$,

$\displaystyle|\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x}+\bm{h},\bm{x}+\bm{h})-\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x}+\bm{h},\bm{x})-\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x},\bm{x}+\bm{h})+\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x},\bm{x})|$
$\displaystyle=(\partial^{\bm{\alpha}}\bm{\phi}(\bm{x}+\bm{h})-\partial^{\bm{\alpha}}\bm{\phi}(\bm{x}))^{T}(\partial^{\bm{\beta}}\bm{\phi}(\bm{x}+\bm{h})-\partial^{\bm{\beta}}\bm{\phi}(\bm{x}))$
$\displaystyle\leq\left(\sum_{i=1}^{m}|\partial^{\bm{\alpha}}\phi_{i}(\bm{x}+\bm{h})-\partial^{\bm{\alpha}}\phi_{i}(\bm{x})|^{2}\right)^{\nicefrac{{1}}{{2}}}\left(\sum_{i=1}^{m}|\partial^{\bm{\beta}}\phi_{i}(\bm{x}+\bm{h})-\partial^{\bm{\beta}}\phi_{i}(\bm{x})|^{2}\right)^{\nicefrac{{1}}{{2}}}$
$\displaystyle=\mathcal{O}(\|\bm{h}\|^{\bm{\epsilon}})\cdot\mathcal{O}(\|\bm{h}\|^{\bm{\epsilon}})=\mathcal{O}(\|\bm{h}\|^{2\bm{\epsilon}})$

as $\bm{h}\to\bm{0}$, locally uniformly in $\bm{x}\in O$, for all
$\epsilon\in(0,\gamma)$, where we used the Cauchy-Schwarz inequality.

###  B.4 Proof of Proposition 17

Proof  We assume $m=2$; the proof for general $m$ can then be done
inductively. Note that by the product rule

$\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x},\bm{y})=\sum_{\bm{\gamma}_{1}+\bm{\gamma}_{2}=\bm{\alpha}}\sum_{\bm{\delta}_{1}+\bm{\delta}_{2}=\bm{\beta}}\binom{|\bm{\alpha}|}{|\bm{\gamma}_{1}|}\binom{|\bm{\beta}|}{|\bm{\delta}_{1}|}\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x},\bm{y})\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x},\bm{y})$

for all $\bm{x},\bm{y}\in O$ and $\bm{\alpha},\bm{\beta}\in\mathbb{N}_{0}^{d}$
with $|\bm{\alpha}|,|\bm{\beta}|\leq n$.

Now for $|\bm{\alpha}|=|\bm{\beta}|=n$,

$\displaystyle\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x}+\bm{h},\bm{x}+\bm{h})-\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x}+\bm{h},\bm{x})-\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x},\bm{x}+\bm{h})+\partial^{\bm{\alpha},\bm{\beta}}k(\bm{x},\bm{x})$
(12)
$\displaystyle=\sum_{\bm{\gamma}_{1}+\bm{\gamma}_{2}=\bm{\alpha}}\sum_{\bm{\delta}_{1}+\bm{\delta}_{2}=\bm{\beta}}\binom{|\bm{\alpha}|}{|\bm{\gamma}_{1}|}\binom{|\bm{\beta}|}{|\bm{\delta}_{1}|}\Big{(}\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x}+\bm{h},\bm{x}+\bm{h})\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x}+\bm{h},\bm{x}+\bm{h})$
$\displaystyle\quad-\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x}+\bm{h},\bm{x})\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x}+\bm{h},\bm{x})-\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x},\bm{x}+\bm{h})\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x},\bm{x}+\bm{h})$
$\displaystyle\quad+\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x},\bm{x})\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x},\bm{x})\Big{)}.$

To show the expression in Equation (12) is
$\mathcal{O}(\|\bm{h}\|^{2\epsilon})$, it is therefore sufficient to show that

$\displaystyle|\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x}+\bm{h},\bm{x}+\bm{h})\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x}+\bm{h},\bm{x}+\bm{h})-\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x}+\bm{h},\bm{x})\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x}+\bm{h},\bm{x})$
(13)
$\displaystyle\quad-\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x},\bm{x}+\bm{h})\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x},\bm{x}+\bm{h})+\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x},\bm{x})\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x},\bm{x})|=\mathcal{O}(\|\bm{h}\|^{2\epsilon})$

as $\bm{h}\to\bm{0}$, locally uniformly in $\bm{x}\in O$, for all
$\epsilon\in(0,\gamma)$, $\bm{\gamma}_{1}+\bm{\gamma}_{2}=\bm{\alpha}$ and
$\bm{\delta}_{1}+\bm{\delta}_{2}=\bm{\beta}$. The inside of the absolute value
on the left hand side of Equation (13) can be written as

$\displaystyle\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x}+\bm{h},\bm{x}+\bm{h})\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x}+\bm{h},\bm{x}+\bm{h})-\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x}+\bm{h},\bm{x})\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x}+\bm{h},\bm{x})$
(14)
$\displaystyle\quad-\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x},\bm{x}+\bm{h})\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x},\bm{x}+\bm{h})+\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x},\bm{x})\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x},\bm{x})$
$\displaystyle=\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x}+\bm{h},\bm{x}+\bm{h})$
$\displaystyle\quad\cdot\underbrace{(\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x}+\bm{h},\bm{x}+\bm{h})-\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x}+\bm{h},\bm{x})-\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x},\bm{x}+\bm{h})+\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x},\bm{x}))}_{(*)}$
$\displaystyle\quad+\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x},\bm{x})$
$\displaystyle\quad\cdot\underbrace{(\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x}+\bm{h},\bm{x}+\bm{h})-\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x}+\bm{h},\bm{x})-\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x},\bm{x}+\bm{h})+\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x},\bm{x}))}_{(**)}$
$\displaystyle\quad+\underbrace{(\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x}+\bm{h},\bm{x})-\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x},\bm{x}))}_{(I)}\underbrace{(\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x}+\bm{h},\bm{x}+\bm{h})-\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x}+\bm{h},\bm{x}))}_{(II)}$
$\displaystyle\quad+\underbrace{(\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x},\bm{x}+\bm{h})-\partial^{\bm{\gamma}_{1},\bm{\delta}_{1}}k_{1}(\bm{x},\bm{x}))}_{(III)}\underbrace{(\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x}+\bm{h},\bm{x}+\bm{h})-\partial^{\bm{\gamma}_{2},\bm{\delta}_{2}}k_{2}(\bm{x},\bm{x}+\bm{h}))}_{(IV)}$

By arguing as in the proof of Theorem 6, Equation (11) we have
$(*),(**)=\mathcal{O}(\|\bm{h}\|^{2\epsilon})$. Moreover, by the Cauchy-
Schwarz inequality,

$\displaystyle|(I)|$
$\displaystyle=\left|\mathbb{E}\big{[}\partial^{\bm{\gamma}_{1}}f_{1}(\bm{x}+\bm{h})\partial^{\bm{\delta}_{1}}f_{1}(\bm{x})-\partial^{\bm{\gamma}_{1}}f_{1}(\bm{x})\partial^{\bm{\delta}_{1}}f_{1}(\bm{x})\big{]}\right|$
(15)
$\displaystyle\leq\mathbb{E}\big{[}\partial^{\bm{\delta}_{1}}f_{1}(\bm{x})^{2}\big{]}^{\nicefrac{{1}}{{2}}}\mathbb{E}\left[(\partial^{\bm{\gamma}_{1}}f_{1}(\bm{x}+\bm{h})-\partial^{\bm{\gamma}_{1}}f_{1}(\bm{x}))^{2}\right]^{\nicefrac{{1}}{{2}}}$
$\displaystyle=\partial^{\bm{\delta}_{1},\bm{\delta}_{1}}k_{1}(\bm{x},\bm{x})^{\nicefrac{{1}}{{2}}}(\partial^{\bm{\gamma}_{1},\bm{\gamma}_{1}}k_{1}(\bm{x}+\bm{h},\bm{x}+\bm{h})-\partial^{\bm{\gamma}_{1},\bm{\gamma}_{1}}k_{1}(\bm{x}+\bm{h},\bm{x})$
$\displaystyle\quad-\partial^{\bm{\gamma}_{1},\bm{\gamma}_{1}}k_{1}(\bm{x},\bm{x}+\bm{h})+\partial^{\bm{\gamma}_{1},\bm{\gamma}_{1}}k_{1}(\bm{x},\bm{x}))^{\nicefrac{{1}}{{2}}}$
$\displaystyle=\mathcal{O}(\|\bm{h}\|^{\epsilon})$

where $f_{1}\sim\mathcal{GP}(0,k_{1})$. Similarly
$(II),(III),(IV)=\mathcal{O}(\|\bm{h}\|^{\epsilon})$. By Equation (14) we
therefore deduce Equation (13), which concludes the proof.

## References

  * Abramowitz and Stegun (1965) Milton Abramowitz and Irene A. Stegun.  _Handbook of Mathematical Functions: With Formulas, Graphs, and Mathematical Tables_.  Courier Corporation, 1965.  ISBN 978-0-486-61272-0.
  * Adler and Taylor (2007) Robert Adler and Jonathan Taylor.  _Random Fields and Geometry_.  Springer Monographs in Mathematics. Springer, 2007.  ISBN 978-0-387-48116-6.  doi:10.1007/978-0-387-48116-6.
  * Azmoodeh et al. (2014) Ehsan Azmoodeh, Tommi Sottinen, Lauri Viitasaari, and Adil Yazigi.  Necessary and sufficient conditions for Hölder continuity of Gaussian processes.  _Statistics & Probability Letters_, 94(C):230–235, 2014.  doi:10.1016/j.spl.2014.07.030.
  * Berg et al. (1984) Christian Berg, Jens Peter Reus Christensen, and Paul Ressel.  _Harmonic Analysis on Semigroups_ , volume 100 of _Graduate Texts in Mathematics_.  Springer, 1984.  ISBN 978-1-4612-1128-0.  doi:10.1007/978-1-4612-1128-0.
  * Borovitskiy et al. (2020) Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky, and Marc Deisenroth.  Matérn Gaussian processes on Riemannian manifolds.  In _Advances in Neural Information Processing Systems_ , volume 33, pages 12426–12437. Curran Associates, Inc., 2020.
  * Chentsov (1956) N. N. Chentsov.  Weak convergence of stochastic processes whose trajectories have no discontinuities of the second kind and the “heuristic” approach to the Kolmogorov-Smirnov tests.  _Theory of Probability & Its Applications_, 1(1):140–144, 1956.  ISSN 0040-585X.  doi:10.1137/1101013.
  * Da Costa et al. (2023) Nathaël Da Costa, Cyrus Mostajeran, Juan-Pablo Ortega, and Salem Said.  Invariant kernels on Riemannian symmetric spaces: a harmonic-analytic approach.  2023\.  doi:10.48550/arXiv.2310.19270.
  * Fernique (1975) Xavier Fernique.  Régularité des trajectoires des fonctions aleatoires gaussiennes.  In _Ecole d’Eté de Probabilités de Saint-Flour IV—1974_ , Lecture Notes in Mathematics, pages 1–96. Springer, 1975.  ISBN 978-3-540-37600-2.  doi:10.1007/BFb0080190.
  * Friston et al. (2023) Karl Friston, Lancelot Da Costa, Noor Sajid, Conor Heins, Kai Ueltzhöffer, Grigorios A. Pavliotis, and Thomas Parr.  The free energy principle made simpler but not too simple.  _Physics Reports_ , 1024:1–29, 2023.  ISSN 0370-1573.  doi:10.1016/j.physrep.2023.07.001.
  * Gneiting (1999) Tilmann Gneiting.  On the derivatives of radial positive definite functions.  _Journal of Mathematical Analysis and Applications_ , 236(1):86–93, 1999.  ISSN 0022-247X.  doi:10.1006/jmaa.1999.6434.
  * Heins et al. (2023) Conor Heins, Beren Millidge, Lancelot Da Costa, Richard Mann, Karl Friston, and Iain Couzin.  Collective behavior from surprise minimization.  2023\.  doi:10.48550/arXiv.2307.14804.
  * Henderson (2024) Iain Henderson.  Sobolev regularity of Gaussian random fields.  _Journal of Functional Analysis_ , 286(3):110241, 2024.  ISSN 0022-1236.  doi:10.1016/j.jfa.2023.110241.
  * Kanagawa et al. (2018) Motonobu Kanagawa, Philipp Hennig, Dino Sejdinovic, and Bharath K. Sriperumbudur.  Gaussian processes and kernel methods: a review on connections and equivalences.  2018\.  doi:10.48550/arXiv.1807.02582.
  * Klenke (2014) Achim Klenke.  _Probability Theory: A Comprehensive Course_.  Universitext. Springer, 2014.  ISBN 978-1-4471-5361-0.  doi:10.1007/978-1-4471-5361-0.
  * Pförtner et al. (2022) Marvin Pförtner, Ingo Steinwart, Philipp Hennig, and Jonathan Wenger.  Physics-informed Gaussian process regression generalizes linear PDE solvers.  2022\.  doi:10.48550/arXiv.2212.12474.
  * Pierce (1970) Albert Pierce.  _Fundamentals of Nonparametric Statistics_.  Dickenson Publishing Company, 1970.  ISBN 978-0822143109.
  * Potthoff (2009) Jürgen Potthoff.  Sample properties of random fields. II. Continuity.  _Communications on Stochastic Analysis_ , 3(3), 2009.  ISSN 0973-9599.  doi:10.31390/cosa.3.3.02.
  * Potthoff (2010) Jürgen Potthoff.  Sample properties of random fields III: differentiability.  _Communications on Stochastic Analysis_ , 4(3), 2010.  ISSN 2688-6669.  doi:10.31390/cosa.4.3.03.
  * Rasmussen and Williams (2005) Carl Edward Rasmussen and Christopher K. I. Williams.  _Gaussian Processes for Machine Learning_.  The MIT Press, 2005.  ISBN 978-0-262-25683-4.  doi:10.7551/mitpress/3206.001.0001.
  * Scheuerer (2010a) Michael Scheuerer.  _A Comparison of Models and Methods for Spatial Interpolation in Statistics and Numerical Analysis_.  PhD thesis, 2010a.
  * Scheuerer (2010b) Michael Scheuerer.  Regularity of the sample paths of a general second order random field.  _Stochastic Processes and their Applications_ , 120(10):1879–1897, 2010b.  ISSN 0304-4149.  doi:10.1016/j.spa.2010.05.009.
  * Schober et al. (2014) Michael Schober, David K Duvenaud, and Philipp Hennig.  Probabilistic ODE solvers with Runge-Kutta means.  In _Advances in Neural Information Processing Systems_ , volume 27, 2014.
  * Talagrand (1987) Michel Talagrand.  Regularity of gaussian processes.  _Acta Mathematica_ , 159:99–149, 1987.  ISSN 0001-5962, 1871-2509.  doi:10.1007/BF02392556.
  * Wang et al. (2021) Junyang Wang, Jon Cockayne, Oksana Chkrebtii, T. J. Sullivan, and Chris. J. Oates.  Bayesian numerical methods for nonlinear partial differential equations.  _Statistics and Computing_ , 31(5):55, 2021.  ISSN 1573-1375.
  * Wendland (2004) Holger Wendland.  _Scattered Data Approximation_.  Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, Cambridge, 2004.  ISBN 978-0-521-84335-5.  doi:10.1017/CBO9780511617539.