# Posterior Contraction Rates for Matérn Gaussian Processes on Riemannian
Manifolds

Paul Rosa
University of Oxford&Viacheslav Borovitskiy
ETH ZürichAlexander Terenin
University of Cambridge
and Cornell University&Judith Rousseau
University of Oxford

###### Abstract

Gaussian processes are used in many machine learning applications that rely on
uncertainty quantification. Recently, computational tools for working with
these models in geometric settings, such as when inputs lie on a Riemannian
manifold, have been developed. This raises the question: can these intrinsic
models be shown theoretically to lead to better performance, compared to
simply embedding all relevant quantities into $\mathbb{R}^{d}$ and using the
restriction of an ordinary Euclidean Gaussian process? To study this, we prove
optimal contraction rates for intrinsic Matérn Gaussian processes defined on
compact Riemannian manifolds. We also prove analogous rates for extrinsic
processes using trace and extension theorems between manifold and ambient
Sobolev spaces: somewhat surprisingly, the rates obtained turn out to coincide
with those of the intrinsic processes, provided that their smoothness
parameters are matched appropriately. We illustrate these rates empirically on
a number of examples, which, mirroring prior work, show that intrinsic
processes can achieve better performance in practice. Therefore, our work
shows that finer-grained analyses are needed to distinguish between different
levels of data-efficiency of geometric Gaussian processes, particularly in
settings which involve small data set sizes and non-asymptotic behavior.

##  1 Introduction

Gaussian processes provide a powerful way to quantify uncertainty about
unknown regression functions via the formulation of Bayesian learning.
Motivated by applications in the physical and engineering sciences, a number
of recent papers [11, 9, 10, 37] have studied how to extend this model class
to spaces with geometric structure, in particular Riemannian manifolds
including important special cases such as spheres and Grassmannians [4],
hyperbolic spaces and spaces of positive definite matrices [5], as well as
general manifolds approximated numerically by a mesh [11].

These Riemannian Gaussian process models are starting to be applied for
statistical modeling, and decision-making settings such as Bayesian
optimization. For example, in a robotics setting, [27] has shown that using
Gaussian processes with the correct geometric structure allows one to learn
quantities such as the orientation of a robotic arm with less data compared to
baselines. The same model class has also been used by [15] to perform Gaussian
process regression on a manifold which models the geometry of a human heart
for downstream applications in medicine.

Given these promising empirical results, it is important to understand whether
these learning algorithms have good theoretical properties, as well as their
limitations. Within the Bayesian framework, a natural way to quantify data-
efficiency and generalization error is to posit a data-generating mechanism
model and study if—and how fast—the posterior distribution concentrates around
the true regression function as the number of observations goes to infinity.

Within the Riemannian setting, it is natural to compare _intrinsic_ methods,
which are formulated directly on the manifold of interest, with _extrinsic_
ones, which require one to embed the manifold within a higher-dimensional
Euclidean space. For example, the two-dimensional sphere can be embedded into
the Euclidean space $\mathbb{R}^{3}$: intrinsic Gaussian processes model
functions on the sphere while extrinsic ones model functions on
$\mathbb{R}^{3}$, which are then restricted to the sphere. Are the former more
efficient than the latter? Since embeddings—even isometric ones—at best only
preserve distances locally, they can induce spurious dependencies, as points
can be close in the ambient space but far away with respect to the intrinsic
geodesic distance: this is illustrated in Figure 1. In cases where embeddings
significantly alter distances, one can expect intrinsic models to perform
better, and it is therefore interesting to quantify such differences.

In other settings, the manifold on which the data lies can be unknown, which
makes using intrinsic methods directly no longer possible. There, one would
like to understand how well extrinsic methods can be expected to perform.
According to the _manifold hypothesis_ [18], it is common for perceptual data
such as text and images to concentrate on a lower-dimensional submanifold
within, for instance, pixel space or sequence space. It is therefore also
interesting to investigate how Gaussian process models—which, being kernel-
based, are simpler than for instance deep neural networks—perform in such
scenarios, at least in the asymptotic regime.

(a) Intrinsic

(b) Extrinsic

(c) Intrinsic

(d) Extrinsic

Figure 1:  Samples from different Matérn Gaussian processes on different
manifolds, namely a one-dimensional dumbbell-shaped manifold and a two-
dimensional sphere. Notice that the values across the dumbbell’s bottleneck
can be very different for the intrinsic process in (a), despite being very
close in the ambient Euclidean distance and in contrast to the situation for
the extrinsic model in (b). On the other hand, there is little qualitative
difference between (c) and (d), since the embedding produces a reasonably-good
global approximation to geodesic distances on the sphere.

In this work, we develop geometric analogs of the Gaussian process posterior
contraction theorems of [56]. More specifically, we derive posterior
contraction rates for three main geometric model classes: (1) the intrinsic
Riemannian Matérn Gaussian processes, (2) truncated versions of the intrinsic
Riemannian Matérn Gaussian processes, which are used in practice to avoid
infinite sums, and (3) the extrinsic Euclidean Matérn Gaussian processes under
the assumption that the data lies on a compact Riemannian manifold. In all
cases, we focus on IID randomly-sampled input points—commonly referred to as
_random design_ in the literature—and contraction in the sense of the
$L^{2}(p_{0})$ distance, defined in Section 2. We focus on _compact_
Riemannian manifolds: this allows one to define Matérn Gaussian processes
through their Karhunen–Loève expansions, which requires a discrete spectrum
for the Laplace-Beltrami operator—see for instance [11] and [13, Chapter
1]—and is a common setting in statistics [39].

#### Contributions.

We show that all three classes of Gaussian processes lead to optimal
procedures, in the minimax sense, as long as the smoothness parameter of the
kernel is aligned with the regularity of the unknown function. While this
result is natural—though non trivial—in the case of intrinsic Matérn
processes, it is rather remarkable that it also holds for extrinsic ones. This
means that in order to understand their differences better, finite-sample
considerations are necessary. We therefore present experiments that compute
the worst case errors numerically. These experiments highlight that intrinsic
models are capable of achieving better performance in the small-data regime.
We conclude with a discussion of why these results—which might at first seem
counterintuitive—are very natural when viewed from an appropriate mathematical
perspective: they suggest that optimality is perhaps best seen as a basic
property or an important guarantee that any sensible model should satisfy.

##  2 Background

_Gaussian process regression_ is a Bayesian approach to regression where the
modeling assumptions are $y_{i}=f(x_{i})+\varepsilon_{i}$, with
$\varepsilon_{i}\sim\operatorname{N}(0,\sigma_{\varepsilon}^{2})$, $x_{i}\in
X$, and $f$ is assigned a Gaussian process prior. A _Gaussian process_ is a
random function $f:X\rightarrow\mathbb{R}$ for which all finite-dimensional
marginal distributions are multivariate Gaussian. The distribution of such a
process is uniquely determined by its _mean function_
$m(\cdot)=\operatorname{\mathbb{E}}(f(\cdot))$ and _covariance kernel_
$k(\cdot,\cdot^{\prime})=\operatorname{Cov}(f(\cdot),f(\cdot^{\prime}))$,
hence we write $f\sim\operatorname{GP}(m,k)$.

For Gaussian process regression, the posterior distribution given the data is
also a Gaussian process with probability kernel
$\Pi(\cdot\mid\boldsymbol{x},\boldsymbol{y})=\operatorname{GP}(m_{\Pi(\cdot\mid\boldsymbol{x},\boldsymbol{y})},k_{\Pi(\cdot\mid\boldsymbol{x},\boldsymbol{y})})$,
see [41],

$\displaystyle
m_{\Pi(\cdot\mid\boldsymbol{x},\boldsymbol{y})}(\cdot)=\mathbf{K}_{(\cdot)\boldsymbol{x}}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}+\sigma_{\varepsilon}^{2}\mathbf{I})^{-1}\boldsymbol{y},$
(1) $\displaystyle
k_{\Pi(\cdot\mid\boldsymbol{x},\boldsymbol{y})}(\cdot,\cdot^{\prime})=\mathbf{K}_{(\cdot,\cdot^{\prime})}-\mathbf{K}_{(\cdot)\boldsymbol{x}}(\mathbf{K}_{\boldsymbol{x}\boldsymbol{x}}+\sigma_{\varepsilon}^{2}\mathbf{I})^{-1}\mathbf{K}_{\boldsymbol{x}(\cdot^{\prime})}.$
(2)

These quantities describe how incorporating data updates the information
contained within the Gaussian process. We will be interested studying the case
where $X$ is a Riemannian manifold, but first review the existing theory on
the asymptotic behaviour of the posterior when $X=[0,1]^{d}$.

###  2.1 Posterior Contraction Rates

Posterior contraction results describe how the posterior distribution
concentrates around the true data generating process, as the number of
observations increases, so that it eventually uncovers the true data-
generating mechanism. The area of _posterior asymptotics_ is concerned with
understanding conditions under which this does or does not occur, with
questions of _posterior contraction rates_ —how fast such convergence
occurs—being of key interest. At present, there is a well-developed literature
on posterior contraction rates, see [20] for a review.

In the context of Gaussian process regression with _random design_ , which is
the focus of this paper, the true data generating process is assumed to be of
the form

$\displaystyle y_{i}\mid x_{i}$
$\displaystyle\sim\operatorname{N}(f_{0}(x_{i}),\sigma_{\varepsilon}^{2})$
$\displaystyle x_{i}$ $\displaystyle\sim p_{0}$ (3)

where $f_{0}\in\mathcal{F}\subseteq\mathbb{R}^{X}$, a class of real-valued
functions, and $\mathrm{N}(\mu,\sigma^{2})$ denotes the Gaussian with moments
$\mu,\sigma^{2}$. Note that, in this particular variant, these equations
exactly mirror those of the Gaussian process model’s likelihood, including the
use of the same noise variance $\sigma_{\varepsilon}^{2}$ in both cases: in
this paper, we focus on the particular case where $\sigma_{\varepsilon}$ is
known in advance. This setting is restrictive, one can extend to an unknown
$\sigma_{\varepsilon}>0$ using techniques that are not specific to our
geometric setting: for instance, the approach of [55] allows to handle an
unknown $\sigma_{\varepsilon}$ if one assumes an upper and lower bound on it
and keep the same contraction rates. In practice, more general priors,
including ones that do not assume an upper or lower bound on
$\sigma_{\varepsilon}$, can be used, such as a conjugate one like in [6]—these
can also be analyzed to obtain contraction rates, albeit with additional
considerations. The generalization error for prediction in such models is
strongly related to the _weighted $L^{2}$ loss_ given by

$\left\lVert
f-f_{0}\right\rVert_{L^{2}\left(p_{0}\right)}=\left(\int_{X}\left|f(x)-f_{0}(x)\right|^{2}\mathop{}\\!\mathrm{d}p_{0}(x)\right)^{1/2}$
(4)

which is arguably the natural way of measuring discrepancy between $f$ and
$f_{0}$, given the fact that the covariates $x_{i}$ are sampled from $p_{0}$.
The posterior contraction rate is then defined as

$\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\operatorname{\mathbb{E}}_{f\sim\Pi(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left\lVert
f-f_{0}\right\rVert_{L^{2}(p_{0})}^{2}$ (5)

where
$\operatorname{\mathbb{E}}_{f\sim\Pi(\cdot\mid\boldsymbol{x},\boldsymbol{y})}(\cdot)$
denotes expectation under the posterior distribution while
$\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}(\cdot)$ denotes
expectation under the true data generating process.111Note that other notions
of posterior contraction can be found in the literature, see [20, 42] for
examples that are slightly weaker than the definition we work with. In the
case of covariates distributed on $[0,1]^{d}$, posterior contraction rates
have been derived under Matérn Gaussian process priors [47] in [56], who
showed the following result.

######  Result 1 (Theorem 2 of [56]).

In the Bayesian regression model, let $f$ be a mean-zero Matérn Gaussian
process prior on $\mathbb{R}^{d}$ with amplitude $\sigma^{2}_{f}$, length
scale $\kappa$, and smoothness $\nu>d/2$. Assume that the true data generating
process is given by (3), where $p_{0}$ has a Lebesgue density on $X=[0,1]^{d}$
which is bounded from below and above by $0<c_{p_{0}}<C_{p_{0}}<\infty$,
respectively. Let $f_{0}\in H^{\beta}\cap\mathcal{CH}^{\beta}$ with
$\beta>d/2$, where $H^{\beta}$ and $\mathcal{CH}^{\beta}$ the Sobolev and
Hölder spaces, respectively. Then there exists a constant $C>0$, which does
not depend on $n$ but does depend on $d$, $\sigma^{2}_{f}$, $\nu$, $\kappa$,
$\beta$, $p_{0}$, $\sigma_{\varepsilon}^{2}$, $\left\lVert
f_{0}\right\rVert_{H^{\beta}\left(\mathcal{M}\right)}$, and $\left\lVert
f_{0}\right\rVert_{\mathcal{CH}^{\beta}\left(\mathcal{M}\right)}$, such that

$\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\operatorname{\mathbb{E}}_{f\sim\Pi(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left\lVert
f-f_{0}\right\rVert_{L^{2}(p_{0})}^{2}\leq
Cn^{-\tfrac{2\min(\beta,\nu)}{2\nu+d}}$ (6)

and, moreover, the posterior mean satisfies

$\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\,\lVert
m_{\Pi(\cdot\mid\boldsymbol{x},\boldsymbol{y})}-f_{0}\rVert_{L^{2}(p_{0})}^{2}\leq
Cn^{-\tfrac{2\min(\beta,\nu)}{2\nu+d}}.$ (7)

Note that $m_{\Pi(\cdot\mid\boldsymbol{x},\boldsymbol{y})}$ is the Bayes
estimator [52] of $f$ associated to the weighted $L^{2}$ loss and that the
second inequality above is a direct consequence of the first. Therefore the
posterior contraction rate implies the same convergence rate for
$m_{\Pi(\cdot\mid\boldsymbol{x},\boldsymbol{y})}$. The best rate is attained
when $\beta=\nu$: that is, when true smoothness and prior smoothness
match—which is known to be minimax optimal in the problem of estimating
$f_{0}$: see [52]. In this paper, we extend this result to the manifold
setting.

###  2.2 Related Work and Current State of Affairs

The formalization of posterior contraction rates of Bayesian procedures dates
back to the work of [46] and [29], but has been extensively developed since
the seminal paper of [19] for various sampling and prior models, see for
instance [20, 42] for reviews. This includes, in particular, work on Gaussian
process priors [54, 56, 57, 43, 49]. Most of the results in the literature,
however, assume Euclidean data: as a consequence, contraction properties of
Bayesian models under manifold assumptions are still poorly understood, with
exception of some recent developments in both density estimation [7, 8, 60]
and regression [63, 60, 60].

The results closest to ours are those of [63, 12]. In the former, the authors
use an extrinsic length-scale-mixture of squared exponential Gaussian
processes to achieve optimal contraction rates with respect to the weighted
$L^{2}$ norm, using a completely different proof technique compared to us, and
their results are restricted to $f_{0}$ having Hölder smoothness of order less
than or equal to two. On the other hand [12] consider, as an intrinsic process
on the manifold, a hierarchical Gaussian process based on its heat kernel and
provide posterior contraction rates. For the Matérn class, [30] presents
results which characterize the asymptotic behavior of kernel hyperparameters:
our work complements these results by studying contraction of the Gaussian
process itself toward the unknown ground-truth function. One can also study
analogous discrete problems: [17, 45] present posterior contraction rates for
a specific graph Gaussian process model in a semi-supervised setting. In the
next section, we present our results on Matérn processes, defined either by
restriction of an ambient process or by an intrinsic construction, and discuss
their implications.

##  3 Posterior Contraction Rates on Compact Riemannian Manifolds

We now study posterior contraction rates for Matérn Gaussian processes on
manifolds, which are arguably the most-widely-used Gaussian process priors in
both the Euclidean and Riemannian settings. We begin by more precisely
describing our geometric setting before stating our key results and discussing
their implications. From now on, we write $X=\mathcal{M}$, to emphasize that
the covariate space is a manifold.

######  Assumption 2.

Assume that $\mathcal{M}\subset\mathbb{R}^{D}$ is a smooth, compact
submanifold (without boundary) of dimension $d<D$ equipped with the standard
Riemannian volume measure $\mu$.

We denote $|\mathcal{M}|=\int_{\mathcal{M}}d\mu(x)$ for volume of
$\mathcal{M}$. With this geometric setting defined, we will need to describe
regularity assumptions in terms of functional spaces on the manifold
$\mathcal{M}$. We work with Hölder spaces
$\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)$, defined using charts via the
usual Euclidean Hölder spaces, the Sobolev spaces
$H^{s}\left(\mathcal{M}\right)$, and Besov spaces
$B_{\infty,\infty}^{s}\left(\mathcal{M}\right)$ which are one of the ways of
generalizing the Euclidean Hölder spaces of smooth functions to manifolds. We
follow [14, 12], and define these spaces using the Laplace-Beltrami operator
on $\mathcal{M}$ in Appendix A.

Recall that the data-generating process is given by (3), with $f_{0}$ as the
true regression function and $p_{0}$ as the distribution of the covariates.

######  Assumption 3.

Assume that $p_{0}$ is absolutely continuous with respect to $\mu$, and that
its density, denoted by $p_{0}$, satisfies $c\leq p_{0}\leq C$ for
$0<c,C<\infty$. Assume the regression function $f_{0}$ satisfies $f_{0}\in
H^{\beta}\left(\mathcal{M}\right)\cap
B_{\infty,\infty}^{\beta}\left(\mathcal{M}\right)$ for some $\beta>d/2$, and
that $\sigma_{\varepsilon}^{2}>0$ is fixed and known.

This setting can be extended to handle unknown variance $\sigma_{\varepsilon}$
by putting a prior on $\sigma_{\varepsilon}$, following the strategy of [44,
36]. Since we are focused primarily on the impact of the manifold, we do not
pursue this here. With the setting fully defined, we proceed to develop
posterior contraction results for different types of Matérn Gaussian process
priors: intrinsic, intrinsic truncated and extrinsic.

(a) Intrinsic

(b) Extrinsic

(c) Intrinsic

(d) Extrinsic

Figure 2: Different Matérn kernels $k(\bullet,x)$ on different manifolds.

###  3.1 Intrinsic Matérn Gaussian Processes

We now introduce the first geometric Gaussian process prior under study—the
Riemannian Matérn kernel of [62, 33, 11]. This process was originally defined
using stochastic partial differential equations: here, we present it by its
Karhunen–Loève expansion, to facilitate comparisons with its truncated analogs
presented in Section 3.2.

######  Definition 4 (Intrinsic Matérn prior).

Let $\nu>0$, and let $(\lambda_{j},f_{j})_{j\geq 0}$ be the eigenvalues and
orthonormal eigenfunctions of the Laplace–Beltrami operator on $\mathcal{M}$,
in increasing order. Define the intrinsic _Riemannian Matérn Gaussian process_
through its Karhunen–Loève expansion to be

$\displaystyle f(\cdot)$
$\displaystyle=\frac{\sigma_{f}^{2}}{C_{\nu,\kappa}}\sum_{j=1}^{\infty}\left(\frac{2\nu}{\kappa^{2}}+\lambda_{j}\right)^{-\frac{\nu+d/2}{2}}\xi_{j}f_{j}(\cdot)$
$\displaystyle\xi_{j}$ $\displaystyle\sim\operatorname{N}\left(0,1\right)$ (8)

where $\nu,\kappa,\sigma_{f}^{2}$ are positive parameters and $C_{\nu,\kappa}$
is the normalization constant, chosen such that
$\frac{1}{|\mathcal{M}|}\int_{M}\operatorname{Var}\left(f(x)\right)\mathop{}\\!\mathrm{d}\mu(x)=\sigma_{f}^{2}$,
where $\operatorname{Var}$ denotes the variance.

The covariance kernels of these processes are visualized in Figure 2. With
this prior, and the setting defined in Section 3, we are ready to present our
first result: this model attains the desired optimal posterior contraction
rate as soon as the regularity of the ground-truth function matches the
regularity of the Gaussian process, as described by the parameter $\nu$.

theoremThmIntrinsic Let $f$ be a Riemannian Matérn Gaussian process prior of
Definition 4 with smoothness parameter $\nu>d/2$ and let $f_{0}$ satisfy 3.
Then there is a $C>0$ such that

$\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\operatorname{\mathbb{E}}_{f\sim\Pi(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left\lVert
f-f_{0}\right\rVert_{L^{2}(p_{0})}^{2}\leq
Cn^{-\tfrac{2\min(\beta,\nu)}{2\nu+d}}.$ (9)

All proofs are given in Appendix B. Our proof follows the general approach of
[56], by first proving a contraction rate with respect to the distance
$n^{-1/2}\left\lVert
f(\boldsymbol{x})-f_{0}(\boldsymbol{x})\right\rVert_{\mathbb{R}^{n}}$ at input
locations $\boldsymbol{x}$, and then extending the result to the true
$L^{2}$-distance by applying a suitable concentration inequality. The first
part is obtained by studying the _concentration function_ , which is known to
be the key quantity to control in order to derive contraction rates of
Gaussian process priors—see [20, 57] for an overview. Given our regularity
assumptions on $f_{0}$, the most difficult part lies in controlling the small-
ball probabilities $\Pi\big{[}\left\lVert
f\right\rVert_{\mathcal{C}(\mathcal{M})}<\varepsilon\big{]}$: we handle this
by using results relating this quantity with the entropy of an RKHS unit ball
with respect to the uniform norm. Since our process’ RKHS is related to the
Sobolev space $H^{\nu+d/2}\left(\mathcal{M}\right)$ which admits a description
in terms of charts, we apply results on the entropy of Sobolev balls in the
Euclidean space to conclude the first part. Finally, to extend the rate to the
true $L^{2}(p_{0})$ norm, following [56], we prove a Hölder-type property for
manifold Matérn processes, and apply Bernstein’s inequality. Together, this
gives the claim.

This result is good news for the intrinsic Matérn model: it tells us that
asymptotically it incorporates the data as efficiently as possible at least in
terms of posterior contraction rates, given that its regularity matches the
regularity of $f_{0}$. An inspection of the proof shows that the constant
$C>0$ can be seen to depend on $d$,$\sigma_{f}^{2}$, $\nu$, $\kappa$, $\beta$,
$p_{0}$,$\sigma_{\varepsilon}^{2}$, $\left\lVert
f_{0}\right\rVert_{H^{\beta}\left(\mathcal{M}\right)}$, $\left\lVert
f_{0}\right\rVert_{B_{\infty\infty}^{\beta}\left(\mathcal{M}\right)}$, and
$\left\lVert
f_{0}\right\rVert_{\mathcal{CH}^{\beta}\left(\mathcal{M}\right)}$. Definition
4 extends to the case where the norm is raised to any power $q>1$ rather than
the second power, with the right-hand side raised to the same power: see
Appendix B for details. We now consider variants of this prior that can be
implemented in practice.

###  3.2 Truncated Matérn Gaussian Processes

The Riemannian Matérn prior’s covariance kernel cannot in general be computed
exactly, since Definition 4 involves an infinite sum. Arguably the simplest
way to implement these processes numerically is to truncate the respective
infinite series in the Karhunen–Loève expansion by taking the first $J$ terms,
which is also optimal in an $L^{2}(\mathcal{M})$-sense.

Note that the truncated prior is a randomly-weighted finite sum of
Laplace–Beltrami eigenfunctions, which have different smoothness properties
compared to the original prior: the truncated prior takes its values in
$\mathcal{C}^{\infty}\left(\mathcal{M}\right)$ since the eigenfunctions of
$\mathcal{M}$ are smooth—see for instance [16]. Nevertheless, if the
truncation level is allowed to grow as the sample size increases, then the
regularity of the process degenerates and one gets a function with
essentially-finite regularity in the limit.

Truncated random basis expansions have been studied extensively in the
Bayesian literature in the Euclidean setting—see for instance [2, 64] or [20,
Chapter 11] for examples with priors based on wavelet expansions. It is known
that truncating the expansion at a high enough level usually allows one to
retain optimality. Instead of truncating deterministically, it is also
possible to put a prior on the truncation level and resort to MCMC
computations which would then select the optimal number of basis functions
adaptively, at the expense of a more computationally intensive method—this is
done, for instance, in [53] in the context of drift estimation for diffusion
processes. Random truncation has been proven to lead in many contexts to
adaptive posterior contraction rates, meaning that although the prior does not
depend on the smoothness $\beta$ of $f_{0}$, the posterior contraction rate—up
to possible $\ln n$ terms—is of order $n^{-\beta/(2\beta+d)}$: see for
instance [2, 43].

By analogy of the Euclidean case with its random Fourier feature
approximations [40], we can call the truncated version of Definition 4 the
_manifold Fourier feature_ model, for which we now present our result.

theoremThmTruncatedIntrinsic Let $f$ be a Riemannian Matérn Gaussian process
prior on $\mathcal{M}$ with smoothness parameter $\nu>d/2$, modified to
truncate the infinite sum to at least $J_{n}\geq
cn^{\frac{d\left(\min(1,\nu/\beta)\right)}{2\nu+d}}$ terms, and let $f_{0}$
satisfy 3. Then there is a $C>0$ such that

$\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\operatorname{\mathbb{E}}_{f\sim\Pi(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left\lVert
f-f_{0}\right\rVert_{L^{2}(p_{0})}^{2}\leq
Cn^{-\tfrac{2\min(\beta,\nu)}{2\nu+d}}.$ (10)

The proof is essentially-the-same as the non-truncated Matérn, but involves
tracking dependence of the inequalities on the truncation level $J_{n}$, which
implicitly defines a sequence of priors rather than a single fixed prior.

This result is excellent news for the intrinsic models: it means that they
inherit the optimality properties of the limiting one, even in the absence of
the infinite sum—in spite of the fact that the corresponding finite-truncation
prior places its probability on
$\mathcal{C}^{\infty}\left(\mathcal{M}\right)$. Again, the constant $C>0$ can
be seen to depend on $d$,$\sigma_{f}^{2}$,
$\nu$,$\kappa$,$\beta$,$p_{0}$,$\sigma_{\varepsilon}^{2}$, $\left\lVert
f_{0}\right\rVert_{H^{\beta}\left(\mathcal{M}\right)}$, $\left\lVert
f_{0}\right\rVert_{B_{\infty\infty}^{\beta}\left(\mathcal{M}\right)}$, and
$\left\lVert
f_{0}\right\rVert_{\mathcal{CH}^{\beta}\left(\mathcal{M}\right)}$. This
concludes our results for the intrinsic Riemannian Matérn priors. We now study
what happens if, instead of working with a geometrically-formulated model, we
simply embed everything into $\mathbb{R}^{d}$ and formulate our models there.

###  3.3 Extrinsic Matérn Gaussian Processes

The results of the preceding sections provide good reason to be excited about
the intrinsic Riemannian Matérn prior: the rates it obtains match the usual
minimax rates seen for the Euclidean Matérn prior and Euclidean data, provided
that we match the smoothness $\nu$ with the regularity of $f_{0}$. Another
possibility is to consider an extrinsic Gaussian process, that is, a Gaussian
process defined over an ambient space. This has been considered by [63] for
instance for the square-exponential process, in an adaptive setting where one
does not assume that the regularity $\beta$ of $f_{0}$ is explicitly known,
but where $\beta\leq 2$. In this section we prove a non-adaptive analog of
this result for the Matérn process.

######  Definition 5 (Extrinsic Matérn prior).

Assume that the manifold $\mathcal{M}$ is isometrically embedded in the
Euclidean space $\mathbb{R}^{D}$, such that we can regard $\mathcal{M}$ as a
subset of $\mathbb{R}^{D}$. Consider the Gaussian process with zero mean and
kernel given by restricting onto $\mathcal{M}$ the standard Euclidean Matérn
kernel

$k_{\nu,\kappa,\sigma_{f}^{2}}(x,x^{\prime})=\sigma_{f}^{2}\frac{2^{1-\nu}}{\Gamma(\nu)}\left(\sqrt{2\nu}\frac{\left\lVert
x-x^{\prime}\right\rVert_{\mathbb{R}^{D}}}{\kappa}\right)^{\nu}K_{\nu}\left(\sqrt{2\nu}\frac{\left\lVert
x-x^{\prime}\right\rVert_{\mathbb{R}^{D}}}{\kappa}\right)$ (11)

where $\sigma_{f},\kappa,\nu>0$ and $K_{\nu}$ is the modified Bessel function
of the second kind [22].

Since the extrinsic Matérn process is defined in a completely agnostic way
with respect to the manifold geometry, we would expect it to be less
performant when $\mathcal{M}$ is known. However, it turns out that the
extrinsic Matérn process converges at the same rate as the intrinsic one, as
given in the following claim.

theoremThmExtrinsic Let $f$ be a mean-zero extrinsic Matérn Gaussian process
prior with smoothness parameter $\nu>d/2$ on $\mathcal{M}$, and let $f_{0}$
satisfy 3. Then for some $C>0$ we have

$\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\operatorname{\mathbb{E}}_{f\sim\Pi(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left\lVert
f-f_{0}\right\rVert_{L^{2}(p_{0})}^{2}\leq
Cn^{-\tfrac{2\min(\beta,\nu)}{2\nu+d}}.$ (12)

Definition 5 is a surprising result because the optimal rates in this setting
only require the knowledge of the regularity $\beta$, but not the knowledge of
the manifold or the intrinsic dimension. More precisely, the prior is not
designed to be an adaptive prior, since it is a fixed Gaussian process, but it
surprisingly adapts to the dimension of the manifold, and thus to the
manifold.

The proof is also based on control of concentration functions. The main
difference is that, although the ambient process has a well known RKHS—the
Sobolev space $H^{s+D/2}\left(\mathbb{R}^{D}\right)$—the restricted process
has a non-explicit RKHS, which necessitates further analysis. We tackle this
issue by using results from [24] relating manifold and ambient Sobolev spaces
by linear bounded trace and extension operators, and from [63] describing a
general link between the RKHS of an ambient process and its restriction. This
allows us to show that the restricted process has an RKHS that is actually
norm-equivalent to the Sobolev space $H^{\nu+d/2}\left(\mathcal{M}\right)$,
which allows us to conclude the result in the same manner as in the intrinsic
case.

As consequence, our argument applies _mutatis mutandis_ in any setting where
suitable trace and extension theorems apply, with the Riemannian Matérn case
corresponding to the usual Sobolev results. In particular, our arguments
therefore apply directly to other processes possessing similar RKHSs, such as
for instance various kernels defined on the sphere—see e.g. [61, Chapter 17]
and [26]. The constant $C>0$ can be seen to depend on
$d$,$D$,$\sigma_{f}^{2}$,
$\nu$,$\kappa$,$\beta$,$p_{0}$,$\sigma_{\varepsilon}^{2}$, $\left\lVert
f_{0}\right\rVert_{H^{\beta}\left(\mathcal{M}\right)}$,$\left\lVert
f_{0}\right\rVert_{B_{\infty\infty}^{\beta}\left(\mathcal{M}\right)}$,$\left\lVert
f_{0}\right\rVert_{\mathcal{CH}^{\beta}\left(\mathcal{M}\right)}$—notice that
here $C$ depends implicitly on $D$ because of the presence of trace and
extension operator continuity constants. We now proceed to understand the
significance of the overall results.

###  3.4 Summary of Results

As a consequence of our previous results, fixing a single common data
generating distribution determined by $p_{0},f_{0}$, under suitable conditions
the intrinsic Matérn process, its truncated version, and the extrinsic Matérn
process all possess the _same_ posterior contraction rate with respect to the
$L^{2}\left(p_{0}\right)$-norm, which depends on $d$, $\nu$, and $\beta$, and
is optimal if the regularities of $f_{0}$ and the prior match. These results
imply the following immediate corollary, which follows by convexity of
$\left\lVert\cdot\right\rVert_{L^{2}(p_{0})}^{2}$ using Jensen’s inequality.

######  Corollary 6.

Under the assumptions of Definitions 4, 3.2 and 5, it follows that, for some
$C>0$

$\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\,\lVert
m_{\Pi(\cdot\mid\boldsymbol{x},\boldsymbol{y})}-f_{0}\rVert_{L^{2}(p_{0})}^{2}\leq
Cn^{-\frac{2\min(\beta,\nu)}{2\nu+d}}$ (13)

where $m_{\Pi(\cdot\mid\boldsymbol{x},\boldsymbol{y})}$ is the posterior mean
given a particular value of $\left(x_{i},y_{i}\right)_{i=1}^{n}$.

When $\nu=\beta$, the optimality of the rates we present in the manifold
setting can be easily inferred by lower bounding the $L^{2}$-risk of the
posterior mean by the $L^{2}$-risk over a small subset of $\mathcal{M}$ and
using charts, which translates the problem into the Euclidean framework for
which the rate is known to be optimal—see for instance [52].

To contextualize this, observe that even in cases where the geometry of the
manifold is non-flat, the asymptotic rates are unaffected by the choice of the
prior’s length scale $\kappa$—in either the intrinsic, or the extrinsic
case—but only by the smoothness parameter $\nu$. Indeed, the RKHS of the
process is only determined—up to norm equivalence—by $\nu$, which plays an
important role in the proofs. This, and the fact that extrinsic processes
attain the same rates, implies that the study of asymptotic posterior
contraction rates _cannot detect geometry_ in our setting, as was already
hinted by [63]. Hence, in the geometric setting, optimal posterior contraction
rates should be thought of more as a basic property that any reasonable model
should satisfy. Differences in performance will be down to constant factors
alone—but as we will see, these can be significant. To understand these
differences, we turn to empirical analysis.

##  4 Experiments

From Definitions 4, 3.2 and 5, we know that intrinsic and extrinsic Gaussian
processes exhibit the same posterior contraction rates in the asymptotic
regime. Here, we study how these rates manifest themselves in practice, by
examining how worst-case errors akin to those of Corollary 6 behave
numerically. Specifically, we consider the pointwise worst-case error

$v^{(\tau)}(t)=\sup_{\left\lVert f_{0}\right\rVert_{H^{\nu+d/2}}\leq
1}\operatorname{\mathbb{E}}_{\varepsilon_{i}\sim\operatorname{N}(0,\sigma_{\varepsilon}^{2})}\big{|}m_{\Pi(\cdot\mid\boldsymbol{x},\boldsymbol{y})}^{(\tau)}(t)-f_{0}(t)\big{|}^{2}$
(14)

where $m_{\Pi(\cdot\mid\boldsymbol{x},\boldsymbol{y})}^{(\tau)}$ is the
posterior mean corresponding to the zero-mean Matérn Gaussian process prior
with smoothness $\nu$, length scale $\kappa$, amplitude $\sigma_{f}^{2}$,
which is intrinsic if $\tau=\operatorname{i}$ or extrinsic if
$\tau=\operatorname{e}$. We use a Gaussian likelihood with noise variance
$\sigma_{\varepsilon}^{2}$ and observations
$y_{i}=f_{0}(x_{i})+\varepsilon_{i}$, and examine this quantity as a function
of the evaluation location $t\in\mathcal{M}$. By allowing us to assess how
error varies in different regions of the manifold, this provides us with a
fine-grained picture of how posterior contraction behaves.

One can show that $v^{(\tau)}$ may be computed without numerically solving an
infinite-dimensional optimization problem. Specifically, (14) can be
calculated, in the respective intrinsic and extrinsic cases, using

$\displaystyle v^{(\operatorname{i})}(t)$
$\displaystyle=k^{(\operatorname{i})}(t,t)-\mathbf{K}^{(\operatorname{i})}_{t\mathbf{X}}\left(\mathbf{K}^{(\operatorname{i})}_{\mathbf{X}\mathbf{X}}+\sigma_{\varepsilon}^{2}\mathbf{I}\right)^{-1}\mathbf{K}^{(\operatorname{i})}_{\mathbf{X}t}$
(15) $\displaystyle v^{(\operatorname{e})}(t)$
$\displaystyle\approx(\mathbf{K}_{t\,\mathbf{X}^{\prime}}^{(\operatorname{i})}-\boldsymbol{\alpha}_{t}\mathbf{K}_{\mathbf{X}\,\mathbf{X}^{\prime}}^{(\operatorname{i})})(\mathbf{K}_{\mathbf{X}^{\prime}\mathbf{X}^{\prime}}^{(\operatorname{i})})^{-1}(\mathbf{K}_{\mathbf{X}^{\prime}\,t}^{(\operatorname{i})}-\mathbf{K}_{\mathbf{X}^{\prime}\,\mathbf{X}}^{(\operatorname{i})}\boldsymbol{\alpha}_{t}^{\top})+\sigma_{\varepsilon}^{2}\boldsymbol{\alpha}_{t}\boldsymbol{\alpha}_{t}^{\top}$
(16)

where, for the extrinsic case,
$\boldsymbol{\alpha}_{t}=\mathbf{K}_{t\mathbf{X}}^{(\operatorname{i})}(\mathbf{K}_{\mathbf{X}\mathbf{X}}^{(\operatorname{e})}+\sigma_{\varepsilon}^{2}\mathbf{I})^{-1}$,
and $\mathbf{X}^{\prime}$ is a set of points sampled uniformly from the
manifold $\mathcal{M}$, the size of which determines approximation quality.
The intrinsic expression is simply the posterior variance
$k_{\Pi(\cdot\mid\boldsymbol{x},\boldsymbol{y})}^{(\operatorname{i})}(t,t)$,
and its connection with worst-case error is a well-known folklore result
mentioned somewhat implicitly in, for instance, [35]. The extrinsic expression
is very-closely-related, and arises by numerically approximating a certain
RKHS norm. A derivation of both is given in Appendix F. To assess the
approximation error of this formula, we also consider an analog of (16) but
instead defined for the intrinsic model, and compare it to (15): in all cases,
the difference between the exact and approximate expression was found to be
smaller than differences between models. By computing these expressions, we
therefore obtain, up to numerics, the pointwise worst-case expected error in
our regression model.

Figure 3:  Worst-case error estimates for the intrinsic and extrinsic
processes, on the _dumbbell_ , _sphere_ , and _dragon_ manifolds (lower is
better, $y$ axis is in the logarithmic scale). We see that, on the dumbbell
and dragon manifold, intrinsic models achieve lower expected errors than
extrinsic models for the ranges considered (top), and that their expected
error consistently varies less as a function of space (bottom). In contrast,
on the sphere, both models achieve similar performance, with differences
between models falling within the range of variability caused by different
random number seeds. We also see that the difference between computing the
pointwise worst-case error exactly and approximately, in the intrinsic case
where computing this difference is possible, is small in all cases.

For $\mathcal{M}$ we consider three settings: a dumbbell-shaped manifold, a
sphere, and the dragon manifold from the Stanford 3D scanning repository. In
all cases, we perform computations by approximating the manifold using a mesh,
and implementing the truncated Karhunen–Loève expansion with $J=500$
eigenpairs obtained from the mesh. We fix smoothness $\nu=\frac{5}{2}$,
amplitude $\sigma_{f}^{2}=1$, and noise variance
$\sigma_{\varepsilon}^{2}=0.0005$, for both the intrinsic and extrinsic Matérn
Gaussian processes. Since the interpretation of the length scale parameter is
manifold-specific, for the intrinsic Gaussian processes we set $\kappa=200$
for the dumbbell, $\kappa=0.25$ for the sphere, and $\kappa=0.05$ for the
dragon manifold. In all cases, this yielded functions that were neither close
to being globally-constant, nor resembled noise. Each experiment was repeated
$10$ times to assess variability. Complete experimental details are given in
Appendix G.222Code available at:
https://github.com/aterenin/geometric˙asymptotics.

The length scales $\kappa$ are defined differently for intrinsic and extrinsic
Matérn kernels: in particular, using the same length scale in both models can
result in kernels behaving very differently. To alleviate this, for the
extrinsic process, we set the length scale by maximizing the extrinsic
process’ marginal likelihood using the full dataset generated by the intrinsic
process, except in the dumbbell’s case where the full dataset is relatively
small, and therefore a larger set of 500 points was used instead. This allows
us to numerically match intrinsic and extrinsic length scales to ensure a
reasonably-fair comparison.

Figure 3 shows the mean, and spatial standard deviation of $v_{\tau}(t)$,
where by _spatial standard deviation_ we mean the sample standard deviation
computed with respect to locations in space, rather than with respect to
different randomly sampled datasets. From this, we see that on the dumbbell
and dragon manifold—whose geometry differs significantly from the respective
ambient Euclidean spaces—intrinsic models obtain better mean performance. The
standard deviation plot reveals that intrinsic models have errors that are
less-variable across space. This means that extrinsic models exhibit higher
errors in some regions rather than others—such as, for instance, regions where
embedded Euclidean and Riemannian distances differ—whereas in intrinsic models
the error decays in a more spatially-uniform manner.

In contrast, on the sphere, both models perform similarly. Moreover, both the
mean and spatial standard deviation decrease at approximately the same rates,
indicating that the extrinsic model’s predictions are correct about-as-often
as the intrinsic model’s, as a function of space. This confirms the view that,
since the sphere does not possess any bottleneck-like areas where embedded
Euclidean distances are extremely different from their Riemannian analogs, it
is significantly less affected by differences coming from embeddings.

In total, our experiments confirm that there are manifolds on which geometric
models can perform significantly better than non-geometric models. This
phenomenon was also noticed in [17], where a prior based on the
eigendecomposition of a random geometric graph, which can be thought as an
approximation of our intrinsic Matérn processes, is compared to a standard
extrinsic Gaussian process. In our experiments, we see this through expected
errors, mirroring prior results on Bayesian optimization performance. From our
theoretical results, such differences cannot be captured through posterior
contraction rates, and therefore would require sharper technical tools, such
as non-asymptotic analysis, to quantify theoretically.

##  5 Conclusion

In this work, we studied the asymptotic behavior of Gaussian process
regression with different classes of Matérn processes on Riemannian manifolds.
By using various results on Sobolev spaces on manifolds we derived posterior
contraction rates for intrinsic Matérn process defined via their Karhunen-
Loeve decomposition in the Laplace–Beltrami eigenbasis, including processes
arising from truncation of the respective sum which can be implemented in
practice. Next, using trace and extension theorems which relate manifold and
Euclidean Sobolev spaces, we derived similar contraction rates for the
restriction of an ambient Matérn process in the case where the manifold is
embedded in Euclidean space. These theoretical asymptotic results were
supplemented by experiments on several examples, showing significant
differences in performance between intrinsic and extrincic methods in the
small sample size regime when the manifold’s geometric structure differs from
the ambient Euclidean space. Our work therefore shows that capturing such
differences cannot be done through asymptotic contraction rates, motivating
and paving the way for further work on non-asymptotic error analysis to
capture empirically-observed differences between extrinsic and intrinsic
models.

## Acknowledgments

The authors are grateful to Mojmír Mutný and Prof. Andreas Krause for fruitful
discussions concerning this work. PR and JR were supported by the European
Research Council (ERC) under the European Union’s Horizon 2020 research and
innovation programme (grant agreement No. 834175). VB was supported by an ETH
Zürich Postdoctoral Fellowship. AT was supported by Cornell University,
jointly via the Center for Data Science for Enterprise and Society, the
College of Engineering, and the Ann S. Bowers College of Computing and
Information Science.

## References

  * [1] Charalambos D Aliprantis and Kim C Border  “Infinite Dimensional Analysis: A Hitchhiker’s Guide”  Springer, 2006
  * [2] Julyan Arbel, Ghislaine Gayraud and Judith Rousseau  “Bayesian Optimal Adaptive Estimation Using a Sieve Prior”  In _Scandinavian Journal of Statistics_ , 2013
  * [3] Thierry Aubin  “Some nonlinear problems in Riemannian geometry”  Springer, 1998
  * [4] Iskander Azangulov, Andrei Smolensky, Alexander Terenin and Viacheslav Borovitskiy  “Stationary Kernels and Gaussian Processes on Lie Groups and their Homogeneous Spaces I: the compact case”  In _arXiv:2208.14960_ , 2022
  * [5] Iskander Azangulov, Andrei Smolensky, Alexander Terenin and Viacheslav Borovitskiy  “Stationary Kernels and Gaussian Processes on Lie Groups and their Homogeneous Spaces II: non-compact symmetric spaces”  In _arXiv:2301.13088_ , 2023
  * [6] Sudipto Banerjee  “Modeling massive spatial datasets using a conjugate Bayesian linear modeling framework”  In _Spatial Statistics_ , 2020
  * [7] Clément Berenfeld, Paul Rosa and Judith Rousseau  “Estimating a density near an unknown manifold: a Bayesian nonparametric approach”  In _arXiv:2205.15717_ , 2022
  * [8] Abhishek Bhattacharya and David B. Dunson  “Nonparametric Bayesian density estimation on manifolds with applications to planar shapes.”  In _Biometrika_ , 2010
  * [9] Viacheslav Borovitskiy et al.  “Matérn Gaussian processes on graphs”  In _International Conference on Artificial Intelligence and Statistics_ , 2021
  * [10] Viacheslav Borovitskiy, Mohammad Reza Karimi, Vignesh Ram Somnath and Andreas Krause  “Isotropic Gaussian Processes on Finite Spaces of Graphs”  In _Artificial Intelligence and Statistics_ , 2023
  * [11] Viacheslav Borovitskiy, Alexander Terenin, Peter Mostowsky and Marc Peter Deisenroth  “Matérn Gaussian processes on Riemannian manifolds”  In _Advances in Neural Information Processing Systems_ , 2020
  * [12] Ismael Castillo, Gerard Kerkyacharian and Dominique Picard  “Thomas Bayes’ walk on manifolds”  In _Probability Theory and Related Fields_ , 2013
  * [13] Isaac Chavel  “Eigenvalues in Riemannian Geometry”  Elsevier, 1984
  * [14] Thierry Coulhon, Gerard Kerkyacharian and Pencho Petrushev  “Heat Kernel Generated frames in the setting of Dirichlet Spaces”  In _Journal of Fourier Analysis and Applications_ , 2012
  * [15] Sam Coveney et al.  “Gaussian Process Manifold Interpolation for Probabilistic Atrial Activation Maps and Uncertain Conduction Velocity”  In _Philosophical Transactions of the Royal Society A_ , 2020
  * [16] Ernesto De Vito, Nicole Mücke and Lorenzo Rosasco  “Reproducing kernel Hilbert spaces on manifolds: Sobolev and Diffusion spaces”  In _Analysis and Applications_ , 2021
  * [17] David B Dunson, Hau-Tieng Wu and Nan Wu  “Graph Based Gaussian Processes on Restricted Domains”  In _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 2021
  * [18] Charles Fefferman, Sanjoy Mitter and Hariharan Narayanan  “Testing the Manifold Hypothesis”  In _Journal of the American Mathematical Society_ , 2016
  * [19] Subhashis Ghosal, Jayanta K. Ghosh and Aad Vaart  “Convergence Rates of Posterior Distributions”  In _The Annals of Statistics_ , 2000
  * [20] Subhashis Ghosal and Aad Vaart  “Fundamentals of Nonparametric Bayesian Inference”, Cambridge Series in Statistical and Probabilistic Mathematics  Cambridge University Press, 2017
  * [21] Evarist Giné and Richard Nickl  “Mathematical Foundations of Infinite-Dimensional Statistical Models”, Cambridge Series in Statistical and Probabilistic Mathematics  Cambridge University Press, 2015
  * [22] Izrail Solomonovich Gradshteyn and Iosif Moiseevich Ryzhik  “Table of Integrals, Series, and Products”  Academic Press, 2014
  * [23] Alexander Grigoryan  “Heat Kernel and Analysis on Manifolds”  American Mathematical Society, 2009
  * [24] Nadine Große and Cornelia Schneider  “Sobolev spaces on Riemannian manifolds with bounded geometry: General coordinates and traces”  In _Mathematische Nachrichten_ , 2013
  * [25] Emmanuel Hebey  “Sobolev spaces on Riemannian manifolds”  Springer, 1996
  * [26] Simon Hubbert, Emilio Porcu, Chris J. Oates and Mark Girolami  “Sobolev Spaces, Kernels and Discrepancies over Hyperspheres”  In _Transactions on Machine Learning Research_ , 2023
  * [27] Noémie Jaquier et al.  “Geometry-aware bayesian optimization in robotics using Riemannian Matérn kernels”  In _Conference on Robot Learning_ , 2022
  * [28] Motonobu Kanagawa, Philipp Hennig, Dino Sejdinovic and Bharath K Sriperumbudur  “Gaussian processes and kernel methods: A review on connections and equivalences”  In _arXiv:1807.02582_ , 2018
  * [29] Lucien Le Cam  “Convergence of Estimates Under Dimensionality Restrictions”  In _The Annals of Statistics_ , 1973
  * [30] Didong Li, Wenpin Tang and Sudipto Banerjee  “Inference for Gaussian Processes with Matern Covariogram on Compact Riemannian Manifolds”  In _Journal of Machine Learning Research_ , 2023
  * [31] Wenbo V. Li and Werner Linde  “Approximation, Metric Entropy and Small Ball Estimates for Gaussian Measures”  In _The Annals of Probability_ , 1999
  * [32] Mikhail Lifshits  “Lectures on Gaussian processes”  Springer, 2012
  * [33] Finn Lindgren, Håvard Rue and Johan Lindström  “An explicit link between Gaussian fields and Gaussian Markov random fields: the stochastic partial differential equation approach”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 2011
  * [34] Sergey V Lototsky and Boris L Rozovsky  “Stochastic partial differential equations”  Springer, 2017
  * [35] Mojmir Mutny and Andreas Krause  “Experimental design for linear functionals in reproducing kernel hilbert spaces”  In _Advances in Neural Information Processing Systems_ , 2022
  * [36] Zacharie Naulet and Éric Barat  “Some Aspects of Symmetric Gamma Process Mixtures”  In _Bayesian Analysis_ , 2018
  * [37] Mu Niu et al.  “Intrinsic Gaussian processes on complex constrained domains”  In _Journal of the Royal Statistical Society: Series B (Statistical Methodology)_ , 2018
  * [38] Thomas Pinder and Daniel Dodd  “GPJax: A Gaussian Process Framework in JAX”  In _Journal of Open Source Software_ , 2022
  * [39] Emilio Porcu, Moreno Bevilacqua, Robert Schaback and Chris J Oates  “The Matérn Model: A Journey through Statistics, Numerical Analysis and Machine Learning”  In _arXiv:2303.02759_ , 2023
  * [40] Ali Rahimi and Benjamin Recht  “Random Features for Large-Scale Kernel Machines”  In _Advances in Neural Information Processing Systems_ , 2007
  * [41] Carl Edward Rasmussen and Christopher K.. Williams  “Gaussian Processes for Machine Learning”  MIT Press, 2006
  * [42] Judith Rousseau  “On the frequentist properties of Bayesian nonparametric methods”  In _Annual Review of Statistics and its Application_ , 2016
  * [43] Judith Rousseau and Botond Szabo  “Asymptotic behaviour of the empirical Bayes posteriors associated to maximum marginal likelihood estimator”  In _The Annals of Statistics_ , 2017
  * [44] Jean-Bernard Salomond  “Testing Un-Separated Hypotheses by Estimating a Distance”  In _Bayesian Analysis_ , 2018
  * [45] Daniel Sanz-Alonso and Ruiyi Yang  “Unlabeled Data Help in Graph-Based Semi-Supervised Learning: A Bayesian Nonparametrics Perspective”  In _Journal of Machine Learning Research_ , 2022
  * [46] Lorraine Schwartz  “On Bayes procedures”  In _Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete_ , 1965
  * [47] Michael L. Stein  “Interpolation of Spatial Data: Some Theory for Kriging”  Springer, 1999
  * [48] Robert S Strichartz  “Analysis of the Laplacian on the complete Riemannian manifold”  In _Journal of functional analysis_ 52.1, 1983, pp. 48–79
  * [49] Wenping Tang, Lu Zhang and Sudipto Banerjee  “On identifiability and consistency of the nugget in Gaussian spatial process models”  In _Journal of the Royal Statistical Society Series B: Statistical Methodology_ , 2021
  * [50] Hans Triebel  “Theory of Function Spaces I”  Birkhäuser, 1983
  * [51] Hans Triebel  “Theory of Function Spaces II”  Birkhäuser, 1992
  * [52] Alexandre B Tsybakov  “Introduction to Nonparametric Estimation”  Springer, 2008
  * [53] Frank Meulen, Moritz Schauer and Harry Zanten  “Reversible jump MCMC for nonparametric drift estimation for diffusion processes”  In _Computational Statistics and Data Analysis_ , 2014
  * [54] Aad Vaart and Harry Zanten  “Adaptive Bayesian estimation using a Gaussian random field with inverse Gamma bandwidth”  In _The Annals of Statistics_ , 2009
  * [55] Aad Vaart and Harry Zanten  “Bayesian inference with rescaled Gaussian process priors”  In _Electronic Journal of Statistics_ , 2007
  * [56] Aad Vaart and Harry Zanten  “Information Rates of Nonparametric Gaussian Process Methods”  In _Journal of Machine Learning Research_ , 2011
  * [57] Aad Vaart and Harry Zanten  “Rates of contraction of posterior distributions based on Gaussian process priors”  In _The Annals of Statistics_ , 2008
  * [58] Aad Vaart and Jon Wellner  “Weak Convergence and Empirical Processes”, Springer Series in Statistics  Springer Verlag, 1996
  * [59] Harry Zanten and Aad Vaart  “Reproducing kernel Hilbert spaces of Gaussian priors”  In _Institute of Mathematical Statistics Collections_ , 2008
  * [60] Xu Wang and Gilad Lerman  “Nonparametric Bayesian regression on manifolds via Brownian motion”  In _arXiv:1507.06710_ , 2015
  * [61] Holger Wendland  “Scattered Data Approximation”  Cambridge University Press, 2004
  * [62] Peter Whittle  “Stochastic Processes in Several Dimensions”  In _Bulletin of the International Statistical Institute_ , 1963
  * [63] Yun Yang and David B. Dunson  “Bayesian manifold regression”  In _The Annals of Statistics_ , 2016
  * [64] William Weimin Yoo, Vincent Rivoirard and Judith Rousseau  “Adaptive supremum norm posterior contraction: Wavelet spike-and-slab and anisotropic Besov spaces”  In _arXiv:1708.01909_ , 2017

##  Appendix A Preliminaries

Here we describe preliminaries necessary for Appendices B, C, D and E. This
includes some basic properties of the Laplace–Beltrami operator on compact
manifolds, partitions of unity subordinate to atlases, function spaces such as
Hölder, Sobolev and Besov spaces, general Gaussian random elements on Banach
spaces, and a certain technical lemma. Hereinafter, the expression $a\lesssim
b$ means $a\leq Cb$ for some constant $C>0$ whose value is irrelevant to our
claims.

###  A.1 Laplace–Beltrami Operator and Subordinate Partitions of Unity

Let $\mathcal{M}$ denote a compact connected Riemannian manifold without
boundary of dimension $d\in\mathbb{Z}_{>0}$. The Laplace–Beltrami operator
$\Delta$ on $\mathcal{M}$ is self-adjoint and positive semi-definite [48,
Theorem 2.4]. Let
$(L^{2}(\mathcal{M}),\left\langle\cdot,\cdot\right\rangle_{L^{2}(\mathcal{M})})$
denote the Hilbert space of square integrable (equivalence classes of)
functions on $\mathcal{M}$ with respect to the standard Riemannian volume
measure.

By standard theory [13, 23], there exists an orthonormal basis
$\left\\{f_{j}\right\\}_{j=0}^{\infty}$ of $L^{2}(\mathcal{M})$ consisting of
the eigenfunctions of $\Delta$, namely $\Delta f_{j}=-\lambda_{j}f_{j}$, with
$\lambda_{j}\geq 0$. We assume that the pairs $(\lambda_{j},f_{j})$ are sorted
such that $0=\lambda_{0}\leq\lambda_{j}\leq\lambda_{j+1}$. The growth of
$\lambda_{j}$ can be characterized as follows.

######  Result 7 (Weyl’s Law).

There exists a constant $C>0$ such that for all $j$ large enough we have

$C^{-1}j^{2/d}\leq\lambda_{j}\leq Cj^{2/d}.$ (17)

###### Proof.

See [13, Chapter 1]. ∎

Following [16, 24] we fix a family
$\mathcal{T}=\left(\mathcal{U}_{l},\phi_{l},\chi_{l}\right)_{l=1}^{L}$, where
$L\in\mathbb{Z}_{>0}$, the local coordinates
$\phi_{l}:\mathcal{U}_{l}\subset\mathcal{M}\to\mathcal{V}_{l}=\phi_{l}\left(\mathcal{U}_{l}\right)\subset\mathbb{R}^{d}$
are smooth diffeomorphisms, and the functions $\chi_{l}$ form a partition of
unity subordinate to $\left\\{\mathcal{U}_{l}\right\\}_{l=1}^{L}$, that is
$\chi_{l}\in\mathcal{C}^{\infty}\left(\mathcal{M}\right)$,
$\operatorname{supp}\left(\chi_{l}\right)\subset\mathcal{U}_{l}$,
$0\leq\chi_{l}\leq 1$ and $\sum_{l}\chi_{l}=1$—here, we can choose $L$ finite
by compactness of $\mathcal{M}$. For convenience and without loss of
generality, we assume that $\mathcal{V}_{l}=(0,1)^{d}$. With this, we can
start defining function spaces on $\mathcal{M}$.

###  A.2 Hölder Spaces $\mathcal{CH}^{\gamma}$ and the spaces of continuous
functions $\mathcal{C}$ and smooth functions $\mathcal{C}^{\infty}$

Consider an arbitrary domain $\mathcal{X}\subseteq\mathbb{R}^{d}$ or
$\mathcal{X}\subseteq\mathcal{M}$. We denote the class of infinitely
differentiable functions on $\mathcal{X}$ by
$\mathcal{C}^{\infty}(\mathcal{X})$. Let $\mathcal{C}^{k}(\mathcal{X})$ denote
the Banach space of $k\in\mathbb{Z}_{>0}$ times continuously differentiable
functions on $\mathcal{X}$ with finite norm

$\left\lVert
f\right\rVert_{\mathcal{C}^{k}(\mathcal{X})}=\sup_{x\in\mathcal{X}}\left|\nabla^{k}f(x)\right|$
(18)

where $\nabla^{k}$ is the $k$th covariant derivative, as in [25, Section 2.1]
and [3, Definition 2.2]. We also write
$\mathcal{C}(\mathcal{X})=\mathcal{C}^{0}(\mathcal{X})$ for the space of
continuous functions on $\mathcal{X}$.

The Euclidean Hölder spaces $\mathcal{CH}^{\gamma}(\mathbb{R}^{d})$, where
$\gamma=k+\alpha$ with $k\in\mathbb{Z}_{\geq 0}$ and $0<\alpha\leq 1$, are
defined by333Hölder spaces are often also denoted by $\mathcal{C}^{\gamma}$ as
well, with $\gamma\in\mathbb{R}_{>0}$. Since, using this formulation, they do
_not_ coincide with $\mathcal{C}^{k}(\mathcal{X})$ when
$k=\gamma\in\mathbb{Z}_{\geq 0}$, we use the notation $\mathcal{CH}^{\gamma}$
to avoid confusion.

$\displaystyle\mathcal{CH}^{\gamma}(\mathbb{R}^{d})$
$\displaystyle=\left\\{f\in\mathcal{C}^{k}(\mathbb{R}^{d}):\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}(\mathbb{R}^{d})}<\infty\right\\}$ (19)

where

$\displaystyle\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}(\mathbb{R}^{d})}$
$\displaystyle=\left\lVert
f\right\rVert_{\mathcal{C}^{k}(\mathbb{R}^{d})}+\sup_{\boldsymbol{x},\boldsymbol{y}\in\mathbb{R}^{d},\,\boldsymbol{x}\not=\boldsymbol{y}}\frac{\left|f(\boldsymbol{x})-f(\boldsymbol{y})\right|}{\left\lVert\boldsymbol{x}-\boldsymbol{y}\right\rVert_{\mathbb{R}^{d}}^{\alpha}}.$
(20)

More information on these definitions may be found, for instance, in [21, 51].
We now turn to the manifold versions of the Hölder spaces.

######  Definition 8 (Hölder spaces).

For all $\gamma>0$ we define the Hölder space
$\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)$ on the manifold $\mathcal{M}$
to be the space of all $f:\mathcal{M}\to\mathbb{R}$ satisfying

$\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}=\sum_{l=1}^{L}\left\lVert\left(\chi_{l}f\right)\circ\phi_{l}^{-1}\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathbb{R}^{d}\right)}<\infty.$
(21)

Since the charts $\phi_{l}$ are smooth, Definition 8 can be easily seen to be
independent of the chosen atlas, with equivalence of norms.

###  A.3 Sobolev and Besov Spaces

We now introduce the manifold versions of the Sobolev and Besov spaces, whose
definitions in the standard Euclidean case may be found, for instance, in
[51]. For Sobolev spaces we use the Bessel-potential-based definition,
following [16].

######  Definition 9 (Sobolev spaces).

For any $s>0$ we define the Sobolev space $H^{s}\left(\mathcal{M}\right)$ on
the manifold $\mathcal{M}$ as the Hilbert space of functions $f\in
L^{2}\left(\mathcal{M}\right)$ such that $\left\lVert
f\right\rVert_{H^{s}\left(\mathcal{M}\right)}^{2}=\left\langle
f,f\right\rangle_{H^{s}\left(\mathcal{M}\right)}<\infty$ where

$\left\langle
f,g\right\rangle_{H^{s}\left(\mathcal{M}\right)}=\sum_{j=0}^{\infty}(1+\lambda_{j})^{s}\left\langle
f,f_{j}\right\rangle_{L^{2}(\mathcal{M})}\left\langle
g,f_{j}\right\rangle_{L^{2}(\mathcal{M})}.$ (22)

######  Remark 10.

It is easy to see that substituting $(1+\lambda_{j})^{s}$ in (22) with
$\beta(\alpha+\lambda_{j})^{s}$ or with $(\alpha+\beta\lambda_{j}^{s})$ for
any $\alpha,\beta>0$ results in the same set of functions and an equivalent
norm. The former follows from [11, eq. (109)]. The latter follows from the
Binomial Theorem.

For Besov spaces we follow [14, 12] and define them in terms of approximations
by low-frequency functions. We fix a function
$\Phi\in\mathcal{C}^{\infty}\left(\mathbb{R}_{\geq 0},\mathbb{R}_{\geq
0}\right)$ such that $K=\operatorname{supp}\left(\Phi\right)\subset[0,2]$ and
$\Phi(x)=1$ for $x\in[0,1]$. We also define the functions
$\Phi_{j}(x)=\Phi\left(2^{-j}x\right)$.

[14, Corollary 3.6] shows that the operators $\Phi_{j}(\sqrt{\Delta})$ defined
by

$\Phi_{j}\left(\sqrt{\Delta}\right)f=\sum_{j\geq
0}\Phi_{j}\left(\sqrt{\lambda_{j}}\right)\left\langle f_{j},f\right\rangle
f_{j}$ (23)

are bounded in the space $L^{p}\left(\mathcal{M}\right)$ for all $1\leq
p\leq\infty$.444The space $L^{p}\left(\mathcal{M}\right)$ is the Banach space
of functions (or rather their equivalence classes) that are integrable when
raised to the power $p<\infty$ or essentially bounded for $p=\infty$. See for
instance [50] for details. Moreover, the same result also shows that we can
express any $f\in L^{p}\left(\mathcal{M}\right)$ as
$f=\lim_{j\to\infty}\Phi_{j}(\sqrt{\Delta})f$ in
$L^{p}\left(\mathcal{M}\right)$. $\Phi_{j}(\sqrt{\Delta})f$ can intuitively be
considered as a version of $f$ filtered by a low-pass filter. The next
definition introduces the Besov spaces $B_{p,q}^{s}\left(\mathcal{M}\right)$,
which are formulated in terms of quality-of-approximation by low-frequency
functions.

######  Definition 11 (Besov spaces).

For any $s>0$ and $1\leq p,q\leq\infty$ we define the Besov space
$B_{p,q}^{s}\left(\mathcal{M}\right)$ on the manifold $\mathcal{M}$ as the
space of functions $f\in L^{p}\left(\mathcal{M}\right)$ such that $\left\lVert
f\right\rVert_{B_{p,q}^{s}\left(\mathcal{M}\right)}<\infty$ where

$\displaystyle\left\lVert
f\right\rVert_{B_{p,q}^{s}\left(\mathcal{M}\right)}=\begin{cases}\left\lVert
f\right\rVert_{L^{p}(\mathcal{M})}+\left(\sum_{j\geq
0}\left(2^{js}\lVert\Phi_{j}(\sqrt{\Delta})f-f\rVert_{L^{p}(\mathcal{M})}\right)^{q}\right)^{1/q}&\text{if
}q<+\infty\\\ \left\lVert f\right\rVert_{L^{p}(\mathcal{M})}+\sup_{j\geq
0}2^{js}\lVert\Phi_{j}(\sqrt{\Delta})f-f\rVert_{L^{p}(\mathcal{M})}&\text{if
}q=+\infty.\end{cases}$ (24)

The classical Besov spaces $B_{2,2}^{s}$ coincide with the Sobolev spaces
$H^{s}$ on $\mathbb{R}^{d}$, in the sense that they define the same set of
functions and equivalent norms—see for instance [21] section 4.3.6—and even on
manifolds if one follows the construction of [51, 7.3-7.4] for Besov spaces.
Since our definition of the Besov spaces is somewhat non-standard, we present
a proof.

######  Proposition 12.

For all $s>0$,
$H^{s}\left(\mathcal{M}\right)=B_{2,2}^{s}\left(\mathcal{M}\right)$ as sets
and there exist two constants $C_{1},C_{2}>0$ such that for all $f\in
H^{s}\left(\mathcal{M}\right)=B_{2,2}^{s}\left(\mathcal{M}\right)$ we have

$C_{1}\left\lVert f\right\rVert_{H^{s}\left(\mathcal{M}\right)}\leq\left\lVert
f\right\rVert_{B_{2,2}^{s}\left(\mathcal{M}\right)}\leq C_{2}\left\lVert
f\right\rVert_{H^{s}\left(\mathcal{M}\right)}.$ (25)

###### Proof.

It is enough to prove (25), the rest will follow automatically. The main
technical tools used in the proof are Result 7 and summation by parts. First,
we prove the upper bound. Denote the support set of $\Phi$ by
$K=\operatorname{supp}(\Phi)$. Notice that

$\displaystyle\big{(}\left\lVert
f\right\rVert_{B_{2,2}^{s}(\mathcal{M})}-\left\lVert
f\right\rVert_{L^{2}(\mathcal{M})}\big{)}^{2}=$ $\displaystyle\sum_{j\geq
0}2^{2js}\big{\lVert}\Phi_{j}\left(\sqrt{\Delta}\right)f-f\big{\rVert}_{L^{2}(\mathcal{M})}^{2}$
(26) $\displaystyle=$ $\displaystyle\sum_{j\geq
0}2^{2js}\sum_{l:\sqrt{\lambda_{l}}\notin 2^{j}K}\big{|}\left\langle
f_{l},f\right\rangle_{L^{2}(\mathcal{M})}\big{|}^{2}$ (27) $\displaystyle\leq$
$\displaystyle\sum_{j\geq
0}2^{2js}\sum_{l:\sqrt{\lambda_{l}}>2^{j}}\big{|}\left\langle
f_{l},f\right\rangle_{L^{2}(\mathcal{M})}\big{|}^{2}.$ (28)

The last inequality results from the fact that $[0,1]\subset K$. By Weyl’s Law
(Result 7), there exists a constant $C>0$ such that $\lambda_{l}\leq
Cl^{2/d}$. Without loss of generality, we can assume that $C=2^{2r}$,
$r\in\mathbb{Z}_{>0}$. Hence $\sqrt{\lambda_{l}}>2^{j}$ implies
$l>2^{d(j-r)}$, thus we have

$\displaystyle\big{(}\left\lVert
f\right\rVert_{B_{2,2}^{s}(\mathcal{M})}\\!-\\!\left\lVert
f\right\rVert_{L^{2}(\mathcal{M})}\big{)}^{2}\\!\\!\leq$
$\displaystyle\sum_{j\geq 0}2^{2js}\sum_{l>2^{d(j-r)}}\big{|}\left\langle
f_{l},f\right\rangle_{L^{2}(\mathcal{M})}\big{|}^{2}$ (29) $\displaystyle=$
$\displaystyle\sum_{0\leq j\leq
r}\\!\\!2^{2js}\\!\\!\\!\\!\\!\\!\sum_{l>2^{d(j-r)}}\\!\\!\\!\\!\\!\big{|}\left\langle
f_{l},f\right\rangle_{L^{2}(\mathcal{M})}\big{|}^{2}+\sum_{j>r}2^{2js}\\!\\!\\!\\!\\!\\!\sum_{l>2^{d(j-r)}}\\!\\!\\!\\!\big{|}\left\langle
f_{l},f\right\rangle_{L^{2}(\mathcal{M})}\big{|}^{2}$ (30) $\displaystyle\leq$
$\displaystyle r2^{2rs}\left\lVert
f\right\rVert_{L^{2}(\mathcal{M})}^{2}+2^{2rs}\sum_{j>0}2^{2js}\sum_{l>2^{dj}}\big{|}\left\langle
f_{l},f\right\rangle_{L^{2}(\mathcal{M})}\big{|}^{2}.$ (31)

Now let $R_{j}=\sum_{l>2^{dj}}\big{|}\left\langle
f_{l},f\right\rangle_{L^{2}(\mathcal{M})}\big{|}^{2}$ and
$S_{J}=\sum_{j=1}^{J}2^{2js}\leq\frac{2^{2s}}{2^{2s}-1}2^{2Js},S_{0}=0$. Write

$\displaystyle\sum_{j>0}2^{2js}\sum_{l>2^{dj}}\big{|}\left\langle
f_{l},f\right\rangle_{L^{2}(\mathcal{M})}\big{|}^{2}=$
$\displaystyle\sum_{j>0}\left(S_{j}-S_{j-1}\right)R_{j}$ (32) $\displaystyle=$
$\displaystyle\sum_{j>0}S_{j}\left(R_{j}-R_{j+1}\right)$ (33) $\displaystyle=$
$\displaystyle\sum_{j>0}S_{j}\sum_{2^{dj}<l\leq 2^{d(j+1)}}\big{|}\left\langle
f_{l},f\right\rangle_{L^{2}(\mathcal{M})}\big{|}^{2}$ (34) $\displaystyle\leq$
$\displaystyle\frac{2^{2s}}{2^{2s}-1}\sum_{j>0}2^{2js}\sum_{2^{dj}<l\leq
2^{d(j+1)}}\big{|}\left\langle
f_{l},f\right\rangle_{L^{2}(\mathcal{M})}\big{|}^{2}$ (35) $\displaystyle\leq$
$\displaystyle\frac{2^{2s}}{2^{2s}-1}\sum_{j>0}\,\,\,\,\sum_{2^{dj}<l\leq
2^{d(j+1)}}l^{2s/d}\big{|}\left\langle
f_{l},f\right\rangle_{L^{2}(\mathcal{M})}\big{|}^{2}$ (36) $\displaystyle\leq$
$\displaystyle\frac{c^{s}2^{2s}}{2^{2s}-1}\sum_{j>0}\,\,\,\,\sum_{2^{dj}<l\leq
2^{d(j+1)}}\lambda_{l}^{s}\big{|}\left\langle
f_{l},f\right\rangle_{L^{2}(\mathcal{M})}\big{|}^{2}$ (37) $\displaystyle=$
$\displaystyle\frac{c^{s}2^{2s}}{2^{2s}-1}\sum_{l>2^{d}}\lambda_{l}^{s}\big{|}\left\langle
f_{l},f\right\rangle_{L^{2}(\mathcal{M})}\big{|}^{2}$ (38) $\displaystyle\leq$
$\displaystyle\frac{c^{s}2^{2s}}{2^{2s}-1}\sum_{l\geq
0}\lambda_{l}^{s}\big{|}\left\langle
f_{l},f\right\rangle_{L^{2}(\mathcal{M})}\big{|}^{2}$ (39)

where we have used Result 7 to get existence of a $c$ such that $l^{2/d}\leq
c\lambda_{l}$. This proves the upper bound. The proof for the lower bound is
similar. ∎

Proposition 12 provides a characterization of the Sobolev spaces
$H^{s}(\mathcal{M})$. There is yet another characterization of these spaces
that will be useful later, in terms of charts. We present this
characterization as part of the following result.

######  Theorem 13.

On the Sobolev space $H^{s}\left(\mathcal{M}\right)$, the following norms are
equivalent:

$\displaystyle\left\lVert f\right\rVert_{\,H^{s}\,(\mathcal{M})}$
$\displaystyle=\Bigg{(}\sum_{j=0}^{\infty}(1+\lambda_{j})^{s}\left\langle
f,f_{j}\right\rangle_{L^{2}(\mathcal{M})}^{2}\Bigg{)}^{1/2}$ (40)
$\displaystyle\left\lVert f\right\rVert_{B_{2,2}^{s}(\mathcal{M})}$
$\displaystyle=\left\lVert
f\right\rVert_{L^{2}(\mathcal{M})}+\Bigg{(}\sum_{j\geq
0}\left(2^{js}\lVert\Phi_{j}(\sqrt{\Delta})f-f\rVert_{L^{2}(\mathcal{M})}\right)^{2}\Bigg{)}^{1/2}$
(41) $\displaystyle\left\lVert
f\right\rVert_{H_{\mathcal{T}}^{s}(\mathcal{M})}$
$\displaystyle=\left(\sum_{l=1}^{L}\left\lVert\left(\chi_{l}f\right)\circ\phi_{l}^{-1}\right\rVert_{H^{s}\left(\mathbb{R}^{d}\right)}^{2}\right)^{1/2}$
(42)

###### Proof.

The equivalence between $\left\lVert\cdot\right\rVert_{H^{s}(\mathcal{M})}$
and $\left\lVert f\right\rVert_{B_{2,2}^{s}(\mathcal{M})}$ is given by
Proposition 12. The proof of equivalence between
$\left\lVert\cdot\right\rVert_{H^{s}(\mathcal{M})}$ and $\left\lVert
f\right\rVert_{H_{\mathcal{T}}^{s}(\mathcal{M})}$ can be found in [16]. ∎

###  A.4 Gaussian Random Elements

Here we recall the definition of a Gaussian process as a Banach-space-valued
random variable, following for instance [59].

######  Definition 14 (Gaussian random element).

Let $\left(\mathbb{B},\left\lVert\cdot\right\rVert_{\mathbb{B}}\right)$ be a
Banach space, and $f$ be a Borel random variable with values in $\mathbb{B}$
almost surely. We say that $f$ is a Gaussian random element if
$b^{*}\left(f\right)$ is a univariate Gaussian random variable for every
bounded linear functional $b^{*}:\mathbb{B}\to\mathbb{R}$.

Random variables of this kind are also sometimes called _Gaussian in the sense
of duality_. One should think of a Gaussian random element as a generalization
of a Gaussian process, but which is better-behaved from a function-analytic
point of view and in particular does not require the process to be an actual
function—as opposed to, for instance, a measure or a distribution. Many
connections between the usual Gaussian processes and Gaussian random elements
exist, see [32], [20, Appendix I], [57] for details. The following observation
about Gaussian random elements will be useful later.

######  Lemma 15.

A Gaussian process $f$ on the manifold $\mathcal{M}$ with almost surely
continuous sample paths is a Gaussian random element in the Banach space
$\left(\mathcal{C}\left(\mathcal{M}\right),\left\lVert\cdot\right\rVert_{\mathcal{C}(\mathcal{M})}\right)$
of continuous functions on $\mathcal{M}$.

###### Proof.

Since $\mathcal{C}\left(\mathcal{M}\right)$ is separable, this follows from
Lemma I.6 in [20]. ∎

###  A.5 A Technical Lemma

In order to apply Bernstein’s inequality when going from the error at input
locations to the $L^{2}(p_{0})$ error, we will use the following technical
extrapolation lemma.

######  Lemma 16.

For any function $g:\mathcal{M}\to\mathbb{R}$, a number
$\gamma\in\mathbb{R}_{>0}\setminus\mathbb{Z}_{>0}$ and a density
$p_{0}:\mathcal{M}\to\mathbb{R}_{>0}$ with $1\lesssim p_{0}$, we have

$\left\lVert g\right\rVert_{L^{\infty}(\mathcal{M})}\lesssim\left\lVert
g\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}^{\frac{d}{2\gamma+d}}\left\lVert
g\right\rVert_{L^{2}(p_{0})}^{\frac{2\gamma}{2\gamma+d}}.$ (43)

###### Proof.

We use Lemma 15 from [56] and push it through charts. More precisely we have,
using
$B_{\infty,\infty}^{\gamma}\left([0,1]^{D}\right)=\mathcal{CH}^{\gamma}\left([0,1]^{D}\right)$
for $\gamma\notin\mathbb{Z}_{>0}$, that

$\displaystyle\left\lVert g\right\rVert_{L^{\infty}(\mathcal{M})}\leq$
$\displaystyle\sum_{l}\left\lVert\left(\chi_{l}g\right)\circ\phi_{l}^{-1}\right\rVert_{L^{\infty}\left(\mathcal{V}_{l}\right)}$
(44) $\displaystyle\lesssim$
$\displaystyle\max_{l}\left\lVert\left(\chi_{l}g\right)\circ\phi_{l}^{-1}\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{V}_{l}\right)}^{\frac{d}{2\gamma+d}}\left\lVert\left(\chi_{l}g\right)\circ\phi_{l}^{-1}\right\rVert_{L^{2}\left(\mathcal{V}_{l}\right)}^{\frac{2\gamma}{2\gamma+d}}.$
(45)

By definition of the the manifold Hölder spaces, this gives

$\left\lVert g\right\rVert_{L^{\infty}(\mathcal{M})}\lesssim\left\lVert
g\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}^{\frac{d}{2\gamma+d}}\max_{l}\left\lVert\left(\chi_{l}g\right)\circ\phi_{l}^{-1}\right\rVert_{L^{2}\left(\mathcal{V}_{l}\right)}^{\frac{2\gamma}{2\gamma+d}}.$
(46)

Finally, since $\chi_{l}$ are bounded, the charts are smooth, and $p_{0}$ is
lower bounded, we have

$\left\lVert\left(\chi_{l}g\right)\circ\phi_{l}^{-1}\right\rVert_{L^{2}\left(\mathcal{V}_{l}\right)}^{2}=\int_{\mathcal{V}_{l}}\left|\left(\chi_{l}g\right)\circ\phi_{l}^{-1}(y)\right|^{2}dy\lesssim\int_{\mathcal{U}_{l}}g^{2}(x)p_{0}(x)\mu(dx)\lesssim\left\lVert
g\right\rVert_{L^{2}(\mathcal{M})}^{2}$ (47)

which gives the result. ∎

This concludes the necessary preliminaries. We now turn to the proofs of our
main results.

##  Appendix B Proofs

Throughout this section we use the notation from Appendix A and results from
Appendices C, D and E. We start by defining our notation for Gaussian
likelihood and probability distribution of the sample.

######  Definition 17.

For every $\boldsymbol{x}\in\mathcal{M}^{n}$ and $f:\mathcal{M}\to\mathbb{R}$
we define $p^{(f)}_{\boldsymbol{x},\boldsymbol{y}}$ to be the joint
distribution corresponding to the marginal $p_{\boldsymbol{x}}=p_{0}$ and
conditional
$p^{(f)}_{\boldsymbol{y}\mid\boldsymbol{x}}=\operatorname{N}(f(\boldsymbol{x}),\sigma_{\varepsilon}^{2}\mathbf{I})$,
where $f(\boldsymbol{x})$ is the vector with entries $f(x_{i})$. Expectations
with respect to $p^{(f)}_{\boldsymbol{x},\boldsymbol{y}}$, to
$p_{\boldsymbol{x}}$ and to $p^{(f)}_{\boldsymbol{y}\mid\boldsymbol{x}}$ we
denote by $\operatorname{\mathbb{E}}^{(f)}_{\boldsymbol{x},\boldsymbol{y}}$,
by $\operatorname{\mathbb{E}}_{\boldsymbol{x}}$ and by
$\operatorname{\mathbb{E}}^{(f)}_{\boldsymbol{y}\mid\boldsymbol{x}}$,
respectively. Sometimes, when $f$ is clear from the context, we omit the
subscript (f).

Using [56, Theorem 1], which is valid for any compact metric space and thus
also for $\mathcal{M}$, we can deduce a posterior contraction rate at data
input locactions, with respect to the _empirical $L^{2}$-norm_555This is
actually a seminorm, but we follow the rest of the literature in referring to
it as a norm.

$\left\lVert
f\right\rVert_{n}=\left(\frac{1}{n}\sum_{i=1}^{n}f(x_{i})^{2}\right)^{1/2}$
(48)

by studying the _concentration functions_ with respect to the uniform norm.
This is the purpose of the following lemma, whose proof mainly follows [56],
but crucially relies on two new components: (i) the _small ball asymptotics_
for manifolds we study in Appendix E and (ii) the characterization of the RKHS
corresponding to both intrinsic and extrinsic priors as manifold Sobolev
spaces we obtain in Appendix C. We recall that the prior $\Pi_{n}$ may depend
on $n$ if we consider a truncated intrinsic Matérn process.

######  Theorem 18.

Let $\Pi_{n}$ denote the prior in either Definition 4, Section 3.2 or
Definition 5 with smoothness parameter $\nu>d/2$. Let $\mathbb{H}_{n}$ denote
the corresponding RKHS. Define the _concentration function_ for $f_{0}\in
C(\mathcal{M})$ and $\varepsilon>0$ by

$\varphi_{f_{0}}\left(\varepsilon\right)=-\ln\operatorname{\mathbb{P}}_{f\sim\Pi_{n}}\left[\left\lVert
f\right\rVert_{L^{\infty}(\mathcal{M})}<\varepsilon\right]+\inf_{f\in\mathbb{H}_{n}:\,\left\lVert
f-f_{0}\right\rVert_{L^{\infty}(\mathcal{M})}<\varepsilon}\left\lVert
f\right\rVert_{\mathbb{H}_{n}}^{2}.$ (49)

Then if $f_{0}\in H^{\beta}\left(\mathcal{M}\right)\cap
B_{\infty,\infty}^{\beta}\left(\mathcal{M}\right),\beta>0$ we have
$\varphi_{f_{0}}\left(\varepsilon_{n}\right)\leq n\varepsilon_{n}^{2}$ for
$\varepsilon_{n}$ a multiple of $n^{-\frac{\min(\nu,\beta)}{2\nu+d}}$.

###### Proof.

The first term on the right-hand side of Equation 49 is bounded by
$C\varepsilon^{-d/\nu}$ by Lemma 30. To bound the second term, we assume,
without loss of generality,666If $\beta>\nu$ then $f_{0}\in
H^{\beta}\left(\mathcal{M}\right)\cap B_{\infty,\infty}^{\beta}\subseteq
H^{\nu}\left(\mathcal{M}\right)\cap
B_{\infty,\infty}^{\nu}\left(\mathcal{M}\right)$ gives $\varepsilon_{n}\propto
n^{-\frac{\nu}{2\nu+d}}=n^{-\frac{\min(\beta,\nu)}{2\nu+d}}$. that
$\nu\geq\beta$. Consider an approximation $f=\Phi_{j}(\sqrt{\Delta})f_{0}$ of
$f_{0}$, where $c\varepsilon\leq 2^{-\beta j}\leq\varepsilon$ and $c>0$ does
not depend on $j$. Since we assume $f_{0}\in
B_{\infty,\infty}^{\beta}(\mathcal{M})$, by definition of
$B_{\infty,\infty}^{\beta}(\mathcal{M})$ we have

$\left\lVert f_{0}-f\right\rVert_{L^{\infty}(\mathcal{M})}\leq\left\lVert
f_{0}\right\rVert_{B_{\infty,\infty}^{\beta}(\mathcal{M})}2^{-\beta
j}\lesssim\varepsilon$ (50)

where in the last inequality the $B_{\infty,\infty}^{\beta}(\mathcal{M})$-norm
is the constant implied by notation $\lesssim$. We now show that

$\left\lVert
f\right\rVert_{\mathbb{H}_{n}}^{2}\lesssim\varepsilon^{-\frac{2}{\beta}\left(\nu-\beta+d/2\right)}.$
(51)

First notice that by Lemmas 21 and 24, for any prior considered here we have
$\mathbb{H}_{n}\subseteq H^{\nu+d/2}\left(\mathcal{M}\right)$ and
$\left\lVert\cdot\right\rVert_{\mathbb{H}_{n}}\lesssim\left\lVert\cdot\right\rVert_{H^{\nu+d/2}\left(\mathcal{M}\right)}$
with a constant that does not depend on $n$. Hence using Result 7 (Weyl’s Law)
and properties of $\Phi$, we have

$\displaystyle\left\lVert f\right\rVert_{\mathbb{H}_{n}}^{2}$
$\displaystyle\lesssim\left\lVert
f\right\rVert_{H^{\nu+d/2}\left(\mathcal{M}\right)}^{2}$ (52)
$\displaystyle=\sum_{l\geq
0}(1+\lambda_{l})^{\nu+d/2}\Phi^{2}\left(2^{-j}\sqrt{\lambda_{l}}\right)\left|\left\langle
f_{l},f_{0}\right\rangle_{L^{2}(\mathcal{M})}\right|^{2}$ (53)
$\displaystyle\leq\sum_{l:\sqrt{\lambda_{l}}\leq
2^{j+1}}(1+\lambda_{l})^{\nu+d/2-\beta}(1+\lambda_{l})^{\beta}\left|\left\langle
f_{l},f_{0}\right\rangle_{L^{2}(\mathcal{M})}\right|^{2}$ (54)
$\displaystyle\leq 2^{(j+2)(\nu+d/2-\beta)}\sum_{l:\sqrt{\lambda_{l}}\leq
2^{j+1}}(1+\lambda_{l})^{\beta}\left|\left\langle
f_{l},f_{0}\right\rangle_{L^{2}(\mathcal{M})}\right|^{2}$ (55)
$\displaystyle\leq 2^{(j+2)(\nu+d/2-\beta)}\sum_{l\geq
0}(1+\lambda_{l})^{\beta}\left|\left\langle
f_{l},f_{0}\right\rangle_{L^{2}(\mathcal{M})}\right|^{2}$ (56)
$\displaystyle=2^{(j+2)(\nu+d/2-\beta)}\left\lVert
f_{0}\right\rVert_{H^{\beta}\left(\mathcal{M}\right)}^{2}$ (57)
$\displaystyle\lesssim\varepsilon^{-\frac{1}{\beta}(\nu+d/2-\beta)}\left\lVert
f_{0}\right\rVert_{H^{\beta}\left(\mathcal{M}\right)}^{2}\lesssim\varepsilon^{-\frac{2}{\beta}(\nu+d/2-\beta)}\left\lVert
f_{0}\right\rVert_{H^{\beta}\left(\mathcal{M}\right)}^{2}.$ (58)

Our assumption $\nu\geq\beta$ implies that

$\frac{2}{\beta}\left(\nu-\beta+d/2\right)\geq\frac{d}{\beta}\geq\frac{d}{\nu}.$
(59)

Hence, we have
$\varepsilon^{-d/\nu}\leq\varepsilon^{-\frac{2}{\beta}(\nu-\beta+d/2)}$ which
gives us
$\varphi_{f_{0}}\left(\varepsilon\right)\lesssim\varepsilon^{-\frac{2}{\beta}(\nu-\beta+d/2)}$.
It is then easy to check that $\varepsilon_{n}=Mn^{-\frac{\beta}{2\nu+d}}$
satisfies $\varphi_{f_{0}}\left(\varepsilon_{n}\right)\leq
n\varepsilon_{n}^{2}$ for $M>0$ large enough. ∎

From this, we deduce an upper bound on the error in _the empirical $L^{2}$
norm_ $\left\lVert\cdot\right\rVert_{n}$, that is, on the Euclidean distance
between the posterior Gaussian process $f$ and the ground truth function
$f_{0}$ evaluated at data locations $x_{i}$.

######  Lemma 19.

Let $\Pi_{n}$ denote the prior in either Definition 4, Section 3.2 or
Definition 5 with smoothness parameter $\nu>d/2$. Fix $f_{0}\in
H^{\beta}\left(\mathcal{M}\right)\cap
B_{\infty,\infty}^{\beta}\left(\mathcal{M}\right)$ with $\beta>0$. Then

$\operatorname{\mathbb{E}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left\lVert
f-f_{0}\right\rVert_{n}^{2}\leq\varepsilon_{n}^{2}$ (60)

for $\varepsilon_{n}\propto n^{-\frac{\min(\nu,\beta)}{2\nu+d}}$ with constant
depending on $f_{0},\nu$ but not on $\boldsymbol{x}$.

###### Proof.

By Theorem 18 for $\varepsilon_{n}$ a multiple of
$n^{-\frac{\min(\beta,\nu)}{2\nu+d}}$, we have
$\varphi_{f_{0}}\left(\varepsilon_{n}\right)\leq n\varepsilon_{n}^{2}$. By
virtue of this, the proof of Theorem 1 and Proposition 11 of [56] imply the
result. Indeed, the proof of Theorem 1 relies solely on the fact that
$\varphi_{f_{0}}\left(\varepsilon_{n}/2\right)\leq n\varepsilon_{n}^{2}$ and
an application of [56, Proposition 11]. We have
$\varphi_{f_{0}}\left(\varepsilon_{n}\right)\leq n\varepsilon_{n}^{2}\leq
n\left(2\varepsilon_{n}\right)^{2}$ and hence the condition is satisfied with
$\varepsilon_{n}$ replaced by $2\varepsilon_{n}$. ∎

We now turn to the proofs of our main results, Definitions 4, 3.2 and 5, which
for convenience we restate below. The idea of these proofs is to extend the
result of Lemma 19 from input locations to the whole manifold $\mathcal{M}$
using an appropriate concentration inequality. To this end, the proof closely
follows the one of Theorem 2 in [56], but relies on Lemma 19 proved above and
on the concentration inequality we prove in Appendix D along with some
important sample differentiablity properties of the prior processes.

*

*

*

###### Proof of Definitions 4, 3.2 and 5.

To ease notation, we denote the expectations under the true data generating
process by
$\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}=\operatorname{\mathbb{E}}^{(f_{0})}_{\boldsymbol{x},\boldsymbol{y}}$
and
$\operatorname{\mathbb{E}}_{\boldsymbol{y}\mid\boldsymbol{x}}=\operatorname{\mathbb{E}}^{(f_{0})}_{\boldsymbol{y}\mid\boldsymbol{x}}$,
omitting the superscript ${(\cdot)}^{f_{0}}$. Take $\varepsilon_{n}\propto
n^{-\frac{\min(\beta,\nu)}{2\nu+d}}$ satisfying
$\varphi_{f_{0}}\left(\varepsilon_{n}/2\right)\leq n\varepsilon_{n}^{2}$,
noting that such a rate exists by Theorem 18. Then, for each $n$, there exists
an element $f_{n}\in\mathbb{H}_{n}$, where $\mathbb{H}_{n}$ is the RKHS
corresponding to $\Pi_{n}$, satisfying

$\displaystyle\left\lVert f_{n}\right\rVert_{\mathbb{H}_{n}}^{2}\leq
n\varepsilon_{n}^{2}$ $\displaystyle\left\lVert
f_{n}-f_{0}\right\rVert_{L^{\infty}(\mathcal{M})}\leq\varepsilon_{n}/2.$ (61)

Write

$\displaystyle\varepsilon_{n}^{-2}\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\operatorname{\mathbb{E}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left\lVert
f-f_{0}\right\rVert_{L^{2}\left(p_{0}\right)}^{2}$
$\displaystyle\lesssim\varepsilon_{n}^{-2}\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\operatorname{\mathbb{E}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left\lVert
f_{n}-f_{0}\right\rVert_{L^{2}\left(p_{0}\right)}^{2}$ (62)
$\displaystyle\quad+\varepsilon_{n}^{-2}\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\operatorname{\mathbb{E}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left\lVert
f-f_{n}\right\rVert_{L^{2}\left(p_{0}\right)}^{2}$ (63) $\displaystyle\lesssim
1+\varepsilon_{n}^{-2}\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\operatorname{\mathbb{E}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left\lVert
f-f_{n}\right\rVert_{L^{2}\left(p_{0}\right)}^{2}.$ (64)

Thus, we can work with $f-f_{n}$ instead of $f-f_{0}$. Define
$\mathcal{B}(r)=\big{\\{}\left\lVert
f-f_{n}\right\rVert_{L^{2}\left(p_{0}\right)}>\varepsilon_{n}r\big{\\}}$. Then

$\varepsilon_{n}^{-2}\operatorname{\mathbb{E}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left\lVert
f-f_{n}\right\rVert_{L^{2}\left(p_{0}\right)}^{2}.=\int_{0}^{\infty}2r\operatorname{\mathbb{P}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left(\mathcal{B}(r)\right)\mathop{}\\!\mathrm{d}r.$
(65)

Fix a $\gamma$ such that $d/2<\gamma<\nu,\gamma\notin\mathbb{Z}_{>0}$ and
$s>0,\tau>0$ and define

$\displaystyle\mathcal{B}^{(\operatorname{I})}(r)$
$\displaystyle=\left\\{2\left\lVert
f-f_{n}\right\rVert_{n}>\varepsilon_{n}r\right\\}$ (66)
$\displaystyle\mathcal{B}^{(\operatorname{II})}(r)$
$\displaystyle=\left\\{\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}>\tau\sqrt{n}\varepsilon_{n}r^{s}\right\\}$
(67) $\displaystyle\mathcal{B}^{(\operatorname{III})}(r)$
$\displaystyle=\mathcal{B}(r)\setminus\left(\mathcal{B}^{(\operatorname{I})}(r)\cup\mathcal{B}^{(\operatorname{II})}(r)\right).$
(68)

Then
$\mathcal{B}(r)\subseteq\mathcal{B}^{(\operatorname{I})}(r)\cup\mathcal{B}^{(\operatorname{II})}(r)\cup\mathcal{B}^{(\operatorname{III})}(r)$,
and thus for an indexed family of events $\mathcal{A}_{r}$ to be chosen later,
we have

$\displaystyle\varepsilon_{n}^{-2}\operatorname{\mathbb{E}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left\lVert
f-f_{n}\right\rVert_{L^{2}\left(p_{0}\right)}^{2}\lesssim\int_{0}^{\infty}r\operatorname{\mathbb{P}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left(\mathcal{B}^{(\operatorname{I})}(r)\right)\mathop{}\\!\mathrm{d}r+\int_{0}^{\infty}r\mathbbold{1}_{\mathcal{A}_{r}^{c}}$
$\displaystyle\mathop{}\\!\mathrm{d}r$ (69)
$\displaystyle+\int_{0}^{\infty}r\mathbbold{1}_{\mathcal{A}_{r}}\operatorname{\mathbb{P}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left(\mathcal{B}^{(\operatorname{II})}(r)\right)\mathop{}\\!\mathrm{d}r+\int_{0}^{\infty}r\mathbbold{1}_{\mathcal{A}_{r}}\operatorname{\mathbb{P}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left(\mathcal{B}^{(\operatorname{III})}(r)\right)$
$\displaystyle\mathop{}\\!\mathrm{d}r.$ (70)

For the first term, by Lemma 19 applied conditionally on $\boldsymbol{x}$, for
which we got a bound on the integrated empirical $L^{2}$-norm uniformly on the
design points, we have

$\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\int_{0}^{\infty}r\operatorname{\mathbb{P}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left(\mathcal{B}^{(\operatorname{I})}(r)\right)\mathop{}\\!\mathrm{d}r\lesssim\varepsilon_{n}^{-2}\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\operatorname{\mathbb{E}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left\lVert
f-f_{0}\right\rVert_{n}^{2}\lesssim\varepsilon_{n}^{-2}\varepsilon_{n}^{2}=1.$
(71)

Moreover, by Lemma 14 of [56] applied with $r$ in the notation of the
reference being equal to $\sqrt{n}\varepsilon_{n}r^{s}$, for each $r>0$, the
event

$\displaystyle\mathcal{A}_{r}(\boldsymbol{x})$
$\displaystyle=\left\\{\\!\boldsymbol{u}\in\mathbb{R}^{n}\\!:\\!\\!\\!\int\\!\frac{p_{\boldsymbol{y}\mid\boldsymbol{x}}^{(f)}(\boldsymbol{u})}{p_{\boldsymbol{y}\mid\boldsymbol{x}}^{(f_{0})}(\boldsymbol{u})}\mathop{}\\!\mathrm{d}\Pi_{n}\left(f\right)\geq
e^{-n\varepsilon_{n}^{2}r^{2s}}\operatorname{\mathbb{P}}_{f\sim\Pi_{n}}\left(\left\lVert
f\\!-\\!f_{0}\right\rVert_{L^{\infty}(\mathcal{M})}\\!<\\!\varepsilon_{n}r^{s}\right)\\!\\!\right\\}$
(72)
$\displaystyle\supseteq\left\\{\\!\boldsymbol{u}\in\mathbb{R}^{n}\\!:\\!\\!\\!\int\\!\frac{p_{\boldsymbol{y}\mid\boldsymbol{x}}^{(f)}(\boldsymbol{u})}{p_{\boldsymbol{y}\mid\boldsymbol{x}}^{(f_{0})}(\boldsymbol{u})}\mathop{}\\!\mathrm{d}\Pi_{n}\left(f\right)\geq
e^{-n\varepsilon_{n}^{2}r^{2s}}\operatorname{\mathbb{P}}_{f\sim\Pi_{n}}\left(\left\lVert
f\\!-\\!f_{0}\right\rVert_{n}\\!<\\!\varepsilon_{n}r^{s}\right)\\!\\!\right\\}$
(73)

is such that

$p_{\boldsymbol{y}\mid\boldsymbol{x}}^{(f_{0})}\left[\mathcal{A}_{r}^{c}(\boldsymbol{x})\right]\leq
e^{-n\varepsilon_{n}^{2}r^{2s}/8}.$ (74)

It should be noted that the $\sqrt{n}$ factor disappears because of the
discrepancy between the empirical $L^{2}$ norm
$\left\lVert\cdot\right\rVert_{n}$ and the Euclidean norm used in [56]. By
Fubini’s Theorem and since $n\varepsilon_{n}^{2}\geq n^{\frac{d}{2\nu+d}}\geq
1$, the second term is bounded by

$\displaystyle\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\int_{0}^{\infty}r\mathbbold{1}_{\mathcal{A}_{r}^{c}(\boldsymbol{x})}\mathop{}\\!\mathrm{d}r$
$\displaystyle=\int_{0}^{\infty}r\operatorname{\mathbb{E}}_{\boldsymbol{x}}\left[\operatorname{\mathbb{E}}_{\boldsymbol{y}\mid\boldsymbol{x}}\left[\mathbbold{1}_{\mathcal{A}_{r}^{c}(\boldsymbol{x})}\right]\right]\mathop{}\\!\mathrm{d}r$
(75)
$\displaystyle\leq\int_{0}^{\infty}re^{-n\varepsilon_{n}^{2}r^{2s}/8}\mathop{}\\!\mathrm{d}r$
(76)
$\displaystyle\leq\int_{0}^{\infty}re^{-r^{2s}/8}\mathop{}\\!\mathrm{d}r=C_{s}<\infty.$
(77)

It remains to bound the last two terms. By Bayes’ Rule, we have the equality

$\displaystyle\operatorname{\mathbb{P}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left(\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}>\tau\sqrt{n}\varepsilon_{n}r^{s}\right)$
$\displaystyle=\frac{\int_{\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}>\tau\sqrt{n}\varepsilon_{n}r^{s}}p^{(f)}_{\boldsymbol{y}\mid\boldsymbol{x}}(\boldsymbol{y})\mathop{}\\!\mathrm{d}\Pi_{n}(f)}{\int
p^{(f)}_{\boldsymbol{y}\mid\boldsymbol{x}}(\boldsymbol{y})\mathop{}\\!\mathrm{d}\Pi_{n}(f)}$
(78) $\displaystyle=\frac{\int_{\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}>\tau\sqrt{n}\varepsilon_{n}r^{s}}\frac{p^{(f)}_{\boldsymbol{y}\mid\boldsymbol{x}}(\boldsymbol{y})}{p^{(f_{0})}_{\boldsymbol{y}\mid\boldsymbol{x}}(\boldsymbol{y})}\mathop{}\\!\mathrm{d}\Pi_{n}(f)}{\int\frac{p^{(f)}_{\boldsymbol{y}\mid\boldsymbol{x}}(\boldsymbol{y})}{p^{(f_{0})}_{\boldsymbol{y}\mid\boldsymbol{x}}(\boldsymbol{y})}\mathop{}\\!\mathrm{d}\Pi_{n}(f)}$
(79)

therefore for $\boldsymbol{y}\in\mathcal{A}_{r}(\boldsymbol{x})$, we have

$\displaystyle\operatorname{\mathbb{P}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left(\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}>\tau\sqrt{n}\varepsilon_{n}r^{s}\right)$
(80)
$\displaystyle\quad\leq\frac{e^{n\varepsilon_{n}^{2}r^{2s}}}{\operatorname{\mathbb{P}}_{f\sim\Pi_{n}}\left(\left\lVert
f-f_{0}\right\rVert_{L^{\infty}(\mathcal{M})}<\varepsilon_{n}r^{s}\right)}\int_{\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}>\tau\sqrt{n}\varepsilon_{n}r^{s}}\frac{p^{(f)}_{\boldsymbol{y}\mid\boldsymbol{x}}(\boldsymbol{y})}{p^{(f_{0})}_{\boldsymbol{y}\mid\boldsymbol{x}}(\boldsymbol{y})}\mathop{}\\!\mathrm{d}\Pi_{n}(f).$
(81)

Hence taking expectations, using Tonelli’s Theorem, and
$\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\frac{p^{(f)}_{\boldsymbol{y}\mid\boldsymbol{x}}(\boldsymbol{y})}{p^{(f_{0})}_{\boldsymbol{y}\mid\boldsymbol{x}}(\boldsymbol{y})}=1$
gives

$\displaystyle\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\left[\mathbbold{1}_{\mathcal{A}_{r}(\boldsymbol{x})}\operatorname{\mathbb{P}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left(\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}>\tau\sqrt{n}\varepsilon_{n}r^{s}\right)\right]$
(82)
$\displaystyle\quad\leq\frac{e^{n\varepsilon_{n}^{2}r^{2s}}}{\operatorname{\mathbb{P}}_{f\sim\Pi_{n}}\left(\left\lVert
f-f_{0}\right\rVert_{L^{\infty}(\mathcal{M})}<\varepsilon_{n}r^{s}\right)}\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\int_{\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}>\tau\sqrt{n}\varepsilon_{n}r^{s}}\frac{p^{(f)}_{\boldsymbol{y}\mid\boldsymbol{x}}(\boldsymbol{y})}{p^{(f_{0})}_{\boldsymbol{y}\mid\boldsymbol{x}}(\boldsymbol{y})}\mathop{}\\!\mathrm{d}\Pi_{n}(f)$
(83)
$\displaystyle\quad=\frac{e^{n\varepsilon_{n}^{2}r^{2s}}}{\operatorname{\mathbb{P}}_{f\sim\Pi_{n}}\left(\left\lVert
f-f_{0}\right\rVert_{L^{\infty}(\mathcal{M})}<\varepsilon_{n}r^{s}\right)}\int_{\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}>\tau\sqrt{n}\varepsilon_{n}r^{s}}\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\frac{p^{(f)}_{\boldsymbol{y}\mid\boldsymbol{x}}(\boldsymbol{y})}{p^{(f_{0})}_{\boldsymbol{y}\mid\boldsymbol{x}}(\boldsymbol{y})}\mathop{}\\!\mathrm{d}\Pi_{n}(f)$
(84)
$\displaystyle\quad=\frac{e^{n\varepsilon_{n}^{2}r^{2s}}}{\operatorname{\mathbb{P}}_{f\sim\Pi_{n}}\left(\left\lVert
f-f_{0}\right\rVert_{L^{\infty}(\mathcal{M})}<\varepsilon_{n}r^{s}\right)}\operatorname{\mathbb{P}}_{f\sim\Pi_{n}}\left(\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}>\tau\sqrt{n}\varepsilon_{n}r^{s}\right).$
(85)

Therefore, the third term can be bounded as

$\displaystyle\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\int_{0}^{\infty}r\mathbbold{1}_{\mathcal{A}_{r}}\operatorname{\mathbb{P}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left(\mathcal{B}^{(\operatorname{II})}(r)\right)\mathop{}\\!\mathrm{d}r$
(86)
$\displaystyle\quad\leq\int_{0}^{\infty}r\frac{e^{n\varepsilon_{n}^{2}r^{2s}}}{\operatorname{\mathbb{P}}_{f\sim\Pi_{n}}\left(\left\lVert
f-f_{0}\right\rVert_{L^{\infty}(\mathcal{M})}<\varepsilon_{n}r^{s}\right)}\operatorname{\mathbb{P}}_{f\sim\Pi_{n}}\left(\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}>\tau\sqrt{n}\varepsilon_{n}r^{s}\right)\mathop{}\\!\mathrm{d}r.$
(87)

Now, using Lemma 28, for a possibly small constant $c>0$ independent of $n$,
we have

$\operatorname{\mathbb{P}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left(\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}>\tau\sqrt{n}\varepsilon_{n}r^{s}\right)\leq
e^{-c\tau^{2}n\varepsilon_{n}^{2}r^{2s}}.$ (88)

Moreover, by using the bound on the concentration function in Theorem 18 and
[20, Proposition 11.19], we can assume that

$\operatorname{\mathbb{P}}_{f\sim\Pi_{n}}\left[\left\lVert
f-f_{0}\right\rVert_{L^{\infty}(\mathcal{M})}<\varepsilon_{n}r^{s}\right]\geq
e^{-c^{-1}n\varepsilon_{n}^{2}r^{2s}}.$ (89)

Therefore, the third term is bounded by

$\displaystyle\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\int_{0}^{\infty}r\mathbbold{1}_{\mathcal{A}_{r}(\boldsymbol{x})}\operatorname{\mathbb{P}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left(\mathcal{B}^{(\operatorname{II})}(r)\right)\mathop{}\\!\mathrm{d}r$
$\displaystyle\leq\int_{0}^{\infty}re^{-n(c\tau^{2}-1)\varepsilon_{n}^{2}r^{2s}}e^{c^{-1}n\varepsilon_{n}^{2}r^{2s}}\mathop{}\\!\mathrm{d}r$
(90)
$\displaystyle\leq\int_{0}^{\infty}re^{-r^{2s}}\mathop{}\\!\mathrm{d}r<\infty$
(91)

if $\tau^{2}c>2+c^{-1}$. It remains to bound the last term. We have, by the
same arguments as above, that

$\displaystyle\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\int_{0}^{\infty}r\mathbbold{1}_{\mathcal{A}_{r}(\boldsymbol{x})}\operatorname{\mathbb{P}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left(\mathcal{B}^{(\operatorname{III})}(r)\right)\mathop{}\\!\mathrm{d}r$
(92)
$\displaystyle=\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\int_{0}^{\infty}r\mathbbold{1}_{\mathcal{A}_{r}(\boldsymbol{x})}$
(93)
$\displaystyle\quad\times\operatorname{\mathbb{P}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left(\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}\leq\tau\sqrt{n}\varepsilon_{n}r^{s},2\left\lVert
f-f_{n}\right\rVert_{n}\leq\varepsilon_{n}r\leq\left\lVert
f-f_{n}\right\rVert_{L^{2}(p_{0})}\right)\mathop{}\\!\mathrm{d}r$ (94)
$\displaystyle\leq\int_{0}^{\infty}r\frac{e^{n\varepsilon_{n}^{2}r^{2s}}}{\operatorname{\mathbb{P}}_{f\sim\Pi_{n}}\left(\left\lVert
f-f_{0}\right\rVert_{L^{\infty}(\mathcal{M})}<\varepsilon_{n}r^{s}\right)}$
(95)
$\displaystyle\quad\times\operatorname{\mathbb{E}}_{\boldsymbol{x}}\operatorname{\mathbb{P}}_{f\sim\Pi_{n}}\left(\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}\leq\tau\sqrt{n}\varepsilon_{n}r^{s},2\left\lVert
f-f_{n}\right\rVert_{n}\leq\varepsilon_{n}r\leq\left\lVert
f-f_{n}\right\rVert_{L^{2}(p_{0})}\right)\mathop{}\\!\mathrm{d}r$ (96)
$\displaystyle\leq\int_{0}^{\infty}re^{(c+1)n\varepsilon_{n}^{2}r^{2s}}$ (97)
$\displaystyle\quad\times\int_{\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}\leq\tau\sqrt{n}\varepsilon_{n}r^{s},\varepsilon_{n}r\leq\left\lVert
f-f_{n}\right\rVert_{L^{2}(p_{0})}}\operatorname{\mathbb{E}}_{\boldsymbol{x}}\mathbbold{1}_{\left\lVert
f-f_{n}\right\rVert_{L^{2}(p_{0})}\geq 2\left\lVert
f-f_{n}\right\rVert_{n}}\mathop{}\\!\mathrm{d}\Pi_{n}(f)\mathop{}\\!\mathrm{d}r.$
(98)

As the squared empirical $L^{2}$-norm is a Monte Carlo approximation of the
true $L^{2}$-norm, the probability in the integrand can be controlled via a
concentration inequality. As in [56], we use Bernstein’s inequality [58, Lemma
2.2.9]. For a collection $Y_{1},\ldots Y_{n}$ of random variables such that
$\operatorname{\mathbb{E}}Y_{i}=0$ and $Y_{i}\in[-M,M]$ almost surely for some
constant $M>0$, this inequality asserts

$\operatorname{\mathbb{P}}(\left|Y_{1}+\ldots+Y_{n}\right|>x)\leq
2\exp\left(-\frac{1}{2}\frac{x^{2}}{v+Mx/3}\right)$ (99)

where $v\geq\operatorname{Var}(Y_{1}+\ldots+Y_{n})$. We put
$Y_{i}=\frac{1}{n}\left(f(x_{i})-f_{n}(x_{i})\right)^{2}-\frac{1}{n}\left\lVert
f-f_{n}\right\rVert_{L^{2}(p_{0})}^{2}$ where $x_{i}$ are IID with $x_{i}\sim
p_{0}$. It is easy to check that $\operatorname{\mathbb{E}}Y_{i}=0$ and
$Y_{i}\in[-M,M]$ for $M=\frac{1}{n}\left\lVert
f-f_{n}\right\rVert_{L^{\infty}(\mathcal{M})}$. Furthermore,
$v=\frac{1}{n}\left\lVert
f-f_{n}\right\rVert_{L^{\infty}(\mathcal{M})}^{2}\left\lVert
f-f_{n}\right\rVert_{L^{2}(\mathcal{M})}^{2}$ upper-bounds the variance of the
respective sum, since

$\displaystyle\operatorname{Var}(Y_{1}+\ldots+Y_{n})$
$\displaystyle=\frac{1}{n}\operatorname{Var}_{x\sim
p_{0}}(f(x)-f_{n}(x))^{2}\leq\frac{1}{n}\operatorname{\mathbb{E}}_{x\sim
p_{0}}(f(x)-f_{n}(x))^{4}$ (100) $\displaystyle\leq\frac{\left\lVert
f-f_{n}\right\rVert_{L^{\infty}(\mathcal{M})}^{2}}{n}\operatorname{\mathbb{E}}_{x\sim
p_{0}}(f(x)-f_{n}(x))^{2}=v.$ (101)

Using Bernstein’s inequality with $Y_{i}$, $M$, $v$ as above and with
$x=\frac{3}{4}\left\lVert f-f_{n}\right\rVert_{L^{2}\left(p_{0}\right)}^{2}$,
we have

$\displaystyle\operatorname{\mathbb{E}}_{\boldsymbol{x}}\mathbbold{1}_{\left\lVert
f-f_{n}\right\rVert_{L^{2}(p_{0})}\geq 2\left\lVert f-f_{n}\right\rVert_{n}}$
$\displaystyle=\operatorname{\mathbb{E}}_{\boldsymbol{x}}\mathbbold{1}_{\left\lVert
f-f_{n}\right\rVert_{L^{2}(p_{0})}^{2}\geq 4\left\lVert
f-f_{n}\right\rVert_{n}^{2}}$ (102)
$\displaystyle=\operatorname{\mathbb{E}}_{\boldsymbol{x}}\mathbbold{1}_{\left\lVert
f-f_{n}\right\rVert_{n}^{2}-\left\lVert
f-f_{n}\right\rVert_{L^{2}(p_{0})}^{2}\leq-\frac{3}{4}\left\lVert
f-f_{n}\right\rVert_{L^{2}(p_{0})}^{2}}$ (103)
$\displaystyle\lesssim\exp\left(-\frac{1}{2}\frac{x^{2}}{v+Mx/3}\right)$ (104)
$\displaystyle=\exp\left(-\frac{\frac{9}{32}\left\lVert
f-f_{n}\right\rVert_{L^{2}\left(p_{0}\right)}^{4}}{v+\frac{1}{n}\left\lVert
f-f_{n}\right\rVert_{L^{\infty}(\mathcal{M})}^{2}\frac{3}{4}\left\lVert
f-f_{n}\right\rVert_{L^{2}\left(p_{0}\right)}^{2}/3}\right)$ (105)
$\displaystyle=\exp\left(-\frac{9n}{32}\frac{4}{5}\frac{\left\lVert
f-f_{n}\right\rVert_{L^{2}\left(p_{0}\right)}^{2}}{\left\lVert
f-f_{n}\right\rVert_{L^{\infty}(\mathcal{M})}^{2}}\right)$ (106)
$\displaystyle=\exp\left(-\frac{9n}{40}\frac{\left\lVert
f-f_{n}\right\rVert_{L^{2}(p_{0})}^{2}}{\left\lVert
f-f_{n}\right\rVert_{L^{\infty}(\mathcal{M})}^{2}}\right).$ (107)

Moreover, by Lemma 16, since $\gamma\notin\mathbb{Z}_{>0}$, we have

$\left\lVert f-f_{n}\right\rVert_{L^{\infty}(\mathcal{M})}\lesssim\left\lVert
f-f_{n}\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}^{\frac{d}{2\gamma+d}}\left\lVert
f-f_{n}\right\rVert_{L^{2}(p_{0})}^{\frac{2\gamma}{2\gamma+d}}.$ (108)

Using the Sobolev Embedding Theorem, given in [16, Theorem 4], $\left\lVert
f-f_{n}\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}\lesssim\left\lVert
f_{n}\right\rVert_{\mathbb{H}}+\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}\lesssim\tau\sqrt{n}\varepsilon_{n}r^{s}$
whenever $\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}\leq\tau\sqrt{n}\varepsilon_{n}r^{s}$.
Therefore, for a constant $c>0$, when $\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}\leq\tau\sqrt{n}\varepsilon_{n}r^{s}$
and $\varepsilon_{n}r\leq\left\lVert f-f_{n}\right\rVert_{L^{2}(p_{0})}$, we
have that

$\displaystyle\operatorname{\mathbb{E}}_{\boldsymbol{x}}\mathbbold{1}_{\left\lVert
f-f_{n}\right\rVert_{L^{2}(\mathcal{M})}\geq 2\left\lVert
f-f_{n}\right\rVert_{n}}$ $\displaystyle\lesssim\exp\left(-cn\frac{\left\lVert
f-f_{n}\right\rVert_{L^{2}(p_{0})}^{2}}{\left\lVert
f-f_{n}\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}^{\frac{2d}{2\gamma+d}}\left\lVert
f-f_{n}\right\rVert_{L^{2}(p_{0})}^{\frac{4\gamma}{2\gamma+d}}}\right)$ (109)
$\displaystyle\leq
e^{-c\tau^{-\frac{2d}{2\gamma+d}}n^{\frac{2\gamma}{2\gamma+d}}r^{\frac{2d}{2\gamma+d}(1-s)}}.$
(110)

Hence, we can bound the last term as

$\displaystyle\operatorname{\mathbb{E}}_{\boldsymbol{x},\boldsymbol{y}}\int_{0}^{\infty}r\mathbbold{1}_{\mathcal{A}_{r}}\operatorname{\mathbb{P}}_{f\sim\Pi_{n}(\cdot\mid\boldsymbol{x},\boldsymbol{y})}\left(\mathcal{B}^{(\operatorname{III})}(r)\right)\mathop{}\\!\mathrm{d}r$
(111)
$\displaystyle\quad\lesssim\int_{0}^{\infty}re^{(c+1)n\varepsilon_{n}^{2}r^{2s}}e^{-c\tau^{-\frac{2d}{2\gamma+d}}n^{\frac{2\gamma}{2\gamma+d}}r^{\frac{2d}{2\gamma+d}(1-s)}}\mathop{}\\!\mathrm{d}r.$
(112)

We have
$n^{\frac{2\gamma}{2\gamma+d}}=n\left(n^{-\frac{d/2}{2\gamma+d}}\right)^{2}$.
Since $\varepsilon_{n}\lesssim n^{-\frac{\min(\nu,\beta)}{2\nu+d}}$ and
$\min(\nu,\beta)>d/2$, we have $n\varepsilon_{n}^{2}\lesssim
n^{\frac{2\gamma}{2\gamma+d}}$ for some $\gamma\in(d/2,\nu)$. Moreover, for
this choice of $\gamma$ and $s$ small enough we have
$\frac{2d}{2\gamma+d}(1-s)\geq 2s$, which proves that for some constants
$C,C^{\prime},C^{\prime\prime}>0$ the fourth term is bounded by

$C\int_{0}^{\infty}re^{-C^{\prime}r^{C^{\prime\prime}}}\mathop{}\\!\mathrm{d}r<\infty.$
(113)

This concludes the proof. ∎

######  Remark 20.

Following the proof, it is easy to see that $\left\lVert
f-f_{0}\right\rVert_{L^{2}(p_{0})}^{2}$ on the left-hand side of Equations 9,
10 and 12 can be replaced with $\left\lVert
f-f_{0}\right\rVert_{L^{2}(p_{0})}^{q}$ for any $q>1$, changing the exponent
$-\tfrac{2\min(\beta,\nu)}{2\nu+d}$ on the right-hand side to
$-\tfrac{q\min(\beta,\nu)}{2\nu+d}$.

##  Appendix C Characterizing Reproducing Kernel Hilbert Spaces of Matérn
Kernels

We start by describing the reproducing kernel Hilbert spaces (RKHSs) of the
(truncated) intrinsic Matérn processes, proving that they coincide with
(certain subspaces of) the manifold Sobolev spaces. We follow the ideas of
[11], where the same was shown somewhat implicitly. We consider the more-
involved case of the extrinsic Matérn processes immediately after.

The next lemma describes the RKHS of the intrinsic Matérn processes, including
truncated variants. This result is easy to obtain since we have defined them
in terms of the Karhunen–Loève expansions.

######  Lemma 21.

Let $\mathbb{H}_{J}$ be the RKHS of the intrinsic Matérn Gaussian process with
smoothness parameter $\nu$ truncated at the level
$J\in\mathbb{Z}_{>0}\cup\left\\{\infty\right\\}$, and let
$\left\\{f_{j}\right\\}_{j=0}^{\infty}$ be the orthonormal basis of
Laplace–Beltrami eigenfunctions. The space $\mathbb{H}_{J}$ is norm-
equivalent—with constants depending only on $\nu,\kappa$ and
$\sigma_{f}^{2}$—to the set of functions $f=\sum_{j=1}^{J}b_{j}f_{j}$ with
$b_{j}\in\mathbb{R}$, equipped with the inner product

$\Big{\langle}\sum_{j=1}^{J}b_{j}f_{j},\sum_{j=1}^{J}b_{j}^{\prime}f_{j}\Big{\rangle}_{\mathbb{H}_{J}}=\sum_{j=1}^{J}\left(1+\lambda_{j}\right)^{\nu+d/2}b_{j}b_{j}^{\prime}.$
(114)

In particular, $\mathbb{H}_{J}\subset H^{\nu+d/2}\left(\mathcal{M}\right)$ for
all $J$, and for every $h\in\mathbb{H}_{J}$ we have $\left\lVert
h\right\rVert_{\mathbb{H}_{J}}=\left\lVert
h\right\rVert_{H^{\nu+d/2}\left(\mathcal{M}\right)}$.

###### Proof.

By direct computation, the covariance $k$ of the (truncated) intrinsic
Gaussian process is

$k(x,x^{\prime})=\frac{\sigma_{f}^{2}}{C_{\nu,\kappa}}\sum_{j=1}^{J}\left(\frac{2\nu}{\kappa^{2}}+\lambda_{j}\right)^{-(\nu+d/2)}f_{j}(x)f_{j}(x^{\prime}).$
(115)

Hence, the covariance operator $K:L^{2}\left(\mathcal{M}\right)\to
L^{2}\left(\mathcal{M}\right)$ defined by

$(Kf)(x)=\int_{\mathcal{M}}k(x,x^{\prime})f(x^{\prime})\mathop{}\\!\mathrm{d}x^{\prime}$
(116)

is diagonal in the basis $\left\\{f_{j}\right\\}_{j=1}^{J}$, with
$Kf_{j}=\frac{\sigma_{f}^{2}}{C_{\nu,\kappa}}\left(\frac{2\nu}{\kappa^{2}}+\lambda_{j}\right)^{-(\nu+d/2)}f_{j}$.
Then, [28, Theorem 4.2] implies that $\mathbb{H}_{J}$ consists of functions of
form $f=\sum_{j=1}^{J}a_{j}f_{j}$ satisfying

$\left\lVert
f\right\rVert_{\mathbb{H}_{J}}^{2}=\frac{\sigma_{f}^{2}}{C_{\nu,\kappa}}\sum_{j=1}^{J}\left(\frac{2\nu}{\kappa^{2}}+\lambda_{j}\right)^{\nu+d/2}\left|a_{j}\right|^{2}<\infty.$
(117)

Using the elementary inequality
$\min\left(\frac{2\nu}{\kappa^{2}},1\right)\leq\frac{\frac{2\nu}{\kappa^{2}}+\lambda}{1+\lambda}\leq\max\left(\frac{2\nu}{\kappa^{2}},1\right)$,
we find that this space is norm-equivalent to the space
$H^{\nu+d/2}_{J}(\mathcal{M})$ of functions $f=\sum_{j=1}^{J}a_{j}f_{j}$
satisfying

$\left\lVert
f\right\rVert_{H^{\nu+d/2}_{J}(\mathcal{M})}^{2}=\sum_{j=1}^{J}\left(1+\lambda_{j}\right)^{\nu+d/2}\left|a_{j}\right|^{2}<\infty$
(118)

where the comparison constants
$\sqrt{\frac{\sigma_{f}^{2}}{C_{\nu,\kappa}}\min\left(1,\frac{2\nu}{\kappa^{2}}\right)}$
and
$\sqrt{\frac{\sigma_{f}^{2}}{C_{\nu,\kappa}}\max\left(1,\frac{2\nu}{\kappa^{2}}\right)}$
only depend on the parameters $\nu$, $\kappa$, $\sigma_{f}^{2}$. ∎

The next two theorems will be useful to characterize the RKHS of the extrinsic
Matérn process on $\mathcal{M}$. We start by a lemma relating the RKHS of the
restriction of a Gaussian process to the original one.

######  Lemma 22.

Assume that $k$ is a kernel on $\mathbb{R}^{d}$,
$f\sim\operatorname{GP}\left(0,k\right)$ with almost surely continuous sample
paths and $\widetilde{\mathbb{H}}$ is the RKHS of $k$. If
$\mathcal{M}\subseteq\mathbb{R}^{d}$ is a submanifold, then the RKHS
$\mathbb{H}$ corresponding to the restricted process $f_{|\mathcal{M}}$ is the
set of all restrictions $g_{|\mathcal{M}}$ of functions
$g\in\widetilde{\mathbb{H}}$ equipped with the norm

$\left\lVert
h\right\rVert_{\mathbb{H}}=\inf_{g\in\widetilde{\mathbb{H}},\,\,g_{|\mathcal{M}}=h}\left\lVert
g\right\rVert_{\widetilde{\mathbb{H}}}.$ (119)

Moreover there always exists an element $g\in\widetilde{\mathbb{H}}$ such that
$g_{|\mathcal{M}}=f$ and $\left\lVert
g\right\rVert_{\widetilde{\mathbb{H}}}=\left\lVert
f\right\rVert_{\mathbb{H}}$.

###### Proof.

Lemma 5.1 in [63]. ∎

The last result will be used to characterize the RKHS of the extrinsic Matérn
Gaussian processes using trace and extension operators. The second ingredient
for this is the following.

######  Theorem 23.

If $s>\frac{D-d}{2}$ then the restriction operator extends to a bounded linear
map $\operatorname{Tr}:H^{s}\left(\mathbb{R}^{D}\right)\to
H^{s-\frac{D-d}{2}}\left(\mathcal{M}\right)$. Moreover, for every $u>0$ there
exists a bounded right inverse
$\operatorname{Ex}:H^{u}\left(\mathcal{M}\right)\to
H^{u+\frac{D-d}{2}}\left(\mathbb{R}^{D}\right)$ such that
$\operatorname{Tr}\circ\operatorname{Ex}=\operatorname{id}_{H^{u}\left(\mathcal{M}\right)}$
where $\operatorname{Tr}$ corresponds to $s=u+\frac{D-d}{2}$.

###### Proof.

Theorem 4.10 in [24]. ∎

The last two results allow us to characterize the RKHS of the extrinsic Matérn
process on $\mathcal{M}$.

######  Proposition 24.

The RKHS $\mathbb{H}$ of a restricted extrinsic Matérn process $f$ with
smoothness parameter $\nu$ on $\mathcal{M}$ is norm-equivalent to the Sobolev
space $H^{\nu+d/2}\left(\mathcal{M}\right)$.

###### Proof.

Using Lemma 22, the RKHS $\mathbb{H}$ can be characterized as the set of
functions $f:\mathcal{M}\to\mathbb{R}$ that are the restrictions of some
$g\in\widetilde{\mathbb{H}}$, where $\widetilde{\mathbb{H}}$ is the RKHS of
the ambient Matérn process $\tilde{f}$, with

$\left\lVert
f\right\rVert_{\mathbb{H}}=\inf_{g\in\widetilde{\mathbb{H}},\,\,g_{|\mathcal{M}}=f}\left\lVert
g\right\rVert_{\widetilde{\mathbb{H}}}.$ (120)

Since $\widetilde{\mathbb{H}}$ is norm-equivalent to the Sobolev space777Note
that this norm-equivalence is the only property of the Gaussian process we use
in the proofs. Any other Gaussian process satisfying this, including ones
different from the Matérn processes of [11], would also work. This is of
potential interest since other Euclidean kernels, such as Wendland kernels
[61], are known to possess RKHSs which are norm-equivalent to those of the
Matérn kernel. $H^{\nu+D/2}\left(\mathbb{R}^{D}\right)$—see for instance
[28]—by Theorem 23, for every $f\in\mathbb{H}$ we have

$\left\lVert
f\right\rVert_{\mathbb{H}}\lesssim\left\lVert\operatorname{Ex}\left(f\right)\right\rVert_{H^{\nu+D/2}\left(\mathbb{R}^{D}\right)}\lesssim\left\lVert
f\right\rVert_{H^{\nu+D/2-\frac{D-d}{2}}\left(\mathcal{M}\right)}=\left\lVert
f\right\rVert_{H^{\nu+d/2}\left(\mathcal{M}\right)}.$ (121)

Similarly, for every $g\in\widetilde{\mathbb{H}}$ with $g_{|\mathcal{M}}=f$,
we have

$\left\lVert f\right\rVert_{H^{\nu+d/2}\left(\mathcal{M}\right)}=\left\lVert
g_{|\mathcal{M}}\right\rVert_{H^{\nu+d/2}\left(\mathcal{M}\right)}\lesssim\left\lVert
g\right\rVert_{H^{\nu+D/2}\left(\mathbb{R}^{D}\right)}\lesssim\left\lVert
g\right\rVert_{\widetilde{\mathbb{H}}}.$ (122)

Hence, taking the infimum, we obtain

$\left\lVert
f\right\rVert_{H^{\nu+d/2}\left(\mathcal{M}\right)}\lesssim\inf_{g\in\widetilde{\mathbb{H}},\,\,g_{|\mathcal{M}}=f}\left\lVert
g\right\rVert_{\widetilde{\mathbb{H}}}=\left\lVert
f\right\rVert_{\mathbb{H}}.$ (123)

The claim follows. ∎

##  Appendix D Concentration and Sample Path Regularity

In this section, we prove that intrinsic, truncated intrinsic and extrinsic
Matérn processes are sub-Gaussian, uniformly with respect to the truncation
parameter in the case of the truncated intrinsic Matérn process, and live in
Hölder spaces with appropriate exponents. On our way to proving this, we
characterize sample path regularity of (truncated) intrinsic Matérn processes,
which is of independent interest. A simple way to do this would be to build on
the results of Appendix C, using Driscoll’s Theorem—given in [28, Theorem
4.9]—and the Sobolev Embedding Theorem—[16, Theorem 4]—but that would only
give us that the sample paths are almost surely in
$\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)$ for every
$0<\gamma<\nu-d/2,\gamma\notin\mathbb{N}$, whereas here we improve the range
of index to $\gamma<\nu$.

Kolmogorov’s continuity criterion is a standard tool to show that a given
stochastic process has a Hölder continuous version: we re-prove it here
because we will need a form of the result which gives explicit control of the
Hölder norms, which is not usually included in the respective statement.

######  Lemma 25 (Kolmogorov’s continuity criterion).

If $g\sim\Pi$ is a zero-mean Gaussian process on $[0,1]^{d}$ with

$\operatorname{\mathbb{E}}\,\left|g(x)-g(y)\right|^{2}\leq C\left\lVert
x-y\right\rVert_{\mathbb{R}^{d}}^{2\rho}$ (124)

for some $0<\rho\leq 1$ and $C>0$, then there exists a version of $g$ with
samples paths in $\mathcal{CH}^{\alpha}\left([0,1]^{d}\right)$ for every
$0<\alpha<\rho$. Moreover for every $\alpha<\rho$ this version satisfies
$\operatorname{\mathbb{E}}\,\left\lVert
g\right\rVert_{\mathcal{CH}^{\alpha}\left([0,1]^{d}\right)}^{2}\leq
C^{\prime}$ where $C^{\prime}<+\infty$ depends only on $C,\rho,\alpha$ and
$d$.

###### Proof.

Take $x,y\in[0,1]^{d},M>0$ and $q\in\mathbb{N}$. Since the random variable
$g(x)-g(y)$ is Gaussian we have

$\operatorname{\mathbb{E}}\,\left|g(x)-g(y)\right|^{2q}=\frac{\left(2q\right)!}{2^{q}q!}\left(\operatorname{\mathbb{E}}\,\left|g(x)-g(y)\right|^{2}\right)^{q}\leq
C_{q}\left\lVert x-y\right\rVert_{\mathbb{R}^{d}}^{2\rho q}$ (125)

where we have defined $C_{q}=C^{q}\frac{\left(2q\right)!}{2^{q}q!}$. The
reason for considering the the $2q$th power will become clear later in the
proof. By Markov’s inequality, for every $x,y\in[0,1]^{d}$ we have

$\operatorname{\mathbb{P}}\left(\left|g(x)-g(y)\right|>u\right)\leq
C_{q}u^{-2q}\left\lVert x-y\right\rVert_{\mathbb{R}^{d}}^{2q\rho}.$ (126)

Now, take $X=\cup_{k\geq 0}X_{k},X_{k}=2^{-k}\mathbb{Z}^{d}\cap[0,1]^{d}$.
Then, the previous inequality applied to any adjacent $x,y\in X_{k}$, where we
see $X_{k}$ as a graph where two vertices are connected if they differ by
$2^{-k}$, and $u=M2^{-k\alpha}$, implies

$\operatorname{\mathbb{P}}\left(\left|g(x)-g(y)\right|>M2^{-k\alpha}\right)\leq
C_{q}M^{-2q}2^{-2kq(\rho-\alpha)}.$ (127)

Denote $K=\sup_{x,y\in
X\mathrel{\text{missing}}{adjacent}}\frac{\left|g(x)-g(y)\right|}{\left\lVert
x-y\right\rVert^{\alpha}}$. Summing over $k\geq 1$ and adjacent points in
$X$—and there are at most $C2^{k(d+1)}$ of them where $C>0$ is an absolute
constant—gives us for $q>\frac{d+1}{2(\rho-\alpha)}$ that

$\displaystyle\operatorname{\mathbb{P}}\left(K>M\right)$
$\displaystyle\leq\sum_{k\geq 0}\sum_{x,y\in
X\mathrel{\text{missing}}{adjacent}}\operatorname{\mathbb{P}}\left(\left|g(x)-g(y)\right|>M2^{-k\alpha}\right)$
(128) $\displaystyle\leq C\sum_{k\geq
0}2^{k(d+1)}C_{q}M^{-2q}2^{-2kq(\rho-\alpha)}=\frac{CC_{q}}{1-2^{2q(\rho-\alpha)-d-1}}M^{-2q}.$
(129)

In particular, for all $q>\max\left(1,\frac{d+1}{2(\rho-\alpha)}\right)$, we
have

$\displaystyle\operatorname{\mathbb{E}}(K^{2})$
$\displaystyle=2\int_{0}^{\infty}\\!\\!M\operatorname{\mathbb{P}}\left(K>M\right)\mathop{}\\!\mathrm{d}M=2\int_{0}^{1}\\!\\!M\operatorname{\mathbb{P}}\left(K>M\right)\mathop{}\\!\mathrm{d}M+2\int_{1}^{\infty}\\!\\!M\operatorname{\mathbb{P}}\left(k>M\right)\mathop{}\\!\mathrm{d}M$
(130) $\displaystyle\leq
2+2\int_{1}^{\infty}M\operatorname{\mathbb{P}}\left(K>M\right)\mathop{}\\!\mathrm{d}M\leq
C_{C,\alpha,\rho,d}$ (131)

for some constant $C_{C,\alpha,\rho,d}<+\infty$. This means that $K$ is finite
almost surely. Since $X$ is dense in $[0,1]^{d}$ and $g$ is almost surely
uniformly continuous on $X$ by the classical version of the Kolmogorov’s
continuity criterion—see for instance [34]—$g$ admits a unique continuous
extension to $[0,1]^{d}$ on a probability one event $\mathcal{A}$. Let us
define

$\displaystyle\overline{g}(x)=\begin{cases}\displaystyle\lim_{y\to x,y\in
X}g(y)&\mathrel{\text{missing}}{on}\mathcal{A},\\\
0&\mathrel{\text{missing}}{otherwise},\end{cases}$ (132)

for $x\in[0,1]^{d}$. For any $x,y\in[0,1]^{d}$ and $x_{n}\to x,y_{n}\to
y,x_{n},y_{n}\in X$ we have on $\mathcal{A}$

$\displaystyle\left|\overline{g}(x)-\overline{g}(y)\right|\leq$
$\displaystyle\liminf_{n\to\infty}\left(\left|\overline{g}(x)-\overline{g}(x_{n})\right|+\left|\overline{g}(x_{n})-\overline{g}(y_{n})\right|+\left|\overline{g}(y_{n})-\overline{g}(y)\right|\right)$
(133) $\displaystyle\leq$
$\displaystyle\liminf_{n\to\infty}\left(\left|\overline{g}(x)-\overline{g}(x_{n})\right|+K\left\lVert
x_{n}-y_{n}\right\rVert_{\mathbb{R}^{d}}^{\alpha}+\left|\overline{g}(y_{n})-\overline{g}(y)\right|\right)$
(134) $\displaystyle=$ $\displaystyle K\left\lVert
x-y\right\rVert_{\mathbb{R}^{d}}^{\alpha}.$ (135)

On the complement of $\mathcal{A}$, we have $\overline{g}=0$. Hence,
$\overline{g}$ is $\alpha$-Hölder continuous on $[0,1]^{d}$ with the same
(random) constant $K$, since

$\operatorname{\mathbb{E}}\,\left\lVert\overline{g}\right\rVert_{\mathcal{CH}^{\alpha}\left([0,1]^{d}\right)}^{2}=\operatorname{\mathbb{E}}\,\left(\sup_{x,y\in
X}\frac{\left|\overline{g}(x)-\overline{g}(y)\right|}{\left\lVert
x-y\right\rVert^{\alpha}}\right)^{2}\leq\operatorname{\mathbb{E}}K^{2}\leq
C_{C,\alpha,\rho,d}<\infty.$ (136)

Since $\overline{g}$ is a version of $g$, the claim follows. ∎

The next lemma applies our version of Kolmogorov’s criterion, Lemma 25, to the
intrinsic Matérn processes on $\mathcal{M}$ by considering charts. This allows
us to show that the sample paths are almost surely in
$\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)$ for every $\gamma<\nu$, which
is both used in our arguments and also of independent interest. For the claims
in Appendix B, we need to ensure that this property holds somewhat uniformly
with respect to the truncation parameter, which is why we tracked the
constants in our proof of Kolmogorov’s criterion. As we will see, the main
difficulty in the proof of the next result will be to tackle the case of
regularity strictly larger than $1$.

######  Lemma 26.

Let $f\sim\Pi_{n}$ be an intrinsic Matérn process with smoothness parameter
$\nu>0$ truncated at $J_{n}\in\mathbb{N}\cup\left\\{\infty\right\\}$. Then for
every $\gamma<\nu$ we have

$\sup_{n}\operatorname{\mathbb{E}}_{f\sim\Pi_{n}}\left[\left\lVert
f\right\rVert_{\mathcal{CH}^{\gamma}\left(\mathcal{M}\right)}^{2}\right]<\infty.$
(137)

###### Proof.

We start with the case $\nu\leq 1$. Take $1\leq l\leq L$ and define
$h_{l}=\left(\chi_{l}f\right)\circ\phi_{l}^{-1}$. Then $h_{l}$ is a Gaussian
process with covariance kernel $\widetilde{K}_{l}$ given by

$\widetilde{K}_{l}(x,y)=(\chi_{l}\circ\phi_{l}^{-1})(x)K(\phi_{l}^{-1}(x),\phi_{l}^{-1}(y))(\chi_{l}\circ\phi_{l}^{-1})(y)$
(138)

where $x,y\in\mathcal{V}_{l}$ and
$K(x,y)=\operatorname{Cov}\left(f(x),f(y)\right)$ is the covariance kernel of
$f$. Let $\widetilde{\mathbb{H}}_{l}$ be the RKHS induced by
$\widetilde{K}_{l}$. We seek to apply Lemma 25 to $h_{l}$. For all
$x,y\in\mathcal{V}_{l}$, where we recall that we can assume that
$\mathcal{V}_{l}=(0,1)^{d},$, we have

$\displaystyle\operatorname{\mathbb{E}}_{f\sim\Pi_{n}}\left|h_{l}(x)-h_{l}(y)\right|^{2}$
$\displaystyle=\widetilde{K}_{l}(x,x)+\widetilde{K}_{l}(y,y)-2\widetilde{K}_{l}(x,y)$
(139)
$\displaystyle=\left\lVert\widetilde{K}_{l}(x,\cdot)-\widetilde{K}_{l}(y,\cdot)\right\rVert_{\widetilde{\mathbb{H}}_{l}}^{2}$
(140)
$\displaystyle=\sup_{\left\lVert\varphi\right\rVert_{\widetilde{\mathbb{H}}_{l}}=1}\left|\left\langle\widetilde{K}_{l}(x,\cdot)-\widetilde{K}_{l}(y,\cdot),\varphi\right\rangle_{\widetilde{\mathbb{H}}_{l}}\right|^{2}$
(141)
$\displaystyle=\sup_{\left\lVert\varphi\right\rVert_{\widetilde{\mathbb{H}}_{l}}=1}\left|\varphi(x)-\varphi(y)\right|^{2}$
(142)
$\displaystyle\leq\sup_{\left\lVert\varphi\right\rVert_{\widetilde{\mathbb{H}}_{l}}=1}\left\lVert\varphi\right\rVert_{\mathcal{CH}^{\nu}\left(\mathcal{V}_{l}\right)}^{2}\left\lVert
x-y\right\rVert_{\mathbb{R}^{d}}^{2\nu}$ (143)

where
$\left\lVert\varphi\right\rVert_{\mathcal{CH}^{\nu}\left(\mathcal{V}_{l}\right)}^{2}$
can potentially be infinite or unbounded: we will show this is not the case.
To do this it suffices to show that we have a continuous embedding
$\widetilde{\mathbb{H}}_{l}\hookrightarrow\mathcal{CH}^{\nu}\left(\mathcal{V}_{l}\right)$,
that is
$\left\lVert\cdot\right\rVert_{\mathcal{CH}^{\nu}\left(\mathcal{V}_{l}\right)}\lesssim\left\lVert\cdot\right\rVert_{\widetilde{\mathbb{H}}_{l}}$.
The RKHS $\widetilde{\mathbb{H}}_{l}$ is by definition the completion of

$\displaystyle\left\\{\sum_{i=1}^{p}\alpha_{i}\widetilde{K}_{l}\left(x_{i},\cdot\right):p\geq
1,\alpha_{i}\in\mathbb{R},x_{i}\in\mathcal{V}_{l}\right\\}$ (144)
$\displaystyle=\left\\{\sum_{i=1}^{p}\alpha_{i}\left(\chi_{l}\circ\phi_{l}^{-1}\right)(x_{i})\left(\chi_{l}\circ\phi_{l}^{-1}\right)(\cdot)K\left(\phi_{l}^{-1}(x_{i}),\phi_{l}^{-1}(\cdot)\right):p\geq
1,\alpha_{i}\in\mathbb{R},x_{i}\in\mathcal{V}_{l}\right\\}$ (145)

with respect to the topology induced by the RKHS norm

$\left\lVert\sum_{i=1}^{p}\alpha_{i}\widetilde{K}_{l}\left(x_{i},\cdot\right)\right\rVert_{\widetilde{\mathbb{H}}_{l}}^{2}=\sum_{i,j=1}^{p}\alpha_{i}\alpha_{j}\left(\chi_{l}\circ\phi_{l}^{-1}\right)(x_{i})\left(\chi_{l}\circ\phi_{l}^{-1}\right)(x_{j})K\left(\phi_{l}^{-1}(x_{i}),\phi_{l}^{-1}(x_{j})\right).$
(146)

Denote the RKHS of $K$ by $\mathbb{H}$. By Theorem 13, and by the equality
$\left\lVert\cdot\right\rVert_{\mathbb{H}}=\left\lVert\cdot\right\rVert_{H^{\nu+d/2}\left(\mathcal{M}\right)}$
on $\mathbb{H}$ which follows by Lemma 21, we have

$\displaystyle\left\lVert\sum_{i=1}^{p}\alpha_{i}\widetilde{K}_{l}(x_{i},\cdot)\right\rVert_{H^{\nu+d/2}\left(\mathbb{R}^{d}\right)}^{2}$
(147)
$\displaystyle\quad=\left\lVert\sum_{i=1}^{p}\alpha_{i}\left(\chi_{l}\circ\phi_{l}^{-1}\right)(x_{i})\left(\chi_{l}\circ\phi_{l}^{-1}\right)(\cdot)K(\phi_{l}^{-1}(x_{i}),\phi_{l}^{-1}(\cdot))\right\rVert_{H^{\nu+d/2}\left(\mathbb{R}^{d}\right)}^{2}$
(148)
$\displaystyle\quad\lesssim\left\lVert\sum_{i=1}^{p}\alpha_{i}\left(\chi_{l}\circ\phi_{l}^{-1}\right)(x_{i})K(\phi_{l}^{-1}(x_{i}),\cdot)\right\rVert_{H^{\nu+d/2}\left(\mathcal{M}\right)}^{2}$
(149)
$\displaystyle\quad=\left\lVert\sum_{i=1}^{p}\alpha_{i}\left(\chi_{l}\circ\phi_{l}^{-1}\right)(x_{i})K(\phi_{l}^{-1}(x_{i}),\cdot)\right\rVert_{\mathbb{H}}^{2}$
(150)
$\displaystyle\quad=\sum_{i,j=1}^{p}\alpha_{i}\alpha_{j}\left(\chi_{l}\circ\phi_{l}^{-1}\right)(x_{i})\left(\chi_{l}\circ\phi_{l}^{-1}\right)(x_{j})K\left(\phi_{l}^{-1}(x_{i}),\phi_{l}^{-1}(x_{j})\right)$
(151)
$\displaystyle\quad=\left\lVert\sum_{i=1}^{p}\alpha_{i}\widetilde{K}_{l}(x_{i},\cdot)\right\rVert_{\widetilde{\mathbb{H}}_{l}}^{2}.$
(152)

Therefore, we have a continuous embedding
$\widetilde{\mathbb{H}}_{l}\hookrightarrow
H^{\nu+d/2}\left(\mathbb{R}^{d}\right)$ with
$\left\lVert\cdot\right\rVert_{H^{\nu+d/2}\left(\mathbb{R}^{d}\right)}\lesssim\left\lVert\cdot\right\rVert_{\widetilde{\mathbb{H}}_{l}}$
on $\widetilde{\mathbb{H}}_{l}$. By the Sobolev Embedding Theorem in
$\mathbb{R}^{d}$—see for instance [50, Section 2.7.1, Remark 2]—we have
$B_{2,2}^{\nu+d/2}\left(\mathbb{R}^{d}\right)=H^{\nu+d/2}\left(\mathbb{R}^{d}\right)\hookrightarrow\mathcal{CH}^{\nu}\left(\mathbb{R}^{d}\right)$,
which implies
$\widetilde{\mathbb{H}}_{l}\hookrightarrow\mathcal{CH}^{\nu}\left(\mathbb{R}^{d}\right)$
by composition. Thus, there exists a constant $C=C_{\nu}$ such that

$\operatorname{\mathbb{E}}_{f\sim\Pi_{n}}\left|h_{l}(x)-h_{l}(y)\right|^{2}\leq
C\left\lVert x-y\right\rVert_{\mathbb{R}^{d}}^{2\nu}$ (153)

for $x,y\in\mathcal{V}_{l}$. Hence, by applying Lemma 25, there exists a
version $\tilde{h}_{l}$ of $h_{l}$ with almost surely $\alpha$-Hölder
continuous sample paths for every $\alpha<\nu$. Now consider
$\tilde{f}=\sum_{l=1}^{L}\tilde{h}_{l}\circ\phi_{l}$. Then $\tilde{f}$ is a
version of $f$. We proceed to bound
$\operatorname{\mathbb{E}}_{f\sim\Pi_{n}}\left[\lVert\tilde{f}\rVert_{\mathcal{CH}^{\alpha}\left(\mathcal{M}\right)}^{2}\right]$.
For any $1\leq l,r\leq L$ write

$\left|\tilde{h}_{r}(\phi_{r}(\phi_{l}^{-1}(x)))-\tilde{h}_{r}(\phi_{r}(\phi_{l}^{-1}(y)))\right|\leq
K\left|\phi_{r}(\phi_{l}^{-1}(x))-\phi_{r}(\phi_{l}^{-1}(y))\right|^{\alpha}\leq
CK\left|x-y\right|^{\alpha}$ (154)

where $C$ is the Lipshitz constant of $\phi_{r}\circ\phi_{l}^{-1}$ which is
well defined and finite because this composition is a diffeomorphism and $K$
is a random constant with $\operatorname{\mathbb{E}}_{f\sim\Pi_{n}}K^{2}\leq
C_{\alpha,\nu,d}$. Hence

$\displaystyle\operatorname{\mathbb{E}}_{f\sim\Pi_{n}}\big{\lVert}\tilde{f}\big{\rVert}_{\mathcal{CH}^{\alpha}\left(\mathcal{M}\right)}^{2}$
$\displaystyle=\operatorname{\mathbb{E}}_{f\sim\Pi_{n}}\sum_{l=1}^{L}\left\lVert\left(\chi_{l}\tilde{f}\right)\circ\phi_{l}^{-1}\right\rVert_{\mathcal{CH}^{\alpha}\left(\mathbb{R}^{d}\right)}^{2}$
(155)
$\displaystyle=\operatorname{\mathbb{E}}_{f\sim\Pi_{n}}\sum_{l=1}^{L}\left\lVert\left(\chi_{l}\sum_{r=1}^{L}\tilde{h}_{r}\circ\phi_{r}\right)\circ\phi_{l}^{-1}\right\rVert_{\mathcal{CH}^{\alpha}\left(\mathbb{R}^{d}\right)}^{2}$
(156)
$\displaystyle\lesssim\sum_{l=1}^{L}\sum_{r=1}^{L}\operatorname{\mathbb{E}}_{f\sim\Pi_{n}}\left\lVert\tilde{h}_{r}\circ\phi_{r}\circ\phi_{l}^{-1}\right\rVert_{\mathcal{CH}^{\alpha}\left(\mathbb{R}^{d}\right)}^{2}\lesssim\operatorname{\mathbb{E}}_{f\sim\Pi_{n}}K^{2}\leq
C_{\alpha,\nu,d}.$ (157)

which gives the $\nu\leq 1$ case.

We now turn to the general case. The proof will be similar to the one of [20,
Proposition I.3] although we need to control the Hölder norms and work through
charts since our Gaussian processes are supported on manifolds. Assume for
simplicity that $d=1,1<\nu\leq 2$, otherwise it suffices to introduce
coordinates and to proceed by induction on $\left\lfloor\nu\right\rfloor$. Let
$l\in\left\\{1,\ldots,L\right\\}$, and as before define
$\widetilde{K}_{l}\left(x,y\right)=\left(\chi_{l}\circ\phi_{l}^{-1}\right)(x)\left(\chi_{l}\circ\phi_{l}^{-1}\right)(y)K\left(\phi_{l}^{-1}(x),\phi_{l}^{-1}(y)\right)$
the covariance kernel of $h_{l}=\left(\chi_{l}f\right)\circ\phi_{l}^{-1}$ as
well as $\widetilde{\mathbb{H}}_{l}$ its RKHS.

First, let us construct an $L^{2}(\Omega)$-derivative $\dot{h}_{l}$ of
$h_{l}$, where $L^{2}(\Omega)$ is the space of random variables with finite
variance with $\left\langle
a,b\right\rangle_{L^{2}(\Omega)}=\operatorname{\mathbb{E}}(ab)$. This
derivative is a square integrable process on $\mathcal{V}_{l}$ such that

$\operatorname{\mathbb{E}}_{f\sim\Pi_{n}}\left|\frac{h_{l}(x+\delta)-h_{l}(x)}{h}-\dot{h}_{l}(x)\right|^{2}\to
0$ (158)

as $h\to 0$, for all $x\in\mathcal{V}_{l}$. For this, we will first show that
$\frac{\partial\widetilde{K}_{l}}{\partial
x}(x,\cdot)\in\widetilde{\mathbb{H}}_{l}$ for every $x\in\mathcal{V}_{l}$—here
$\frac{\partial\widetilde{K}_{l}}{\partial x}$ denotes the derivative of the
function $\widetilde{K}_{l}(\cdot,\cdot^{\prime})$ with respect to the first
argument—and that

$\left\lVert\frac{\partial\widetilde{K}_{l}}{\partial
x}(x,\cdot)-\frac{\partial\widetilde{K}_{l}}{\partial
x}(x^{\prime},\cdot)\right\rVert_{\widetilde{\mathbb{H}}_{l}}\leq
C_{\nu}\left|x-x^{\prime}\right|^{\nu-1}.$ (159)

We first show that
$\frac{\widetilde{K}_{l}(x+\delta,\cdot)-\widetilde{K}_{l}(x,\cdot)}{h}$ is a
Cauchy net888See for instance [1] for a review of Cauchy nets. in
$\widetilde{\mathbb{H}}_{l}$. We have

$\displaystyle\left\lVert\frac{\widetilde{K}_{l}(x+\delta,\cdot)-\widetilde{K}_{l}(x,\cdot)}{h}-\frac{\widetilde{K}_{l}(x+\delta^{\prime},\cdot)-\widetilde{K}_{l}(x,\cdot)}{h^{\prime}}\right\rVert_{\widetilde{\mathbb{H}}_{l}}$
(160)
$\displaystyle\quad=\sup_{\left\lVert\varphi\right\rVert_{\widetilde{\mathbb{H}}_{l}}=1}\left\langle\frac{\widetilde{K}_{l}(x+\delta,\cdot)-\widetilde{K}_{l}(x,\cdot)}{h}-\frac{\widetilde{K}_{l}(x+\delta^{\prime},\cdot)-\widetilde{K}_{l}(x,\cdot)}{h^{\prime}},\varphi\right\rangle_{\widetilde{\mathbb{H}}_{l}}$
(161)
$\displaystyle\quad=\sup_{\left\lVert\varphi\right\rVert_{\widetilde{\mathbb{H}}_{l}}=1}\left(\frac{\varphi(x+\delta)-\varphi(x)}{h}-\frac{\varphi(x+\delta^{\prime})-\varphi(x)}{h^{\prime}}\right)$
(162)
$\displaystyle\quad=\sup_{\left\lVert\varphi\right\rVert_{\widetilde{\mathbb{H}}_{l}}=1}\int_{0}^{1}\left[\varphi^{\prime}(x+th)-\varphi^{\prime}(x+th^{\prime})\right]\mathop{}\\!\mathrm{d}t$
(163)
$\displaystyle\quad\leq\sup_{\left\lVert\varphi\right\rVert_{\widetilde{\mathbb{H}}_{l}}=1}\left\lVert\varphi^{\prime}\right\rVert_{\mathcal{CH}^{\nu-1}\left(\mathcal{V}_{l}\right)}\left|h-h^{\prime}\right|^{\nu-1}$
(164)
$\displaystyle\quad\leq\sup_{\left\lVert\varphi\right\rVert_{\widetilde{\mathbb{H}}_{l}}=1}\left\lVert\varphi\right\rVert_{\mathcal{CH}^{\nu}\left(\mathcal{V}_{l}\right)}\left|h-h^{\prime}\right|^{\nu-1}$
(165)

where in (163) the derivative $\varphi^{\prime}$ exists because, exactly as in
the case $\nu\leq 1$, we can show show that
$\widetilde{\mathbb{H}}_{l}\hookrightarrow\mathcal{CH}^{\nu}\left(\mathbb{R}^{d}\right)$.
This also implies that for a constant $C=C_{\nu}$ we have

$\left\lVert\frac{\widetilde{K}_{l}(x+\delta,\cdot)-\widetilde{K}_{l}(x,\cdot)}{h}-\frac{\widetilde{K}_{l}(x+\delta^{\prime},\cdot)-\widetilde{K}_{l}(x,\cdot)}{h^{\prime}}\right\rVert_{\widetilde{\mathbb{H}}_{l}}\leq
C\left|h-h^{\prime}\right|^{\nu-1}.$ (166)

As $\left|h-h^{\prime}\right|^{\nu-1}\to 0$ when $h,h^{\prime}\to 0$, because
$\nu>1$, this proves that
$\frac{\widetilde{K}_{l}(x+\delta,\cdot)-\widetilde{K}_{l}(x,\cdot)}{h}$ is a
Cauchy net in $\widetilde{\mathbb{H}}_{l}$: by completeness of
$\widetilde{\mathbb{H}}_{l}$ it converges in $\widetilde{\mathbb{H}}_{l}$ to a
limit. Since by general properties of RKHSs, convergence in
$\widetilde{\mathbb{H}}_{l}$ implies pointwise convergence, the limit
satisfies

$\lim_{h\to
0}\frac{\widetilde{K}_{l}(x+\delta,y)-\widetilde{K}_{l}(x,y)}{h}=\frac{\partial\widetilde{K}_{l}}{\partial
x}(x,y).$ (167)

Hence the partial derivative $\frac{\partial\widetilde{K}_{l}}{\partial
x}(x,y)$ exists for all $y$ and $\frac{\partial\widetilde{K}_{l}}{\partial
x}(x,\cdot)\in\widetilde{\mathbb{H}}_{l}$. Moreover, by the isometry
$L^{2}(\Omega)\ni
h_{l}(x)\mapsto\operatorname{\mathbb{E}}_{f\sim\Pi_{n}}h_{l}(x)h_{l}(\cdot)=\widetilde{K}_{l}(x,\cdot)\in\widetilde{\mathbb{H}}_{l}$,