Posterior Contraction Rates for Matérn Gaussian Processes on Riemannian Manifolds
Authors: Paul Rosa
Viacheslav Borovitskiy
Alexander Terenin
Judith Rousseau
Venue: NeurIPS 2023 spotlight (2023)
Keywords: Gaussian processes
posterior contraction
manifolds
kernels

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.

Posterior Contraction Rates for Matérn Gaussian Processes on Riemannian Manifolds
Paul Rosa University of Oxford
Viacheslav Borovitskiy
ETH Zürich
Alexander Terenin
Judith Rousseau University Oxford
University of Cambridge and Cornell University
Abstract
Gaussian processes are used in many machine learning applications that rely on uncertainty quantification. Recently, computational tools for working with these models 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 Rd 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 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 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 work
intrinsic processes can achieve better performance practice. Therefore, 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-asymin asympttt behavior.
1 Introduction
Gaussian processes provide powerful way to quantify uncertainty about unknown regression functions via the formulation of Bayesian learning. Motivated by applications in the physical and engineering sciences, 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, a robotics setting. Jaquier al. [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 Coveney et al. [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 algo- rithms have good theoretical properties, as well as their limitations. Within the Bayesian framework, natural way to quantify data-efficiency and generalization enror is to posit data-generating mech-
37th Conference on Neural Information Processing Systems (NeurIPS 2023).

---

(a) Intrinsic
(b) Extrinsic
(c) Intrinsic
(d) Extrinsic
Figure 1: Samples from different Matérn Gaussian processes on different manifolds, namely one-dimensional dumbbell-shaped manifold and 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-gove global approximation to geodesic distances on the sphere.
anism model and study in-and how fast-the posterior distribution concentrates around the true regression function 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 higher-dimensional Euclidean space. For example, the two-dimensional sphere can be embedded into the Euclidean space R intrinsic Gaussian processes model functions on the sphere while extrinsic ones model functions on R&, which are then restricted to the sphere. Are the former more efficient than the latter? Since embeddings even isometric ones distances locally, they can induce spurious dependencies, as points can be close in the ambient space but far away with respect 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 quantify such differences.
best only preserve
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 common for perceptual data such text and images concentrate pixel space or sequence space. therefore also interesting to investigate how Gaussian process being kemel-based based, simpler than for instance deep neural networks perffsrmerom in such scenarios, at least in the asymptotic regime.
diwer-dimeniional submanifold within, for instance,
In this work, we develop geometric analogs of the Gaussian process posterior contraction theorems of van der Vaart and van Zanten [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 sundomly-sample input points commonly referred to random design in the literature- and contraction in the sense of the (po) distance, defined in Section We focus on compact Riemannian manifolds: this allows one define Matérn Gaussian processes through their Karhunen Loevee expansions, which requires discrete spectrum for the Laplace-Beltrami Beltrami operator-sec for instance Borovitskiy et al. [11] and Chavel [13], Chapter l-and is 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-thoug! 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, unite-sample considerations are necessary. We therefore present experiments that compute the worst case emoor numerically. These experiments highlight that intrinsic models are capable of achieving better performance in the small-data regime. We conclude with discussion of why these shulhh-which might at first seem counterintuitive- 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 1/i = [(:) + Et. with & ~ N(0, o?), J; X, and f is assigned a Gaussian process prior. A Gaussian process is a random function S x + R for which all finiteddimensional marginal distributions are multivariate Gaussian. The distribution of such process is uniquely determined by its mean function ml-) (f)(-) and covariance kernel k(., -') Cov(f(), hence we write f ~ GP k).
For Gaussian process regression, the posterior distribution given the data is also a Gaussian process with probability kernel II(- x. y) = (M(M (a..)) (m.)))).)))) see Rasmussen and Williams [41].
ma M) (.) (.)ae ( Rex + +o2
(1)
KTT .. (:, = K( K (o) (Ka + -o 'K*
(2)
These quantities deseribe 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|
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 posterior asymptotics is concerned with understanding conditions under which this does or does not occur, with questions of posterior contraction rales- how fast such convergence occurs-beinn of key interest. At present, there well-developed literature posterior contraction rates, see Ghosal and van der Vaart [20] for 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
yi x; N(fo(z;),
34 ~ Po
where So E F C R-X, class of real-valued functions, and N(2, o²) denotes the Gaussian with moments μe o². 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 o² in both cases: in this paper, we focus on the particular case where ag is known in advance. This setting is restrictive, one can extend an unknown dg > using techniques that are not specific to our geometric setting: for instance, the approach of [55] allows to handle an unknown de 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 de can be used, such as a conjugate one like in Banerjee [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 T,2 loss given by
1/2
(1 fo(x)) (o(=))(dpo()))
I - (ollelaa(p))
which is arguably the natural way of measuring discrepancy between f and for given the fact that the covariates W; are sampled from Po. The posterior contraction rate is then defined as
E EI (w. fI fo (0
where Ey- (m (-) denotes expectation under the posterior distribution while (..y(-) denotes expectation under the true data generating process." In the case of covariates distributed on [o, 17d, posterior contraction rates have been derived under Matérn Gaussian process priors [47] in van der Vaart and van Zanten [56], who showed the following result.
Note that other notions of posterior contraction can be found in the literature, see Ghosal and van der Vaart [20] and Rousseau [42] for examples that are slightly weaker than the definition work with.
3

---

Result 1 (Theorem 2 of van der Vaart and van Zanten [56]). In the Bayesian regression model, let be a mean-zero Matérn Gaussian process prior on R' with amplitude of length scale ii, and smoothness v > d/2.. Assume that the true data generating process is given by (3), where po has a Lebesgue density on X = [0, 1] which bounded from below and above by 0 < Cpo < Cpo < co. respectively. Let fo E H CCH with 3 > d/2, where H and CH the Sobolev and Hölder spaces, respectively. Then there exists constant C > O, which does not depend on n but does depend on d, off U, R, B. Po, o², (aMy and |I [l (My such that
m)) 20+d
Buyy (e IIf sull S SOn
(6)
and, moreover, the posterior mean satisfies
2 2min(3.v)
Ex.y |m..) folli2(m0) I On
2x+d
(7)
Note that ( is the Bayes estimator [52] of associated the weighted L 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 mn(lr.y) The best rate is attained when 3 v: that is, when true smoothness and prior smoothness match-which is known to be minimax optimal in the problem of estimating fo: see Tsybakov [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 Schwartz [46] and Le Cam [29], but has been extensively developed since the seminal paper of Ghosal et al. [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].
The results closest to ours are those of Yang and Dunson [63] and Castillo et al. [12]. In the former, the authors use an extrinsic length-scale-nn scale-mixture of squared exponential Gaussian processes to achieve optimal contraction rates with respect to the weighted L norm, using completely different proof technique compared to us, and their results are restricted to So having Hölder smoothness of order less than or equal to two. On the other hand Castillo et al. [12] consider, as an intrinsic process on the manifold, a hicrarchical Gaussian process based on its heat kernel and provide posterior contraction rates. For the Matérn class, Li et al. [30] presents results which characterize the asymptotic behavior of kernel nyperparameters: 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: Dunson et al. [17] and Sanz-Alonso and Yang [45] present posterior contraction rates for specific graph Gaussian process model in 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 M, to emphasize that the covariate space is a manifold.
Assumption 2. Assume that M C. R a smooth, compuct submanifold (without boundary) of dimension d < D equipped with the standard Riemannian volume measure p.
We denote | = J du(xx) for volume of M. With this geometric setting defined, we will need to describe regularity assumptions in terms of functional spaces on the manifold M. We work with Hölder spaces CHM)), defined using charts via the usual Euclidean Hölder spaces, the Sobolev spaces H²M (M), and Besov spaces B (M) which are of the ways of generalizing

---

(a) Intrinsic
(b) Extrinsic
(c) Intrinsic
(d) Extrinsic
Figure 2: Different Matérn kemels k:( x) on different manifolds.
the Euclidean Hölder spaces of smooth functions to manifolds. We follow Coulhon et al. [14] and Castillo [12], and define these spaces using the Laplace-Beltrame Beltrami operator M in Appendix A.
Recall that the data-genent lata-generating process given by (3), with fo as the true regression function and po the distribution of the covariates.
Assumption& 3. Assume that Po is absolutely continuous with respect to μ., and that its density, denoted by po, satisfies c S Po S for 0 < C, C < cc. Assume the regression function fo satisfies fo H (M) n B (MM for some 3 2/2, and that o² > 0 is fixed and known.
This setting can be extended to handle unknown variance oz by putting a prior on de, following the strategy of Salomond [44] and Naulet and Barat [36]. Since we focused primarily 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.
3.1 Intrinsic Matérn Gaussian Processes
We now introduce the first geometrie Gaussian process prior under study- the Riemannian Matérn kernel of Whittle [62]. Lindgren et al. [33], and Borovitskiy et al. [11]. This process was originally defined using stochastic partial differential equations: here, present it by its Lovve expansion, to facilitate comparisons with its truncated analogs presented Section 3.2.
Definition 4 (Intrinsic Matérn prior). Let v > and let (Ag, f1);20 be the eigenvalues and orthonor- mal eigenfunctions of the Laplace-Beltrami Beltrami operator M, in increasing order. Define the intrinsic RIEMANNIAN MATÉRN GAUSSIAN PROCESS through its Karhunen-Loeve expansion to be
of 10)= Cu.
Σ j=1
2v 12 +>g
E,5s()
Ej N(O,1) 1)
(8)
where 07 are positive parameters and Cux is the normalization constant, chosen such that [M] JM Var(f(x)) du(x)) ozt. where Var denotes the variance.
1
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 the regularity of the ground -tuth function matches the regularity of the Gaussian process, as described by the parameter v.
Theorem Let f be a Riemannian Matérn Gaussian process prior of Defimition 4 with smoothness parameter d/2 and let fo satisfy Assumption 3. Then there is C 0 such that
2 min((..))
Em I foll2A1po) S Cn
204-
All proofs are given in Appendix B. Our proof follows the general approach of van der Vaart and van Zanten [56], by first proving a contraction rate with respect to the distance n [=1/2]]]]()) fo(æ) [) input locations æ, and then extending the result to true L,--distance by applying suitable concentration inequality. The first part is obtained by studying the concentration function, funitioon which known be the key quantity to control in order derive contraction rates Gaussian process priors-sec Ghosal and van der Vaart [20] and van der Vaart and van Zanten [57] for an overview.
5

---

Given our regularity assumptions on fu, the most difficult part lies in controlling the small- ball probabilities 11 [Illffll < E| 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 related to the Sobolev space (=+ 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 (po) norm, following van der Vaart and van Zanten [56]. we prove a Hölder-type property for manifold Matérn processes, and apply Bernstein's 's inequality. Together, this gives the claim.
This result is good news for the intrinsic Matérn model: it tells 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 fo- An inspection of the proof shows that the constant C' > 0 can be seen to depend on d , V, K, G. po.o". Il (M))) and [[follc**( Theorem 5 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 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 Bechumenn-Lovvee expansion by taking the first J terms, which is also optimal in an L (MM)-sense.
Note that the truncated prior is a randomly-weighted finite sum of Laplace- Laplace-Beltrami eigenfunctions, which have different smoothness properties compared to the original prior: the truncated prior takes its values in C (M) since the eigenfunctions of M are smooth-se See for instance De Vito et al. [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 function with essentially-finie regularity in the limit.
Truncated random basis expansions have been studied extensively in the Bayesian literature in the Euclidean setting-see for instance Arbel et al. [2] and Yoo et al. [64] or Ghosal and van der Vaart [20], Chapter 11 for examples with priors based on wavelet expansions. It is known that truncating the expansion at high enough level usually allows one to retain optimality. Instead of truncating deterninistically,
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 -ths is done, for instance, in van der Meulen et al. [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 3 of fo. the posterior contraction rate-up to possible In terms-is of order n. #))28++d; see for instance Arbel et al. [2] and Rousseau and Szabo [43].
By analogy of the Fuclidean 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.
Theorem Let I be a Riemannian Mutérn Gaussian process prior on M with smoothness parameter d/2, modified to truncate the infinite sum to at least Jn, d(min(1.0/)) Assumption 3. Then there is a C >> O such that
25
terms, and let fu satisfy
(B.)
Ex.y mmyilf FollE2(00)
CCn"
204d
(10)
The proof is essentially -thel---themmmeeee as the non-truncated Matérn, but involves tracking dependence of the inequalities on the truncation level Jn, 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 iinit-trunncation prior places its probability on C (M). Again, the constant C' >>> 0 can be seen to depend on d,o²,, V.K.B.po.o", [l (=( Il pg....) (M) and [[oolcce This concludes our results for the intrinsic Riemannian Matérn priors. We now study what happens if, instead of working with a geometrically model, we simply embed everything into Rd and formulate our models there.

---

3.3 Extrinsic Matérn Gaussian Processes
The results of the preceding sections provide good reason 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 V/ with the regularity of fo- Another possibility is to consider an extrinsic Gaussian process, that is, Gaussian process defined over an ambient space. This has been considered by Yang and Dunson [63] for instance for the square- exponential process, in adaptive setting where one does not assume that the regularity / of fo is explicitly known, but where 3 < 2. In this section we prove a non-adaptive analog of this result for the Matérn process.
Definition 7 (Extrinsic Matérn prior). Assume that the manifold M is isometrically embedded in the Euclidean space RDD, such that we can regard M as subset of R.P. Consider the Gaussian process with zero mean and kernel given by restricting onto M the standard Euclidean Matérn kernel
KurNoo3 (2, 2')
of -- T(v)
(220
1
x'l =2
)
Kv V2V
x' FRR
(II)
where o, V O and Ky 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, would expect it to be less performant when 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.
Theorem 8. Let f be mean-zero extrinsic Matérn Gaussian process prior with smoothness parame- ter > d/2 on M, and let fo satisfy Assumption 3. Then for some C G> we have
2 2.n((,.)
Eg E, f-1T( [m foll S Cn
20++d
(12)
Theorem 8 is a surprising result because the optimal rates in this setting only require the knowledge of the regularity 3, 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 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 Sobolev space H (T) DT= restricted process has non-explicit RKHS, which necessitates further analysis. We tackle this issue by using results from Große and Schneider [24] relating manifold and ambient Sobolev spaces by linear bounded trace and extension operators, and from Yang and Dunson [63] describing 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 (T++d/((M). which allows us conclude the result in the same manner as in the intrinsic case.
As consequence, our argument applics 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 kernlls defined on the sphere-see e.g. Wendland [61]. Chapter 17 and Hubbert et al. [26]. The constant C > O can be seen to depend on d,D.o?, V.R.S.po.o? [l Ha B& ..yylllo (M)) that here C' depends implicitly on 1 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 Po: for 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 (po) -norm, which depends on d, v, and 8, and is optimal if the regularities of fo and the prior match. These results imply the following immediate corollary, which follows by convexity of (1))) using Jensen's 's inequality.
7

---

Corollary 9. Under the assumptions of Theorems 5, 6 and 8, it follows that, for some C' > 0
2mind B.v)
Eny | - foll S Cn
2044
(13)
where mll(-a the posterior mean given a particular value of (si, M:)
When v = 3. the optimality of the rates we present in the manifold setting can be easily inferred by lower bounding the 1,7-risk of the posterior mean by the r'"risk over a small subset of M and using charts, which translates the problem into the Euclidean framework for which the rate is known to be optimal-see for instance Tsybakov [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 R-in in either the intrinsic, the extrinsic case case-but only by the smoothness parameter v. Indeed, the RKHS of the process is only determined- up norm equivalence- by v, which plays 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 geomeiry in our setting, as was already hinted by Yang and Dunson [63]. Hence, in the geometric setting, optimal posterior contraction rates should be thought of more basic property that any reasonable model should satisfy. Differences in performance will be down to constant factors alone-bus as we will see, these can be significant. To understand these differences, we turn to empirical analysis.
Experiments
From Theorems 5, 6 and 8, 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 9 behave numerically. Specifically, we consider the pointwise worst-case error
(T)
sup
(T) E Ee ([ (t) fo(t))
(14)
[oflll 11
where m (7) (. is the posterior mean corresponding to the zero-mean Matérn Gaussian process prior with smoothness v, length scale π, amplitude oy. which intrinsic if 7 = i or extrinsic if r = e. We use a Gaussian likelihood with noise variance oe and observations y; = fo (z,) + and examine this quantity as a function of the evaluation location E 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 wl()) may be computed without numerically solving an optimization problem. Specifically, (14) can be caleulated, in the respective intrinsic and extrinsic cases, using
() xt K
KXx
xx +
K
Xt
(15)
(+) (K{) x R
x
X
[ (i) X
K i1)
xax x oa
(16)
where, for the extrinsic case, at = (2 (K xx + and X' is set of points sampled uniformly
from the manifold M, the size of which determines approximation quality. The intrinsic expression
simply the posterior variance km]() (t. t), and its connection with worst-case error is a well-known folklore result mentioned somewhat implicitly in, for instance, Mutny and Krause [35]. The extrinsic and arises by numerically approximating a certain RKHS norm. A
expression is cery-closely-reatte derivation of both is given in Appendix F. To assess the approximation eиror of this formula, we also consider an analog of (16) but instead defined for the intrinsic model, and compare 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 our regression model.
For M we consider three settings: dumbbell-shaped manifold, a sphere, and the dragon manifold from the Stanford 3D scanning repository. In all cases, perform computations by approximating
8

---

Dumbbell
Sphere
Dragon
100
10
10
Mean
10-2
10 2
10-1
10 10-4
10-4
10
2
100
10
10
Std. dev.
10-2
10-2 2
10
10-4
10 10-4
10
0
25
50
250
500
500
1000
Data size
Intrinsic
Intrinsic (approx. formula)
Extrinsic
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 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 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.
the manifold using a mesh, and implementing the truncated Karhunen-Loève expansion with J 500 eigenpairs obtained from the mesh. We fix smoothness 3. amplitude of 1, and noise variance of = 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 200 for the dumbbell, x = 0.25 for the sphere, and λ = 0.05 for the dragon manifold. In all cases, this yielded functions that were neither close to being globally-constee experiment was repeated 10 times to assess variability. Complete experimental details are given in Appendix 2
resembled noise. Each
The length scales is 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 s case where the full dataset relatively small, and therefore a larger set of 500 points was used instead. This allows to numerically match intrinsic and extrinsic length scales to ensure reasonably-fun comparison.
Figure 3 shows the mean, and spatial standard deviation of uz (t),,, where by spatial standard deviation we mean the sample standard deviation computed with respect to locations space, rather than with respect to different randomly sampled datasets. From this, we see that on the dumbbell and dragon sanifold-whoor 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 wherer- intrinsic models the error decays a more ppatiall--unifone 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 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, significantly less affected by differences coming from embeddings.
"oode available at: HITPS://GITHUB http://///STHU....OA/AAEEEEEEEEE WTW.:C/TTTTEEENNEEEEEEEECC SWW.:OMMMETRCCAAASPTTTIIOIIC ASYMPTOTICS.

---

In total, our experiments confirm that there are manifolds on which geometric models can perform significantly better than non-geoneeee -eoomeettiiic models. This phenomenon was also noticed in Dunson et al. [17], where prior based on the eigendecomposition of 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 mon-asymptotic amalysis, to quantify theoretically.
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- Lovve 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 non-asymptotic error analysis to capture empiriiallly-oserveed differences between extrinsic and intrinsic models.
Acknowledgments
The authors are grateful to Mojmír Mutný and Prof. Andreas Krause for fruitful discussions con- cerning 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 Comell Univer- sity, 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] C. D. Aliprantis and K. C. Border. Infinite Dimensional Analysis: Hitchhiker's Guide. Springer, 2006. Cited on page 29.
[2] J. Arbel, G. Gayraud, and J. Rousseau. Bayesian Optimal Adaptive Estimation Using a Sieve Prior. Scandinavian Journal of Statistics, 2013. Cited on page 6.
[3] T. Aubin. Some nonlinear problems in Riemannian geometry. Springer, 1998. Cited on page 14.
[4] I. Azangulov, A. Smolensky, A. Terenin, and V. Borovitskiy. Stationary Kernels and Gaussian Processes on Lie Groups and their Homogeneous Spaces the compact case. arXiv:2208. 14960, 2022. Cited page 1.
[5] I. Azangulov, Smolensky, A. Terenin, and Borovitskiy. Stationary Kernels and Gaussian Processes on Lie Groups and their Homogeneous Spaces II: non non-compact symmetric spaces. arXiv:2301. 13088 2023. Cited on page 1.
[6] S. Banerjee. Modeling massive spatial datasets using a conjugate Bayesian linear modeling framework. Spatial Statistics, 2020. Cited page 3.
[7] C. Berenfeld, Rosa, and Rousseau. Estimating a density near an unknown manifold: Bayesian nonparametric approach. arXiv:2205.1571 2022. Cited on page 4.
[8] A. Bhattacharya and D. B. Dunson. Nonparametric Bayesian density estimation on manifolds with applications to planar shapes. Biometrika, 2010. Cited on page 4.
[9] V. Borovitskiy, I. Azangulov, A. Terenin, P. Mostowsky, M. P. Deisenroth, and N. Durrande. Matérn Gaussian processes graphs. In International Conference Artificial Intelligence and Statistics, 2021. Cited on page 1.
10

---

[10] V. Borovitskiy, M. R. Karimi, V. R. Somnath, and A. Krause. Isotropic Gaussian Processes on Finite Spaces of Graphs. In Artificial Intelligence and Statistics, 2023. Cited on page 1.
[11] V. Borovitskiy, A. Terenin, P. Mostowsky, and M. Deisenroth. Matérn Gaussian processes Riemannian manifolds. Advances in Neural Information Processing Systems, 2020. Cited on pages 1, 1,2, 5, 15, 23, 25.
[12] I. Castillo, G. Kerkyacharian, and D. Picard. Thomas Bayes' walk on manifolds. Probability Theory and Related Fields, 2013. Cited on pages 4, 5, 15, 32, 33.
[13] 1. Chavel. Eigenvalues in Riemannian Geometry. Elsevier, 1984, Cited on pages 2, 14.
[14] T. Coulhon, G. Kerkyacharian, and P. Petrushev. Heat Kernel Generated frames in the setting of Dirichlet Spaces. Journal of Fourier Analysis and Applications, 2012. Cited on pages 5, 15,
32.
[15]
S. Coveney, C. Corrado, C. H. Roney. D. O' O'Hare, S. Williams, M. D. O'Neill, S. A. Niederer, R. H. Clayton, J. E. Oakley, and R. D. Wilkinson. Gaussian Process Manifold Interpolation for Probabilistic Atrial Activation Maps and Uncertain Conduction Velocity. Philosophical Transactions of the Royal Society A, 2020. Cited on page
[16] E. De Vito, N. Mücke, and L. Rosasco. Reproducing kernel Hilbert spaces on manifolds: Sobolev and Diffusion spaces. Analysis and Applications, 2021. Cited pages 6, 14, 15, 17, 23, 25, 32.
[17] D. B. Dunson, H.-T. Wu, and N. Wu. Graph Based Gaussian Processes on Restricted Domains. Journal of the Royal Statistical Society Series Statistical Methodology, 2021. Cited on pages 4, 10.
[18] C. Fefferrman S. Mitter, and H. Narayanan. Testing the Manifold Hypothesis. Journal of the American Mathematical Society, 2016. Cited on page 2.
[19] S. Ghosal, J. K. Ghosh, and A. van der Vaart. Convergence Rates of Posterior Distributions. The Annals of Statistics, 2000. Cited on page 4.
[20] S. Ghosal and van der Vaart. Fundamentals of Nonparametric Bayesian Inference. Cam- bridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press, 2017.
Cited on pages 3-6, 17, 21, 28,31. 31.
[21] E. Giné and R. Nickl. Mathematical Foundations of Infinite-Dimenssioiee Statistical Models. Cambridge Series in Statistical and Probabilistic Mathematics. Cambridge University Press,
2015. Cited on pages 14, 15, 32.
[22] I. S. Gradshteyn and I. M. Ryzhik. Table of Integrals. Series, and Products. Academic Press,
7th edition, 2014. Cited on page 7.
[23] A. Grigoryan. Heat Kernel and Analysis on Manifolds. American Mathematical Society, 2009. Cited on page 14.
[24] N. GroBe and C. Schneider. Sobolev spaces Riemannian manifolds with bounded geometry: General coordinates and traces. Mathematische Nachrichten, 2013. Cited on pages 7, 14, 24,
32.
[25] E. Hebey. Sobolev spaces on Riemannian manifolds, volume 1635. Springer, 1996. Cited on page 14.
[26] S. Hubbert, E. Porcu, C. J. Oates, and M. Girolami. Sobolev Spaces, Kernels and Discrepancies over Hyperspheres. Transactions on Machine Learning Research, 2023. Cited on page 7. [27] N. Jaquier, V. Borovitskiy, A. Smolensky. A. Terenin, T. Asfour, and L. Rozo. Geometry- bayesian optimization in robotics using Riemannian Matérn kernels. In Conference on Robot Learning, 2022. Cited on page 1.
[28] M. Kanagawa, P. Hennig, D. Sejdinovic, and B. K. Sriperumbudur. Gaussian processes and kernel methods: A review on connections and equivalences. arXiv:1807.02582 2018 Cited on pages 24, 25.
[29] L. Le Cam. Convergence of Estimates Under Dimensionality Restrictions. The Annals of Statistics, 1973. Cited on page 4.
[30] D. Li, W. Tang, and S. Banerjee. Inference for Gaussian Processes with Matern Covariogram on Compact Riemannian Manifolds. Journal of Machine Learning Research, 2023. Cited on page 4.
[31] W. V. Li and W. Linde. Approximation, Metric Entropy and Small Ball Estimates for Gaussian Measures. The Annals of Probability, 1999, Cited on pages 32, 33.

---

[32] M. Lifshits. Lectures on Gaussian processes. Springer, 2012. Cited on page 17.
[33] F. Lindgren, H. Rue, and J. Lindström. An explicit link between Gaussian fields and Gaussian Markov random fields; the stochastic partial differential equation approach. Journal of the
Royal Statistical Society: Series B (Statistical Methodology), 2011. Cited on page 5.
[34] S. V. Lototsky and B. I.. Rozovsky. Stochastic partial differential equations. Springer, 2017.
Cited on page 26.
[35] M. Mutny and A. Krause. Experimental design for linear functionals in reproducing kernel hilbert spaces. Advances in Neural Information Processing Systems, 2022. Cited on page 8,
[36] Z. Naulet and É. Barat. Some Aspects of Symmetric Gamma Process Mixtures. Bayesian Analysis, 2018. Cited on page 5.
[37] M. Niu, P. Cheung, L. Lin, Z. Dai, N. D. Lawrence, and D. B. Dunson. Intrinsic Gaussian processes complex constrained domains. Journal of the Royal Statistical Society: Series B (Statistical Methodology), 2018. Cited on page 1.
[38] T. Pinder and D. Dodd. GPJax: A Gaussian Process Framework in JAX. Journal of Open Source Software, 2022. Cited on page 35.
[39] E. Porcu, M. Bevilacqua, R. Schaback, and J. Oates. The Matérn Model: A Journey through Statistics, Numerical Analysis and Machine Learning. arXiv:2303.02759 2023. Cited on page 2.
[40] A. Rahimi and B. Recht. Random Features for Large-Scale Kernel Machines. Advances in Neural Information Processing Systems, 2007. Cited on page 6.
[41] C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning. MIT Press,
2006. Cited on page 3.
[42] J. Rousseau. On the frequentist properties of Bayesian nonparametric methods. Annual Review of Statistics and its Application, 2016. Cited on pages 3, 4.
[43] J. Rousseau and B. Szabo. Asymptotic behaviour of the empirical Bayes posteriors associated
to maximum marginal likelihood estimator. The Annals of Statisties, 2017. Cited on pages 4, 6.
[44] J.-B. Salomond. Testing Un Un-Separated Hypotheses by Estimating a Distance. Bayesian Analy- sis, 2018. Cited on page 5.
[45] D. Sanz-Alonso and R. Yang. Unlabeled Data Help in Graph-Based Semi-Supervised I Learning:
A Bayesian Nonparametries Perspective. Journal of Machine Learning Research, 2022. Cited page 4.
[46] L. Schwartz. On Bayes procedures. Zeitschrift für Wahrschhinnlichkettttheen und Verwandte Gebiete, 1965. Cited on page 4.
[47] M. L. Stein. Interpolation of Spatial Data: Some Theory for Kriging. Springer, 1999. Cited on page 3.
[48] S. Strichartz, Analysis of the Laplacian on the complete Riemannian manifold. Journal of functionul analyvis, (2(1):48-79, 1983. Cited on page 14.
[49] W. Tang, L. Zhang, and S. Banerjee. On identifiability and consistency of the nugget in Gaussian spatial process models. Journal of the Royal Statistical Society Series B: Statistical Methodology, 2021. Cited on page 4.
[50] H. Triebel. Theory of Function Spaces I. Birkhäuser. 1983. Cited on pages 15, 27.
[51]
H. Triebel. Theory of Function Spaces II. Birkhäuser, 1992. Cited on pages 14, 15.
[52]
A. B. Tsybakov. Introduction to Nonparametric Estimation. Springer, 2008 Cited on pages 4,
8.
[53] F. van der Meulen, M. Schauer, and H. van Zanten. Reversible jump MCMC for nonparametric drift estimation for diffusion processes. Computational Statistics and Data Analysis, 2014. Cited on page 6.
[54] van der Vaart and H. van Zanten. Adaptive Bayesian estimation using a Gaussian random field with inverse Gamma bandwidth. The Annals of Statistics, 2009. Cited on page 4.
[55] A. van der Vaart and H. van Zanten. Bayesian inference with rescaled Gaussian process priors. Electronic Journal of Statistics, 2007. Cited on page 3.
[56] A. van der Vaart and H. van Zanten. Information Rates of Nonparametric Gaussian Process Methods. Journal of Machine Learning Research, 2011. Cited on pages 2-6, 17-22, 30, 31.
[57] van der Vaart and H. van Zanten. Rates of contraction of posterior distributions based on Gaussian process priors. The Annals of Statistics, 2008. Cited pages 4, 5, 17.
12

---

[58] A. van der Vaart and J. Wellner. Weak Convergence and Empirical Processes. Springer Series
in Statistics. Springer Verlag, 1996. Cited on page 22.
[59] van Zanten and A. van der Vaart. Reproducing kernel Hilbert spaces of Gaussian priors. Institute of Mathematical Statistics Collections, 2008. Cited pages 17, 30, 32.
[60] X. Wang and G. Lerman. Nonparametrie Bayesian regression manifolds via Brownian motion. arXiv:1507.06710, 2015. Cited on page 4.
[61] H. Wendland. Scattered Data Approximation. Cambridge University Press, 2004. Cited on pages 7, 25.
[62] P. Whittle. Stochastic Processes in Several Dimensions. Bulletin of the International Statistical Institute, 1963. Cited page 5. [63] Y. Yang and B. Dunson. Bayesian manifold regression. The Annals of Statistics, 2016. Cited
on pages 4, 7, 8, 24.
[64] W. W. Yoo, V. Rivoirard, and J. Rousscau. Adaptive supremum posterior contraction:
Wavelet spike-and-slab and anisotropic Besov spaces, arXiv: 9199,, 2017. Cited on page 6.
13

---

A Preliminaries
Here we describe preliminaries necessary for Appendices B to E. This includes some basic properties of the Bopllce--Beltrame operator 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 certain technical lemma. Hereinafter, the expression some constant C > 0 whose value is irrelevant to our claims.
means I Cb for
A.1 Laplace- Operator and Subordinate Partitions of Unity
Let M denote a compact connected Riemannian manifold without boundary of dimension d E Z The Laplace- saplace-Beltrami Beltrami operator A on M is self-adjoint and positive semi-definite [48, Theorem 2.4]. Let (L²(M), ( ((-))) denote the Hilbert space of square integrable (equivalence classes of) functions on M with respect to the standard Riemannian volume measure.
By standard theory [13, 23], there exists orthonormal basis {زز of (M) consisting of the eigenfunctions of A, namely Aʃ₃ - with X 2 We assume that the pairs (Aj. sorted such that 0 = À S X S Aj+1- The growth of X, can be characterized as follows.
Result 10 (Weyl's Law). There exists a constant C > 0 such that for all j large enough we have
c-1 >b s C/2la.
(17)
Proof. See Chavel [13], Chapter 1.
Following De Vito al. [16] and Große and Schneider [24] we fix a family T (14, ét, Xi) where Zo the local coordinates Φ: 24 C M = (24) C Rd are smooth diffeomorphisms, and the functions form partition of unity subordinate to (2t) that is xi E C (M), supp(xi) C Ui 0 S XI S 1 and Σ, Xi 1-here, can choose finite by compactness of M. For convenience and without loss of generality, assume that Vi- = (0, 1)d, With this, we can start defining function spaces on M.
A.2 Hölder Spaces CH2 and the spaces of continuous functions C and smooth functions C
Consider an arbitrary domain A' C Rd or S M. We denote the class of infinitely differentiable functions on X' by (x)). Let (X)) denote the Banach space of Z,o times continuously differentiable functions on X' with finite norm
1flees sup V f(x) ZEX
(18)
where v is the kth covariant derivative, as in Hebey [25], Section 2.1 and Aubin [3], Definition 2.2. We also write C(X) (x)) for the space of continuous functions on X.
The Euclidean Hölder spaces CH" (IRd), where = k + with E Zza and 0 < a S 1, are defined b
CH²(R = {s E CR(R)) ll [[ <00
(19)
where
lf(x)) f(y) | y R
l[fleee((9) = IIffll +
sup
(20)
B.ERR afy
More information on these definitions may be found, for instance, in Giné and Nickl [21] and Triebel [51]. We now turn to the manifold versions of the Hölder spaces.
'Hölder spaces are often also denoted by C as well, with E Ryo. Since, using this formulation, they do coincide with C* (X') when & Ezo. the notation CH" avoid confusion.
14

---

Definition 11 (Hölder spaces), For all > O we define the Hölder space CH. (M) on the manifold M be the space of all f M + R satisfying
1 Ellann o7 (ICHH 80..
(21)
Since the charts @ are smooth, Definition 11 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 Triebel [51]. For Sobolev spaces we use the poassll-otentiil-beee definition, following De Vito al. [16].
Definition 12 (Sobolev spaces), For any we define the Sobolev space H (M) on the manifold M the Hilbert space of functions f (MM such that If = (f. ff: oo where
2
H-
(22)
1=0
Remark 13. It is easy to see that substituting (1 (1+x,)² in (22) with B(a + >3)² or with (a: + 3X3) for any ax B results in the same set of functions and an equivalent norm. The former follows from Borovitskiy et al. [11], eq. (109). The latter follows from the Binomial Theorem.
For Besov spaces we follow Coulhon et al. [14] and Castillo et al. [12] and define them in terms of approximations by low-frequency functions. We fix function & E C (R2,, R2o) such that K #b) [0, 2] and ((x) for [0, 1]. We also define the functions (x) $ (22-2z)
Coulhon al. [14]. Corollary 3.6 shows that the operators Φ (v defined by
B (vA
Σ th (V ff. J20
(23)
are bounded in the space (P(M) for all 1 1sp5 S oo. Moreover, the same result also shows that express any / (MM) as f -imm--oe (v(V)) in LP(M). (VAD) can intuitively be considered as a version of f filtered by a low-pass filter. The next definition introduces the Besov spaces Be (,,()))) which are formulated in terms of ofppit--o by low-frequency functions.
Definition 14 (Besov spaces). For any & > O and 1 p.g S co we define the Besov space B (M) on the manifold M as the space of functions f E LP(M) such that |lf Il BL < 3o where
1/9
[Ifller(M) (E₁z0 2 4 (VA) ffll-s)) [11fleecm)++ + sup 21 (V A)f fll
if t00
[Ifllpp,,,
(24)
if = +00.
The classical Besov spaces B2,, coincide with the Sobolev spaces H on Fd in the sense that they define the same set of functions and equivalent norms see for instance Giné and Nickl [21] section 4.3.6-and even on manifolds if one follows the construction of Triebel [51], pages 7.3-7.4 for Besov spaces. Since our definition of the Besov spaces is somewhat non-standard, we present proof.
Proposition 15. For all 8 > 0, H* = B (MM) as sets and there exist two constants C1, C2 > 0 such that for all f (M) (M) we have
Gllifll S 11 (1fIlB I Callfflx=<MM>
(25)
'The space (MM) is the Banach space of functions (or rather their equivalence classes) that are integrable when raised to the power p co essentially bounded for p See for instance Triebel [50] for details.
15

---

Proof. It is enough to prove (25), the rest will follow automatically. The main technical tools used in the proof are Result 10 and summation by parts. First, we prove the upper bound. Denote the support set of + by K = supp(()) Notice that
(I (M)) =Ifllemm))
= >> 2230
(vĀ)s
f =1((2())
(26)
$20
Σ
2250
Σ Kfi.ILLANN
(27)
j20
by/igg2KK
$
Σ 120
22j
Σ
Kf. Kra.rr
(28)
The last inequality results from the fact that [0. C K. By Weyl's Law (Result 10), there exists
a constant C > 0 such that X S C12/d. Without loss of generality, we can assume that C' = 22r, r E Zso. Hence VA > 2 implies 1 >> 24li-r), thus we have
((Vll)))u - 11 -1 2 Σ 22j
(29)
120
(>2dG-r)
Σ 223 Σ
+ +E 22 E DinDam
2 (30)
Usisr
1>20(5-))
s2r
1>240-9
111 + 22re
Σ2
(fi. fhrmml
(31)
120
1>26
Now let R₁ = E1>zil Kfi. f) (o and s5 = =E 22je <
1 So = 0. Write
22
22Js
2--11
Σ
22js
(11 1 = E(s, - -)
(32)
1>0
1>24
130
Σ S, - RBzt1)
(33)
1>0
Σ>S Sj
Σ
(34)
1>0
26 <1<122011)
22
22js
Σ
Kf., KisDeam ?
(35)
22
5>0
20 <1
2
Ra/d (fi, (fi
(36)
2
3>0
24 23225241
22
Σ
AlKfi 2
(37)
2
320
2d
c1g246+1
22
(fl flemm 2
(38)
2
222
Ai
(39)
22s
120
where we have used Result 10 to get existence of a e such that 12/d S cd₁, This proves the upper bound, The proof for the lower bound is similar.
Proposition 15 provides characterization of the Sobolev spaces Ha HMM). There is yet another charac- terization of these spaces that will be useful later, in terms of charts. We present this characterization as part of the following result.
16

---

Theorem 16. On the Sobolev space (M), the following norms are equivalent:
1/2
I/ n- (M)) Σ (1+
(1
(40)
1/2
|I (2
[Ifll +
E(2
Σ
2 2 1 (VA) (lllaama)
(41)
120
Σ (xa)) oz "1 (1
(42)
Proof. The equivalence between (1-)) and II (MM))) is given by Proposition 15, The proof of equivalence between (=))) and ( (M4) can be found in De Vito al. [16].
A.4 Gaussian Random Elements
Here we recall the definition of a Gaussian process as a Banach-space--alle valued random variable, following for instance van Zanten and van der Vaart [59].
Definition 17 (Gaussian random element). Let (B, I1-llR) be Banach space, and J be Borel random variable with values almost surely. We say that f is Gaussian random element if (J) is univariate Gaussian random variable for every bounded linear functional b* BB.R.
Random variables of this kind are also sometimes called Gaussian in the sense of duality. One should think of 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 be an actual function opposed to, for instance, a measure or a distribution. Many connections between the usual Gaussian processes and Gaussian random elements exist, see Lifshits [32]. Ghosal and van der Vaart [20], Appendix van der Vaart and van Zanten [57] for details. The following observation about Gaussian random elements will be useful later.
Lemma 18. A Gaussian process f on the manifold M with almost surely continuous sample paths is a Gaussian random element the Banach space (C(M), ll- (ellla)))) of continuous functions on M.
Proof. Since C(M) is separable, this follows from Lemma I.6 in Ghosal and van der Vaart [20].
A.5 A Technical Lemma
In order to apply Bernstein's 's inequality when going from the error at input locations to the L (po) error, we will use the following technical extrapolation lemma.
Lemma 19. For any function g M + R, a number y E R>o Z and density : + R,o with S Po, we have
I & 11 (M) Ilg
(43)
Proof. We use Lemma 15 from van der Vaart and van Zanten [56] and push it through charts. More precisely we have, using 30, s0 (10, 1] CH CH*((O, 1) for y $ Z>o. that
11allz~ S (x11g) o क, ' (1)~()))
(44)
[max][(x])) o db
of "1
(45)
By definition the the manifold Hölder spaces, this gives
|Illl (M)) S ( max||(X19) 07 "1111220))
(46)
17

---

Finally, since XI are bounded, the charts are smooth, and is lower bounded, we have
((xig) of t 12:000 -
(rrto))//,,, ((x0s) o (w) dy S / (())))(())((d)) & (lalliz
(47)
which gives the result.
This concludes the necessary preliminaries. We now turn to the proofs of our main results.
B
Proofs
Throughout this section the notation from Appendix A and results from Appendices C to E. We start by defining our notation for Gaussian likelihood and probability distribution of the sample.
Definition 20. For every x M and f -M R. we define Paly to be the joint distribution (f) = N(J(a), where [xx) is the vector with entries (xx Expectations with respect to Pr.y. pa and p (f) we denote by B by Eg and by [E)) "/w² respectively. Sometimes, when f clear from the context, we omit the subscript (f),
corresponding to the marginal pa po and conditional Pyle
(f).
Using van der Vaart and van Zanten [56]. Theorem 1, which is valid for any compact metric space and thus also for M, we can deduce a posterior contraction rate at data input locactions, with respect to the empiricul L2"-norm
1/2
1 Σ f(z))2 1=1
Ifll,
(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 II, may depend on if we consider truncated intrinsic Matérn process.
Theorem 21. Let II, denote the prior either Theorem Theorem 6 or Theorem with smooth-
ness parameter >> d/2. Let HI, denote the corresponding RKHS. Define the CONCENTRATION FUNCTION for fo E C(M) and > 0 by
inf (5) = - In [[lll==(M) < <e] <0]+ JEH.: I15 i
(49)
<e
Then if fo E HP(M) n B& (M), > 0 we have ( S ne, ner en multiple o
min(u. 2r12
Proof. The first term the right-hand side of Equation (49) is bounded by C'-d/v by Lemma 33. To bound the second term, we assume, without loss of generality," that v 2 3. Consider an approximation f (v(A)) fo fo of fu. where es < 2-j < E and c does not depend on Since we assume fo E 00. (M), by definition of B ox we have
Il flle=(() S [l B (MM) 35 Se
(50)
where the last inequality the o (M)-norm is the constant implied by notation S. We now show that
| -5-A(u--+ap7
(51)
First notice that by Lemma 24 and Proposition 27, for any prior considered here have H. S H=+4/2 (M) and I-I1,, < |1-Il with constant that does not depend Hence using
*This actually seminorm, but we follow rest of the literature referring it as norm. IIf @ > then fo H HM) osc C (MM) oo (M) gives en 2u-4
min(8. 20t4
18

---

Result 10 (Weyl's Law) and properties of Q, we have
|| I
2 || (M2
(52)
= -1 / (0000/2222211// (2 (fi, 50) 1
(53)
120
$ Σ
(1+ (1+1) (1 + 21 (fr. fo) Lo(M))
(54)
1V11522t1
< 2/
Σ (1+ (fr. fo) 1
(55)
lVAig21+1
S <20t2)o4d/2-51 8)
(1+
(56)
20+2) (M)
(57)
S Se-A (15o (M S $0
(58)
Our assumption 2 8 implies that
2 n.
/ 2 B ".
d d
(59)
Hence, we have c d/v d- = 3 (v--8( 84d/2) which gives us & (E) 2 c - is then easy to check that e, = Mn 2o1d satisfies 1 (Ex) I nE for M > 0 large enough.
From this, we deduce an upper bound on the error in the empirical L² norm II-I that is, on the Euclidean distance between the posterior Gaussian process f and the ground truth function fo evaluated at data locations xi.
Lemma 22. Let 1I, denote the prior in either Theorem 5, Theorem 6 or Theorem 8 with smoothness parameter v > d/2. Fix fo E HMM) n B o0 (M) with 3 > 0. Then
EE lol.. S e
(60)
minier.9)
for En O Ti
with constant depending on So, but not on x.
mint0,.)
Proof. By Theorem 21 for e, a multiple of n
we have 4. (E,)) S. ne?, By virtue of this, the proof of Theorem and Proposition 11 of van der Vaart and van Zanten [56] imply the result.
20-td
Indeed, the proof of Theorem 1 relies solely on the fact that (en/ S ne? and application of
van der Vaart and van Zanten [56], Proposition 11. We have 9 (en) S ne S m((an) ) and hence the condition is satisfied with replaced by 2e..
D
We now turn to the proofs of our main results, Theorems 5, and 8, which for convenience we restate below, The idea of these proofs is to extend the result of Lemma 22 from input locations to the whole manifold A using an appropriate concentration inequality. To this end, the proof closely follows the one of Theorem 2 in van der Vaart and van Zanten [56], but relies on Lemma 22 proved above and on the concentration incquality we prove in Appendix D along with some important sample differentiablity properties of the prior processes.
Theorem 5. Let f be a Riemannian Matérn Gaussian process prior of Definition 4 with smoothness parameter d/2 and let fo satisfy Assumption 3. Then there is C > 0 such that
2 min[B..)
Ev (e folli Cn
20+d
Theorem Let be Riemannian Matérn Gaussian process prior on M with smoothness parameter > d/2, modified to truncate the infinite sum to at least Jn ((min(1100)) Assumption 3. Then there is a C > O such that
20-4 terms, and let fo satisfy
2 2min
Ex,y E EX-I((- wilf follEAmm)
C'n
20+d
(10)
19

---

Theorem 8. Let f be mean-zero extrinsic Matérn Gaussian process prior with smoothness parame- ter d/2 on M, and let fo satisfy Assumption Then for some C >>> we have
2 (B,n)
Ex,y E (. f Cn
20+d
(12)
Proof of Theorems 5, 6 and 8. To ease notation, we denote the expectations under the true data generating process by Fey Fo(oo) and Fyya = TPloo) min(3 2442
omitting the superscript (-)% Take
"y/="
En satisfying (5n/2) I ne?, noting that such rate exists by Theorem 21. Then, for each n, there exists an element fn E Hi,, where H,, the RKHS corresponding to II, satisfying
Il li. S ne
Il ( S En/2.
(61)
Write
on Buy R,
foll2 (aw
S en Em,y R. (n Sol (p0))=()))
(62)
+ E Exr E E (1 Ils - [(2(p0)
(63)
$ 1 - en Euy R II, - fm (22))
(64)
Thus, we can work with f - fn instead of f - fo- Define B(r) = {llf fall > Enr}. Then
en (I II S- (2: = [0 2r P fr (B(r)) dr.
(65)
Fix such that d/2 < 7 v $ Z and 0, and define
(r) = {2 - fulln > Enr}
(66)
B(l (r) = {IlIllewsem) > Vnenr=
(67)
B (r) -B(r) (BW)() U (
(68)
Then B(r) S BdD) (r) (r) (())BBBU (Bll)) and thus for an indexed family of events Ap to be chosen later, have
en By By-n.(Hew)) -n 2 ..00 r r (r) dr Lo/ r dr
Iullia % rlA, PP P (B dr+/ + 1/ r1A, PM.. (BII)()) dr. (70)
(69)
For the first term, by Lemma 22 applied conditionally on æ, for which we got a bound on the integrated empirical L -o-m uniformly on the design points, have
Ex,yy 10 (P (r) dr S = Ex.y (1mm)) foll? < 2 = (71)
Moreover, by Lemma 14 of van der Vaart and van Zanten [56] applied with the notation of the reference being equal to Vnewr" ", for each >> O, the event
Pyl
(
(w)
n PMNI, (IIf fff-follmmmn (a <enre)
Ax)
R
dII,(f)
(72)
(1))) (w)
Pyl (w)
alIIn() > e 2 pT-I,, (1IS- <enr%)
5
E Rh,
(73)
(T)
Pylæ (u)
is such that
P Pg/ar [x)(z)) S e n 20/8
(fo)
(74)
20

---

It should be noted that the v factor disappears because of the discrepancy between the empirical L norm I1-1l, and the Euclidean norm used van der Vaart and van Zanten [56]. By Fubini's Theorem
and since ne 2042 2 1, the second term is bounded by
Ex. ! (l) dr
r EE Bul [LAp(m)]] dr
(75)
- dr
(76)
0
re 20 18 dr C. < co.
(77)
It remains to bound the last two terms. By Bayes' Rule, we have the equality
StVne, Paja (3) (f),
Py- (Ilfllenr(s) > T (nenn")
(78)
S S (f)
P (w)
Juflleeee STvee
Pylm Ln (l),
(79)
Paie (v) dIT,(())
(7)
S
P (w)
therefore for (x)), have
IP ( (IIfllcH TVEN))))
(80)
I Pywdn (Ilf foll < EnTe) /IllIATTIMN> 27 veen
en
Pyla (2)
(1)
dlT,(().
(81)
P Pyj (v)
(J0)
(f2()))
Hence taking expectations, using Tonelli's Theorem, and Emsy
gives
Pala (w)
Ex,y [AA(m) PSNTa(()a (irfeei(IIIlleeeeee (11flleneres) > TVnent")
(82)
2
Ev, J
y (v)
dlI,,()
(83)
Prwel (115 Sollz=(Mn) < <Enre)
py (z)
Pywtt, (115 (o < <enre) J IIfilenrunnirveee.... Vnenn
p (f)
(
Ex
dII,(f)
(84)
Py (3)
e
PyweI (I1Illenrcs) > =Vven)))
(85)
PyLII (1f--oolza < <enre)
Therefore, the third term can be bounded
(B() (r) dr
Ex. r1A. Powtia(-
(86)
Jn r PPIII (115 foll ( <enre) PA-m, (IIAleww > TVTIS dr.
(87)
Now, using Lemma 31, for possibly small constant independent of have
Fywm.((1..)) (IIfIlcCET()) > TVTEn)" S cT² ?
(88)
Moreover, by using the bound on the concentration function Theorem and Ghosal and van der Vaart [20], Proposition 11.19, assume that
P P [|I- - -foll=<(u) < Enri 2 ze-c -e-e-ine Ineg
21

---

Therefore, the third term is bounded by
Ev,y / rlA.()) PyTT, (r) (B dr - jo re me- e nner dr (90)
0
re- dr 00
91))
if r2c > 2+ It remains to bound the last term. We have, by the same arguments as above, that
Eay 6 (l) Py (ynln(-e))))) (B(I)( (r)) dr
(92)
-Exy fo (1))
(93)
P (Ifllexr(m)) S &TVTEN"",11 full, I enr S Iff - fullep(po)) dr (94)
enegre
1 Pyne, (IS Sollew(M) < <Ear))
(95)
Ex pyWt (IIlleezeemm S STVTENt" 2 fulla SET 1if-sslfflllll fflll dr
(96)
relot relot1)ne2r
(97)
x
Ex 11f (2)1) full, dII, (f) dr.
(98)
ST S.
As the squared empirical -orm is Monte Carlo approximation of the true L the probability in the integrand can be controlled via concentration inequality. As in van der Vaart and van Zanten
[56], we use Bernstein's 's inequality [58, Lemma 2.2.9]. For collection Yi.... Yn of random variables
such that Y. 0 and Y; E [- M, M] almost surely for some constant M > 0, this inequality asserts
x²
P(Y + ................................ + Y.l x) S 2 exp 2 Mx/3
(99)
where v 2 Var(Y1 +... + Yn). We put Y, = ( - (x.))) IID with H; Po. It is easy check that IE Y, and Y; [-M, for M = =lll full (M)) Furthermore, llf fn -f upper-bounds the variance of the respective sum, since
are 111 - full22(p))) where x,
1
Var(Y; + + Yn) Var Varz=po (f(xx) fn(x))² = n
Ez-ro (f(x) ((=))_f(z))) fn(z))*
(100)
|| fulli- Eaapo (f(x) [n(x))² = u.
(101)
n Using Bernstein's inequality with Y, M, v as above and with x 211 - S'll 2. ((()))> we have
Exz 11ss- 2211 Sall,
F. IIs fn (211)11,,,010111100 24115-
(102)
=E, 1 1- 11 -111--u-1 - (m)))
(103)
x²
S exp
(104)
2 + Mx/3
9
|l
fn ((pp)))))))
exp
32
(105)
μ+ +1ll - fn 11. (M)
2|l
Sulli
On 4 Ilf fm IIEAPEO))
exp
(106)
32 5 || 12122- (M)
9n ||/ - fn -m.112A(m))
exp
(107)
40 [| 1
22

---

Moreover, by Lemma 19, since & Z>u. we have
V/ - La (allam( - - In (M |15 fr (p()))
(108)
Using the Sobolev Embedding Theorem, given in De Vito et al. [16], Theorem 4, ll. fn fullerm( & 1l + (M)))))))))) & TVVEEr" whenever [[flCcт(.) I TyneEr" Therefore, for a constant > O, when ( S Tvvne and enr SII Jnll Jnll we have that
1l fu Fnll12(00)
Eg P1 -00 exp
(109)
47
|I fm (m)) - fn 12(
2d 27
24 (1-2)
S -
(110)
Hence, we can bound the last term as
(B (r)) dr
Exsy
rIA P P)-ta((MM))
(111)
re (ot1)ne -rr 25 - 2 dr.
24 27
00-0
S 0
(112)
We have n 3. n 1542
Since En S min(r.9)
and min(v, ()) > d/2, we have ne 2
n 2--a for some 7 E (d/2, v). Moreover, for this choice of y and s small enough we have 2y+d (1
2d
s) 2 2s, which proves that for some constants C, C, C" > the fourth term bounded by
C cf re -0, dr
(113)
This concludes the proof.
Remark 23, Following the proof, easy to see that [If foll on the left-hand side of Equations (9), (10) and (12) can be replaced with ||5 - foll for any >> 1, changing the exponent 2in((5.) the right-hand side q qmin(8,v) 20-d
2v-d
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 Borovitskiy et al. [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 24. Let HH be the RKHS the the intrinsic Matérn Gaussian process with smoothness parameter truncated at the level J E Z,z0 U {oo), and let {fs) be the orthonormal basis of Laplace- Beltrami eigenfunctions. The space Hi, is norm-equivalent-vvie with constants depending only on v, and o to the set of functions f = E bySz with by E TR, equipped with the inner product
Sa E6,))
(1+ 2 22/2
(114)
3=1
In particular, Hl, H=+d/2 for all J, for every h E Hiy we have ||AlI
[ a
Proof. By direct computation, the covariance k of the (truncated) intrinsic Gaussian process is
of
2v
+>)
-(v+d/2) fxx)
k(xxx')
(115)
Cu.
R2
j=1
23

---

Hence, the covariance operator K (M) - (M) defined by
(Kf)(x)= / k(xxz')
(116)
is diagonal in the basis Us/ with Kʃ [28]. Theorem 4.2 implies that H. consists of functions of form Σj= azsz satisfying
C..e ( 2g) (u+d/2) fz. Then, Kanagawa al.
v+d/2
2v
21
lay| < <00.
(117)
Cu,x ==1
Using the elementary inequality min k2 1) norm-equivalent to the space H²/4/22
1+1 I max 2v 1), we find that this space is (M) of functions f Eラ=ı k2
S
agf, satisfying
| my
(M)
Σ
(1
2,) <
(118)
j=1
where the comparison constants parameters K, oj.
min
(1, 2
and v
max 24) only depend the Cz.a
of
Cor
The next two theorems will be useful to characterize the RKHS of the extrinsic Matérn process on M. We start by a lemma relating the RKHS of the restriction of Gaussian process to the original one.
Lemma 25. Assume that l is a kernel on R4, f ~ GP k) with almost surely continuous sample paths and THI is the RKHS of k. If M C TR is submanifold, then the RKHS TH corresponding to the restricted process fiM is the set of all restrictions 9IM of functions g E H. equipped with the norm
||al| = inf Ig
(119)
gefi, g/Mem
Moreover there always exists an element E IIIl such that 9IM f and ||ll = |1/
Proof. Lemma 5.1 in Yang and Dunson [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 26. If 8
D-d then the restriction operator extends to a bounded linear map Tr. 2 H*(R) + H=- D (M). Moreover, for every 0 there exists bounded right inverse Ex HM(M) + H D-d (R.P) such that Tr Ex idH=(M) where Tr corresponds to u u+ D-d
Proof. Theorem 4.10 in GroBe and Schneider [24].
The last two results allow us to characterize the RKHS of the extrinsic Matérn process on M. Proposition 27. The RKHS H. of a restricted extrinsic Matérn process f with smoothness parameter on M is corm-equivalent to the Sobolev space 4/2
Proof. Using Lemma 25, the RKHS Hi can be characterized the set of functions f M - that are the restrictions of some g E HI, where HI is the RKHS of the ambient Matérn process f. with
| = gefi, fIM==
inf
(120)
24

---

Since H. is norm norm-equivalee equivalent to the Sobolev space' H (R.) see for instance Kanagawa et al. [28] by Theorem 26, for every f HH we have
|| x [[Ex(f]] DD+ S Il 1r1 2 (M) = 1
(121)
Similarly, for every & E H. with gIM - we have
I| a = ||
S 111111ee 10 S |1ll|
(122)
Hence, taking the infimum, we obtain
II g(tes((M X inf
|| = ll
(123)
gEA. 9IM=
The claim follows.
D Concentration and Sample Path Regularity
In this section, prove that intrinsic, truncated intrinsic and extrinsic Matérn processes sub- Gaussian, uniformly with respect to the truncation parameter in the case 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 of independent interest. A simple way to do this would be build the results of Appendix C, using Driscoll's Theorem- given in Kanagawa et al. [28], Theorem 4.9- and the Sobolev Embedding Theorem- De Vito et al. [16], Theorem 4 but that would only give that the sample paths are almost surely in (M) for every 0 <> <<<< < d/2, E N. whereas here we improve the range of index y V.
Kolmogorov's continuity criterion a standard tool to show that given stochastic process has a Hölder continuous version: 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 28 (Kolmogorov's continuity criterion). Ifg ~ 11 is a zero-mean Gaussian process on [0, 1]d with
E [9(x)) -(v) - C'I2= y
(124)
for some O < p S 1 and C > O, then there exists version of g with samples paths in CH (l0, 1)()) for every <, a < p. Moreover for every a < p this version satisfies E. [IgllccH" S C where C" < too depends only on C, and d.
Proof. Take 2, E [0, 1] 1 M.M > and q E N. Since the random variable g(x) - -y()) is Gaussian
we have
E (g(x)) ((y) 2q (2q)!
[E [[)(()) (((((c))gvv)) -(w)))" < Cg y v1
(125)
299!
where we have defined C,2 = C4 (2g) The reason for considering the the 2yth power will become clear later in the proof. By Markov's inequality, for every 2, y E [0, 1]d we have
2
P((((z) |(6)و > u) S Cqp 2012
yll22 2gp
(126)
Now, take = Ll200 X2 = 2- - n [0, 1]dd Then, the previous inequality applied to any adjacent 2, y E X4. where see Xx graph where two vertices are connected if they differ by 2-k and u = M2-ka, ka, implics
P((g(=) m(v)l > M2 k2) < C.
202 ())
(127)
"Note that this norm- euuieellivleeeee the only property of the Gaussian process we use proofs. Any other Gaussian process satisfying this, including ones different from the Matém processes of Borovitskiy et al. [11], would also work. This is of potential interest since other Euclidean kernels, such as Wendland kemels [61]. are known possess RKHSs which norm-equivalent to those of the Matérn kernel.
25

---

[-(=)~9(v]] Summing over à: 2 and adjacent points in X-and there Denote K = sup z.ye Xadjacent 2--91
are at most C2k(d+1) of them where C > is an absolute constant- gives us for > 2(p-a)
d+1
that
P(K > M) $ Σ
Σ
P((g(z) - - > M2 ka)
(128)
k20 x. Xadjacent
CCg g
C Σ
20(d+1)
CyM-292 2ko(o-a) 20M
1 229(222pp -d- M 2g
(129)
k20
In partieular, for all 9 > max (1,
d+1 )
2(p=a)))
we have
MP(K> M) M)dM=1 =2 MMP(K (IP(K>>)dMM M) E(K3) = =2
MP(k> > M) aM (130)
2 /
MP(K> M)dM Ccapp
(131)
for some constant Cc < +00. This means that K is finite almost surely. Since X dense in [0, 1]d and g is almost surely uniformly continuous X by the classical version of the Kolmogorov's continuity criterion- see for instance Lototsky and Rozovsky [34]-g admits a unique continuous extension to (0, 1] on probability one event A. Let us define
lim g(v) on A, J-
g(x)
y
(132)
0
otherwise,
for x E [0, 1jd. For any x, y E [0, 1] and In - v, + tn, Yo E X we have on A
[g(z) - g(v)l - lim (Ig(() g(En)) + (g(xn)) - F(u)) + (ovn) ((8)))
(133)
n-000
< liminf ((g(x)) g(an)l + K K||
yn Il (9(2n) s(vil)
(134)
n-000
= =K| vllga-
(135)
On the complement of A, we have g = O. Hence, is a- a-Hölder continuous on [0, 1]a with the same (random) constant K, since
g(x)) g(v)l
E BlooleeH ((0.)4)
E sup
E TEK2 & C...ppd 60.
(136)
<x,ex | vl""
Since ÿ is version of g, the claim follows.
The next lemma applies our version of Kolmogorov's criterion, Lemma 28, to the intrinsic Matérn processes on M by considering charts. This allows us show that the sample paths are almost surely (H) for every 7 < v, which is both used our arguments and also of independent interest. For the claims in Appendix B, 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 '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
Lemma 29. Let f ~ II, be an intrinsic Matérn process with smoothness parameter v > truncated at Jn NU (oo). Then for every 7 < have
sup Es [[I < <c.
(137)
Proof. We start with the case v < 1. Take 1 < 1 5 L and define hy = (xif) Φ- Then by is Gaussian process with covariance kernel K, given by
kı(z.)) (x1 (z) (z), (v)) (2)
(138)
26

---

where 2, y E Vi and (2 y) (ov(f(z), f())) is the covariance kernel of f. Let Hi be the RKHS
induced by Kt- We seek to apply Lemma 28 h4- For all 2, y V1, where recall that assume that Vi = (0, 1)",, we have
E Ey-u, (x)()) hilv)) k)(x,)) Ri(y.v) - 2K,
(139)
|[(((,,)) - ky((, 12
-Ki(,,,)KK((v,))) >1e.
(140)
2
sup
[i( - Kily. -), -).4))
(141)
I|*lla, =1
sup 1x(xx) - (((1))2
(142)
I|+lla, =1
sup I| (v)) || v1lgdd 20
(143)
I|ela,
where (1)))czee can potentially be infinite or unbounded: we will show this is not the case. To do this it suffices show that we have continuous embedding H. - CH"(V), that is Illl x
Illlis The RKHS H4 is by definition the completion of
Σ a Ki(,, : 1, L.R.ER,,,EE
(144)
(x o7 (2:) o7 ()) (a)))) (T) :p> 1, ER (145)
with respect to the topology induced by the RKHS norm
Σ Ca,K(az
(x
(zx), o7 (z,)). (146)
Denote the RKHS of K by H. By Theorem 16, and by the equality ||-1|x = I1-1 on H which follows by Lemma 24, we have
2
Σ
0: K((x,)))
(147)
(B4d/2(Bd)
o7 (.))
(148)
H H1 d/2
(149)
Hutd/a(M)
Σ o (z:)) (x,),.
(150)
aga (xI (xx ((
(151)
a: Ki(ss,)
(152)
Therefore, we have continuous embedding HH + H²+d/2 (Ted) with |1:1| (4/ |1:11
on Hl. By the Sobolev Embedding Theorem in d see for instance Triebel [50], Section 2.7.1.
27

---

Remark 2-we have Bg,2 2g (RR) H - CH" (R²), which implies H, CH' (RR)
by composition. Thus, there exists a constant = Cy such that
Ey-u, Ih((2) - h2(1)) 2 C/|22 ylla.
(153)
for I, y Vt- Hence, by applying Lemma 28, there exists a version h of h4 with almost surely a-Hölder continuous sample paths for every a < u. Now consider -i=1 h or. Then f is version of f. We proceed to bound Ey-nn [1 For any S 1, L write
(0r((07) - hr (0,(0) (())) & (0) (o))))) - ( (x)) < CKK y1" (154)
where is the Lipshitz constant of o, o which is well defined and finite because this composition diffeomorphism and K is random constant with Ey-m1 K2 Ca.ud Hence
Ey Ey~mn, (.f..llfleeemmm = EsW,,
Σ
(xi)) o o
(155)
(HH)(R)))
1=1
= Eslli
Σ
Xi Σ
hr Dr
o oi
(156)
1m1
r=1
CRd)
Eystl,
o o ,7
(H S EyMI, K² Cav,d (157)
which gives the v < case.
We now turn to the general case. The proof will be similar to the one of Ghosal and van der Vaart [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, < 2, otherwise suffices introduce coordinates and to proceed by induction on [v]. Let E [I,...,4] and as before define Ki(x, y) = (xi o, (=) o o, (x))))))))))))) (w)) the covariance kernel of h (xi)) Φ5 well as H. its RKHS.
First, let construct Le(())deerivatve hy of hy, where (2( is the space of random variables with finite variance with (a, b) 1,(0)) E(ab). This derivative is square integrable process on V/ such that
h1(x + 6) - h((x)
E f~
hn((r) - 0
(158)
0Ki as + 0, for all x E Vj. For this, we will first show that ak (xx E Hy, for every x E Vi-here äz denotes the derivative of the function Ki(: with respect to the first argument- and that
the
0K, Ox (xx,))
0k: (x',
< C.læ 1ale--(--- 2'1011
(159)
0x
28

---

We first show that R1(x+8,)- (0(2++,)-R(((..) is a Cauchy net in IH,, We have
h Ri(x + 0, -K((x)))) Kí(x + 8), (z,
(160)
h
Ki(E + +.,)) -R((a,))
Ri(x + 8', Ki
= sup 1| =1
(161)
h
h'
B1
(0x + 8) - -p(x))
l0x + 8') sc(xr)
sup ll<e llly
(162)
h
h'
[ (s0 + th) 4 + 1h
= sup Ilella,
(163)
< sup Il li,
140'lcx=-1022) 1611-11
(164)
< sup 11 I1eele 1h --1
(165)
where in (163) the derivative of exists because, exactly as in the case S 1, can show show that H. + CH" (R2). This also implies that for constant C Cv have
Ki + 5, Ki(x,.) Kí(x + 8',) - (o)--Kip (x,))))
I Cch--h1
(166)
h'
H.
As Ih W/1-1 when h, h² + because 1, this proves that R1(x+6.)-Rin Ri(oe)) is a Cauchy net in Hr: by completeness of H, it converges in HH, to limit. Since by general properties of RKHSs, convergence in H2 implies pointwise convergence, the limit satisfies
Ki(x + d,y) Ri(a,w) ak, (x. v).
lim h-+0
(167)
h
Or
Hence the partial derivative
Hr. Moreover, by the isometry (2)) e hx(2) + Ey-HH hx(x))(() = Rı(x, E Iill, we deduce that by is actually T.2 (i))eeiffereettabble, with an d:eedd-dervva denoted as hú, and that the derivative process by is Gaussian, as it is an (2()))-limit of Gaussian random variables, satisfying Eyw11, hn(=))u())
0kL (z, v) exists for all and ak ak (z.-)
8e
Gx
a8 Ox
(x, )), aK, (y,
Having established the existence of an I,2 (S))deerivaetve hy of the process hy, we would like now to show that hi possesses (7- 1) reeular version for every y < v. For this, we would like to apply Lemma 28 to hr. Notice that, by isometry, for all h > 0 we have
Ey-T, [h((x)) ha(e) 2
0k,
0k:
(y.
(x, :)
(168)
dr
dr
i
Ki(y + 8. - Ki(y.-)
oK,
3
Or (v)
(169)
h
Ki 5, - Ki(x..)
0k,,
+
(st,
(170)
Or
11.
Ki 8, K((x.-)
Kily + 5.-) Kr(y.))
+
(171)
SSee for instance Aliprantis and Border [1] for review of Cauchy nets.
29

---

Therefore, by the same arguments as above, we have
1/2
Rí(x + -,(0,,
Ri + 6,)) Ke(W,.)
lim inf
(172)
h
h
h
( 150' + th) (0 + th]lde
< liminf
sup
(173)
h-50
Ieell, =1
S Iimmin inf Culæ - y|² = Col or - yl"
(174)
h=0
where the transition from the second-to-last line to the last line is similar to (163)-(165).
Therefore, we can apply Lemma 28 to hy, and find a version 1/g of iu with sample paths in CH almost surely for all u, such that
E f~ | '(M) Cua < to0
(175)
This gives Hölder regularity of the respective derivatives: we now integrate these to obtain Hölder- regular version of the process itself. Take any q E (0, 1) and consider hy = hε(ci) + Sa Mt(t) dt. Then since in is almost surely in CV)))(VV), hy is has sample paths almost surely in CH (V:). Moreover, it is easy to check using our previous results that hy has an dderivvviv given by h/. This implies that hy is version of hy
To conclude the argument, we construct f from the obtained parts, by pulling hy back from the charts to the manifold. Consider f Σ/= hu di. Arguing as in the case v I 1, we see that version of f with CHM(M) sample paths for every a < v. and for every < we have Ey-n, I (22) S Ca, < +00. This gives the claim.
With this, it is easy to prove that all Matérn Gaussian processes considered in this paper can be as Gaussian random elements in the Banach space (C(M), Il-llet,)) of continuous functions on M. This allows us to use the same proof scheme as in van der Vaart and van Zanten [56] through the control of the concentration functions defined in Appendix B. It is also important that work with Gaussian random elements in C(M)-and not only with the classical notion of Gaussian process, as the concentration functions are defined using the Gaussian random element RKHS defined in van Zanten and van der Vaart [59], which can potentially be different from the classical RKHS. Fortunately, when the process is a Gaussian random element in C(M), van Zanten and van der Vaart [59], Theorem 2.1 implies that the two notions of RKHS coincide.
Corollary 30. The intrinsic Matérn Gaussian processes of Definition 4, their truncated versions as in Theorem 6 as well as the extrinsic Matérn Gaussian processes of Desinition 7 are Gaussian random elements in (lll)))))))
Proof. By Lemma 18, suffices to show that the processes have almost surely continuous sample paths. Since Fuclidean Matérn Gaussian processes have continuous sample paths, this implies the same for their restrictions, the extrinsic Matérn Gaussian processes on M. For the intrinsic Matérn process, we apply Lemma 29.
Using Lemma 29 and known properties of Euclidean Matérn processes, now show, in sense, that all of the Matérn processes presented in this paper are sub-Gaussian, in a manner which holds 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. This result used to control Hölder norms when going from the error at input locations to the L (Po))-eero emphasize that the prior depends on n when we consider a truncated intrinsic Matérn process.
We use the notation II, to
Lemma 31. For f ~ II, with II, the prior in either Definition 4, Theorem or Definition 7, for every >> 0 and r < v, & Zzo. there exists constant a (f)) =0,(f) independent of n we have
for x > 0 that
P P(Il/lloto > >x +1)(()))) S 2e -222
(176)
30

---

Proof. We start with the restriction f of an extrinsic Matérn process f to M. as in Definition 7. By van der Vaart and van Zanten [56], Section 3.1, for every < v we have f E CHO ((0, 1)2) almost surely. By Ghosal and van der Vaart [20]. Lemma 1.7, for every Y v, f is Gaussian random element in the Banach space CH2 (10, 1)P). In particular, by the Borell-TIS inequality [20, Proposition I.8] we have for x > that
P(IiilCc > +x+ 1)o((f) & 2e 12
(177)
60.. Since M is smooth, the restriction f also satisfies for where a (f) = TE [|]|//lc
1/2
x > 0 the expression
p
(IIflleee())) > (I 1)o(()) (1)())(f)) s 2e-
(178)
perhaps for a possibly larger constant (f). Finally, the case of the intrinsic Matérn process f~ Il, truncated at Jn E 250 U {xo} follows in the same manner, as we have shown in Lemma 29 that supa21 Ey-n, Ey ll
(M) I Ca.
E
Small Ball Asymptotics
Here, bound the probability that a Matérn Gaussian process lies in the s-ball with respect the -nor for small e. For a Banach space B, an element æ E B and number E R>0, let us denote the closed r-ball around by B (z)). We start by an upper bound on the metric entropy of Soboley balls on M with respect to the uniform norm.
Lemma 32 (Entropy of Sobolev balls). For any > 0 define the E-covering number A C (M(M) with respect to the norm (M)) by
(e,,A,lllea arg min JEZ50
Th1 Thi,-, A:AC UBL (M)
(179)
j=1
Then for any > d/2, there exist C, > 0 such that for every Eo
In N ((,(,, (M(((M) (1)1 CC
(180)
where the left-hand side the inequality above, as a function of is called the METRIC ENTROPY of the Soboley ball B R (M) with respect to the uniform norm [I-])
Proof. Using charts we will reduce the problem to the entropy of the unit ball of the Sobolev space Ha ([0, 1]d) for which the upper bound known. Take f E B (MM) (0) and consider approximations of f by of the form
g xa o di)
(181)
for some functions gi Vi IR where Vi Rdd We have
|| - g 0012~(M0)
Σ
xila o 01 fll
S
Σ 1xlg o or - fllzz
(182)
1
Ig &i /11
Σ
|1 f a -5 "1
(183)
-=1
1=1
L max || f 1
(184)
15KLL
This means that approximate f by g uniformly M, need to choose the functions gr that approximate fo o,- well with respect the uniform norm [0, 1]d.
31

---

Next, we show that the functions f of are contained in an Euclidean Sobolev ball of radius R. with R depending only on v and the atlas, To do this we use Große and Schneider [24], Lemma 2.1 to write
7 Σ ox
| o 07 (-
=
(..(101)4))
Σ Pm1
(xr o 07 (HH((0.11))
(185)
(xv f) o Φr or o 1,
(186)
(=((0.11)))
P=1
(xe f) 0. (=((0,11") < Ilf H
(187)
Note, importantly, that the remark just above Große and Schneider [24], Lemma 2.1 allows us to consider Besov spaces B2, coinciding with the Sobolev spaces Ha instead of the Besov spaces Bixx- get from the second line to the third. Note also that the constant hidden behind the notation S in the last line where we use Theorem 16 the radius R. Without loss of generality we assume R 1. By the Euclidean counterpart of the result we are proving [21, Theorem 4.3.36), we have
In N c, B [0.1)* (0), (0110-1011-001111)1 SE-
(188)
Let h,,,,,,. hy E (101 (0) be such that B (=-"(0.11))) (0) c U BeT
(1.,1111))
(hx). Then for any
8"" n" (M) there exists sequence tirlin S [} such that
||s Σ
xhhh, o villemmm < _
(189)
1
This shows that N (e, B (M proving the claim.
(0), (H)) S LJ, where L is just the number of charts, thereby
For the related diffusion spaces [16], the RKHS corresponding to the heat (diffusion) kernels, Castillo et al. [12] uses the results of Coulhon et al. [14] bound the entropy in terms of wavelet frame instead of relying on charts. We believe this alternative proof scheme should work in our case as well. However, could not, to the best of our effort, get tight-enough bound for the Sobolev spaces by directly using the results of Coulhon et al. [14] and therefore we chose to rely on charts instead.
Having established regularity properties for our prior processes, we now turn to the small ball problem: we want to find sharp lower bounds on P( (l < <) where f ~ will be crucial in order to control the concentration functions used in Appendix B. In fact, it is well-known that this problem is closely related to the estimation of the metric entropy of the unit ball of the RKHS of ff with respect to the uniform norm: see Li and Linde [31] for details. Since we have already characterized the RKHS of our processes in Proposition 27 and Lemma 24, we are able to lower bound the small-ball probabilities. The technicality here involves getting a bound uniform in the truncation parameter for the truncated intrinsic Matérn process, as the truncated Matérn process is a sequence of priors rather than a fixed prior.
our prior process. This
Lemma 33. If II, where IL, is the prior in either Definition 4 and Theorem 6 or Definition 7 with smoothness parameter v > d/2, then there exist two constants C, > O that do not depend on n such that for all E S Eg we have In P (111( <E) S CCe
Proof. Because the processes are Gaussian random elements in C(M) by Lemma 18, their stochastic process RKHS given by Proposition 27 coincide with their Gaussian random element RKHS defined in van Zanten and van der Vaart [59]. Hence, for the non-truncated intrinsic and the extrinsic Matérn processes, the result follows by a direct application of Lemma 32 and Li and Linde [31], Theorem 1.2. For the intrinsic Matérn process truncated at Jn, it is not immediately clear that the constants C, so can be taken independent of n, so we go through the proof of Li and Linde [31], Proposition 3.1 see that this is, in fact, the case. We first need crude upper bound of the form
- In -InP((11// < <e) S cs
(190)
32

---

for some possibly large constant > 0. To get such a bound, we use Castillo et al. [12], Proposition 3 which shows the existence of universal constant C' such that for all E I min(1,40(11)
- In (MPP(IIll <e) - Cn(c) In (n(e) max (1, (1,(())
(191)
1/2
where of) = (Ey~n, 11 (mn)
and ra()) is defined in the following way in Castillo et al. [12), page 684 using auxiliary quantities L that are defined in Li and Linde [31], page 1562, namely n(E) = max 2 0: (1)(1) 2 E),
(192)
1/2 (+()) = inf | sy (mm())))) :f Σ2 Eghu
(193)
with (d) standing for the equality in distributions and the infimum being taken over all possible decompositions Σ,21 Eyhy with hy E C(M), e, being a sequence of IID N(0, 1) random variables, and the series being required converge uniformly almost surely."
The function f = Jo +1 (
X1-1))
w44/2
Ejfز-1 is valid decomposition. Therefore
1/2
Ju+1 Σ
2-4/2
4() S
2v
+ )-)1)1)
e (
(194)
R2
By the Sobolev Embedding Theorem and by Weyl's ' Law, given in Result 10, for every d/2 < 7 < v there exists constant C Cz.M such that for all J Z>o. allowing C' to change from line to line, we have
Ju +1
2
Ja+1
E Br-llp
Σ j=J
(1 (1+A)-1) Egfs-1
EC
Σ j=
(1+ (1+))-1)
vid/2
Egjs-1
(195)
(M)
(MM
Jn+1
2
(1 + \yー1))
(14d/22-- (u+d/2-7)
(196)
|=
Jn+1
C
j (1+2(u-))/1))
(197)
|=1
S 0 (1+2(2--7144 13/4
(198)
12J
< -
(199)
By choosing = this gives us of) < C independent of n. Moreover, by choosing . 2
again for comparison constant C' independent of n. this gives n(e) S Cs -, for C' independent of n. This implies using Castillo et al. [12], Proposition 3 that
In In (1 < <$) S cE
(200)
for c independent of n.
With this crude bound, we can now continue the proof of Li and Linde [31], Proposition 3.1. For this, we need metric entropy estimate. For this notice that for all J Z U foo} we have
BF" (0) B (0) Bll (MM)
(0), and therefore using Lemma 32, we have the metric entropy
estimate
In N (e,, BJf (0), (1112~))m))) S Ce w427z
(201)
for a constant C > independent of J. Therefore following the proof of proposition 3.1 in Li
and Linde [31] (with J 1) find - In TP (IIfIl <e) I Ce- - for every S co. where C, 5o > 0 are constants independent of n. We consider Σ₁21 unlike E5z0 frequently used above, in order follow the respective references.
33

---

Expressions for Pointwise Worst-case Errors
Let be kernel some abstract input domain X, and let Hg be the respective RKHS. Consider input values X S X and let of > O be the noise variance. Define
x. fic (1) Xxx(Xxx +z ( +e),
(202)
e (t) x - k(t, t) Kix + +o2) Kxe-
(203)
Proposition 34. With notation above
e (t)
sup /EH
E oo- mx,x,f,((t)
(204)
$1
Proof. To simplify notation, shorten (E to E and denote Kxx(Xxx First all, by direct computation,
o21)-
E Emx,x.f:st) (x(),, (m =x(x) La
(205)
(206)
Write
[B)((1) maxx,f,c((] S(1)2 +1
(207)
=((1)² ( + (x) T toaaa
(208)
= -(((t) - af(x))) + o0,0a
(209)
-K(()))) E Eack(zs))), 2
+ oa
(210)
j=1
As |1l1| suppem...lflll. 51 (1. implying sup fet 2 (g. f).
Ila we have
sup (I)(()) m.. =0) [|b(_]()))) FEH2 If 1fll,, 21
Σ ayk(oos
+ + | Es
(211)
k(t,t) -aKxx + aKxxa 0
(212)
akx:
= t(t,t) - aKxi = k(t, - KIx(XXX + 001) Kxi
(t,x(1)
(213)
We now move to the misspecified case. Consider the RKHS He for some other kernel e XxxX instead of Hg. Then, continuing from (210), write
R
sup Blf(t) - (t) ||ctt.)) Σ Eajeles SEH
| + +o
(214)
1. S1
The next question how compute the norm on the right-hand side. There is much hope of doing this exactly in the misspecified case: thus, we consider approximations. To this end, we take some large set of locations X' C X. Then we use -| ~ 9(X') C-1 Cxlx-9(X') for
g(-) c(t, Σ azc(I,,) -). As result, obtain the approximation
Etf(t) (x)x x(X) (x') + + ikcx((t) ul))(t) (215)
sup fe 11
where (()((t) was first introduced in Section 4
34

---

To compute spatial averages of this quantity, let 9c(+) = c(t,)) - EE= a)((zy,)), the same
before, but now with explicit dependence on Similarly, put ar K,x (Kxx -2 Then
9r (x(X)) Cx - ax = Cx Cxxx(Xxx + o² Kxx
(216)
gr (x') T Ix = (Crx- - xx -xxx- From here we can also deduce that
(217)
(k.x((t)=
Σ ((')) 9X
(218)
X
X
TX
(r(gx.(X)) Cxxx
(219)
X
where gx (x') = Xx-x XxxxxXxxxXXXXXXXXXX + - -Kxx"
G Full Experimental Details
All of our kernels computed using GPJAX [38] and the GEOMETRIC KERNELS library. 10 We use three manifolds, each represented by a mesh: (i) a dumbbell-shaped manifold represented as a mesh with 1556 nodes, (ii) sphere represented by icosahedral mesh with 2562 nodes, and (iii) the Stanford dragon mesh, preprocessed to keep only its largest connected component, which has 100179 nodes. For the sphere, we also considered finer icosahedral mesh with 10242 nodes, but this was found to have virtually no effect on the computed pointwise expected errors.
We use extrinsic Matérn and Riemannian Matérn kernels with the following hyperparameters: of 1 and of 0.0005. For the truncated Karhunen- expansion, we used 500 eigenpairs obtained from the mesh. We selected smoothness values to ensure corm-equivalence of the intrinsic and extrinsic kernels' reproducing kernel Hilbert spaces, which was 5/2 for the intrinsic model, and μ = 5/2 + d/2 for the extrinsic model, where d is the manifold's ' dimension. We used different length scales for each manifold: x 200 for the dumbbell, x 0.25 for the sphere, and 0.05 for the dragon, selected to ensure that the Gaussian processes were neither approximately constant, nor white-noise-like. We considered data sizes of N = 50, N 500, and for the dumbbell, sphere, and dragon, sampled uniformly from the mesh 's nodes, which in each case distribution of points across the manifold. Finally, for the extrinsic pointwise error approximation, we used a subset X' uniformly sampled from each mesh's 's nodes, of size equal to the data size. For each respective test set, we used the full mesh. Each experiment was repeated for 10 different seeds.
1000, respectively,
resulted
To set the length scales for the extrinsic process, we used maximum marginal likelihood optimization on the full data, except for the dumbbell whose full data size is small and for which we instead generated larger set consisting of 500 points. We optimzied only the length scale, leaving all other hyperparameters fixed. We used ADAM with learning rate of 0.005, and an initialization equal to the length scale k the intrinsic model, except for the dumbbell where this lead to divergence and we instead used an initial value of x/4. We ran the optimizer for total of 1000 steps. With these settings, we found empirically that maximum marginal likelihood optimization always converged.
See READTHEDOCS and HITIPS//GEEMETRC KERNELS GITHUB.IO.
35


Reviews:
[{"review_id": "TZDCEEMtol", "reviewer": "Reviewer_yf8T", "summary": "This paper investigates the theoretical properties and performance of Gaussian processes in machine learning, particularly when applied in geometric settings such as Riemannian manifolds. It compares intrinsic and extrinsic methods, with the former directly formulated on the manifold of interest and the latter requiring a higher-dimensional Euclidean space embedding. The research derives posterior contraction rates for three primary geometric model classes and shows that all three can lead to optimal procedures, given certain conditions. Empirical experiments support these theoretical findings, demonstrating better performance by intrinsic models in small-data regimes.", "questions": "NA", "limitations": "NA", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "NA", "weaknesses": "NA", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for the accurate summary of our work and for your review! If you have any further questions or comments, please feel free to post them - we are happy to provide further information or clarification where needed. If not, thank you very much for your time in reading our work!"}, {"review_id": "maKz1zon0a", "reviewer": "Reviewer_qZfm", "summary": "This paper establishes bounds on the contraction rate of matern processes on Riemannian manifolds. The authors study three variants:\n1) intrinsic matern process\n2) truncated intrinsic matern process\n3) extrinsic matern process\nand show that in each case, the optimal contraction rate can be achieved, which matches the Euclidean case. \n\n", "questions": "Regarding the smoothness parameter nu:\n1) Can the authors provide intuition for why nu (as used in (7) and (10)) is the smoothness parameter (i.e. smoothness in what sense)?\n2) The Riemannian matern process (7) and the Euclidean matern process (10) appear quite different, can the authors explain the analogy between these two processes, and specifically, why is nu, as used in (7), comparable to nu as used in (10)? In particular, can the authors explain the comment from line 275 \"the intrinsic Matérn process, its truncated version and the extrinsic Matérn process all possess the same posterior contraction rates\", and how these two rates are comparable under the different assumptions?\n3) All theorems assume nu > d/2. Is this a standard assumption? How necessary is this assumption, and the theorems not work when nu is small? (and I have the a similar question for the beta > d/2 assumption in assumption 3 as well)\n\nRegarding the intrinsic matern process:\n4) expression (7) involves an sum over eigenfunctions of the laplace beltrami operator. The authors do mention that this sum can be truncated, but for a general manifold, it seems like even computing a single eigenfunction can be quite expensive. Can the authors comment on how this is done computationally (and what is the cost), when the manifold is an arbitrary one, e.g. the dragon?\n\n5) In figure 2, the authors give a dumbbell example which highlights the difference between intrinsic and extrinsic kernels -- one important difference seems to be that two points can be far away in manifold distance, but close in euclidean distance (under the embedding). Consequently, a function may have a small lipschitz constant wrt manifold distance, but a huge lipschitz constant wrt euclidean embedding distance. Intuitively, why is this not reflected in the contraction rates? Is it because the assumptions made do not care about things like lipschitz smoothness? (related to my earlier question of what is the meaning of nu?)\n\n\nRegarding the bound:\n\n6) What is contained in the constant C in the theorems, is this a universal constant? Or does it depend polynomially/exponentially on any problem parameters?", "limitations": "n/a", "rating": 7, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 3, "strengths": "This is a fundamental problem, and it is remarkable that the authors are able to prove the same optimal contraction rates for both the intrinsic and extrinsic Riemannian matern processes. Furthermore, the manifold hypothesis has been receiving increasing attention recently, and so the analysis of Gaussian processes on manifolds is well motivated.\n\nI also appreciate the examples that the authors provided for illustrating the difference between intrinsic and extrinsic processes.", "weaknesses": "see questions below", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for taking the time to review our work! Thank you for the encouraging comment that work addresses a fundamental research problem! Below we address the questions, some of which we had to partially quote due to character limits:\n\n---\n\n*\".. intuition for why nu is the smoothness parameter ..\"*\n* Thank you for the question! Indeed, it may not appear obvious that $\\nu$ in the definitions of the priors can be interpreted as a smoothness parameter. This stems from the sample path regularity of the processes: it is shown in [44] that the samples of the Euclidean Matérn process (and therefore also those of the extrinsic process) are $\\alpha$-Holder for every $\\alpha<\\nu$, whereas this property for the intrinsic Matern process is **precisely the content of our Lemma 27**.\n\n*\".. explain the analogy between these two processes  (Riemannian vs. Euclidean) ..\", \".. why is nu, as used in (7), comparable to nu as used in (10)\"*\n* Thank you for bringing this important point of confusion, which is likely to be shared by other readers, to our attention. To compare the two processes, note that both can be represented as **(weak) solutions of the same stochastic partial differential equation**, defined over Euclidean space and the manifold, respectively. This is the viewpoint from which the intrinsic Matérn process was studied in [8] and other prior works, and is why both are called Matérn processes. \n* Regarding $\\nu$: see the remark above, but also note that the Euclidean Matérn process on R^d has an RKHS norm equivalent to the Sobolev space $H^{nu+d/2}(R^d)$, and the same is true for the intrinsic Matern by our Lemma 23 where the Sobolev space $H^{nu+d/2}(M)$ is defined using the Bessel potential space formulation.\n\n*\".. how these two (contraction) rates are comparable under the different assumptions?\"*\n* Thank you for the comment! Let us clarify this: we fix a **single, common data generating process**, namely nonparametric regression with random design and a fixed unknown regression function $f_0$, and, for **three different Bayesian models**, we compare the asymptotic contraction/convergence rate of the posterior distribution towards the true regression function $f_0$. The rates all depend on the intrinsic dimensionality of the data $d$, the smoothness parameter of our prior processes $nu$, and the smoothness of the true regression function $f_0$, and for the models considered end up equal up to constants. To improve clarity, we will update the manuscript to further emphasize this.\n\n*\".. nu > d/2. Is this a standard assumption? .. *\n* This is a good question! **Yes, this is relatively standard:** $\\nu>d/2$ is also used in the prior work [44], where the Euclidean counterparts of our results are proved. From a technical standpoint, this is needed in order to get a convergence rate under the $L^2(p_0)$ norm from a convergence rate at the input locations. Removing this assumption would be an interesting research question, but we suspect that it does not particularly involve geometry specific to the manifold setting which is our focus here.\n\n*\".. comment on how this (obtaining Laplace-Beltrami eigenfunctions) is done computationally ..\"*\n* Thank you for the comment - this is actually the main computational challenge in working with kernels on manifolds. In practice one can rely on two different sets of techniques. The first one is to **discretize/mesh the manifold and solve for the eigenpairs of a large sparse matrix** as is done in [8], although this inevitably introduces numerical errors on top of the asymptotic contraction rates that we present here. The second way is to rely on **algebraic techniques based on symmetries which makes an exact computation possible for a large class of manifolds** of interest - see [2,3]. \n\n*\".. a function may have a small lipschitz constant wrt manifold distance, but a huge lipschitz constant wrt euclidean embedding distance. Intuitively, why is this not reflected in the contraction rates? ..\"*\n* This is an *extremely relevant point*: curvature and more general distortions between geodesic and Euclidean distances introduce undesirable behaviors of the extrinsic processes when compared to the intrinsic ones. In our analysis we are however incorporating the **same smoothness assumptions** on the regression function in both the intrinsic and the extrinsic case and nonetheless find **equivalent contraction rates** for both priors: the crucial fact is that the contraction rates presented here are all asymptotic. Hence, we highly suspect a **very bad dependence of the constant on the curvature and embedding** in the extrinsic case which would explain the differences in performance that we can witness in the small data regime - in fact, we believe our work provides an **excellent motivation to initiate future work on non-asymptotic contraction analysis** that we believe may be needed to capture these differences.\n\n*\".. constant C in the theorems, is this a universal constant? Or does it depend polynomially/exponentially on any problem parameters?\"*\n* Thank you for pointing that out! As discussed above, the value of $C$ is precisely what makes the extrinsic and intrinsic processes non-equivalent in practice when it comes to the performances. An inspection of the proof shows that $C$ depends on (1) $d,D$ (for the extrinsic prior), (2) the prior process hyperparameters, (3) $\\beta$, (4) $M$, (5) the distribtuion $p_0$ over $x$-values, (6) $\\sigma_\\epsilon^2$, and (7) the Sobolev/Holder norms of the true regression function $f_0$. As mentioned in line 295 the constant in the case of the extrinsic process is expected to be bigger than the one for the intrinsic process, especially when the distortion between geodesic and Euclidean distance is high. We suspect in particular that the operator norm of the trace and extension operators between Sobolev spaces to be big in this case. We will add further clarifications on the constant in the revised version."}, {"review_id": "agHcAQabC2", "reviewer": "Reviewer_woEb", "summary": "This paper studies the contraction rate(s) for both the intrinsic and extrinsic Mat\\'ern Gaussian process in compact Riemannian manifold. The authors proved that the (optimal) rate in both cases is $\\frac{2 \\min(\\beta, \\nu)}{2 \\nu +d}$, where $\\nu$ is the smooth parameter of the Mat\\'ern process, $\\beta$ is the regression function smoothness class, and $d$ is the dimension. The authors also showed with examples that the geometric models outperform the non-geometric ones through empirical error analysis. \n", "questions": "While the paper is generally well written, I feel the authors need to address the following points:\n\n(1) I catch the idea of the intrinsic vs extrinsic modeling quickly. But the authors may point out a reference for this terminology (is this inspired from [15])? \n\n(2) p.3, line 80: I guess the authors mean $f \\sim GP(m, k)$ instead of $GP(0,k)$.\n\n(3) p.3, line 115-116: the notation $C^{\\beta}$ for the H\\\"older class is nonstandard. Does this mean $\\mathcal{C}^{0, \\beta}$? I guess not since $\\beta > \\frac{d}{2}$ meaning that it can be larger than $1$. The authors may need to clarify this.\n\n(4) In Result 1, it is assumed that $f$ is mean zero. Is this also assumed in Theorems 5 and 6? It is known that for the Mat\\'ern processes in the Euclidean space, the mean function may raise the identifiability issue (see Stein's Interpolation of Spatial Data, or Tang, Zhang and Banerjee's paper On identifiability and consistency of the nugget in Gaussian spatial process models). There are also related discussions in [24]. The authors may need to clarify, and add a few sentences in the paper.\n\n(5) p.4, Assumption 3: the authors may \"indicate\" $\\sigma_\\epsilon$ is known earlier, as this is important for the Bayesian workers.\n\n(6) Regarding the use of Theorems 5 and 6: it is clear from the rates that one should take $\\nu = \\beta$ (i.e. if one can identify the function class, and set the same smoothness parameter in the Mat\\'ern process). However, we never know exactly $\\beta$. The question is whether there is an adaptive way to select $\\nu$. (Of course, this may be discussed in another paper. I only want to bring this question to the authors.)\n\n(7) Regarding Theorem 8: the rate is $\\frac{2 \\min(\\beta , \\nu)}{2 \\nu +d} = \\frac{2 \\nu}{2 \\nu +d}$ if $\\beta = \\nu$. This rate also appears in the prediction based on BLUE (best linear unbiased predictor) on the context of MLE. See e.g. Tang, Zhang and Banerjee's paper On identifiability and consistency of the nugget in Gaussian spatial process models, JRSS-B, 2021, page 1055. The authors may want to mention this as well.\n\n(8) There are a few more references that the authors may want to add.\n\n(a) Stein's book Interpolation of Spatial Data is one of the main references on the Mat\\'ern Gaussian processes. \n\n(b) Tang, Zhang and Banerjee's paper On identifiability and consistency of the nugget in Gaussian spatial process models studies the Mat\\'ern process (with nugget) in the Euclidean space, where the identifiability issue occurs as pointed out in (7). This is related to the setting of Theorem 8. Arafat, Porcu, Bevilacqua and Mateu's paper Equivalence and orthogonality of Gaussian measures on spheres, Journal of Multivariate Analysis, 2018 is also relevant. \n\n(c) Regarding the prior on $\\sigma^2_\\epsilon$, one such model is the conjugate Bayesian linear model in Banerjee's paper Modeling massive spatial datasets using a conjugate bayesian linear modeling framework, Spatial Statistics, 2020; this model was analyzed in Zhang, Tang and Banerjee's Exact Bayesian geostatistics using predictive stacking, arXiv:2304.12414, Section 2. I believe it should be possible to generalize the results in this paper to the conjugate bayesian linear model.", "limitations": "NA", "rating": 8, "confidence": 5, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "The results in this paper are novel and enlightening. Up to minor typos, I enjoy reading the paper. It is concerned with the fine topic of Gaussian processes on manifolds, and on why (and how) this kind of modeling is valuable through quantitative posterior contraction analysis. The main contribution is the (optimal) contraction rate of both intrinsic and ambient Mat\\'ern processes on compact manifolds in the nonparametric setting, the analysis of which is different from that in the Euclidean setting as the definition of the Mat\\'ern processes on manifolds is subtle (though the rates are the same). It is also shown how the underlying geometric analysis outperforms by numerical experiments. ", "weaknesses": "As for every (good) paper, there are always plenty of things remained to be done. For instance,\n\n(1)  In the current analysis it is assumed that the \"nugget\" $\\sigma_\\epsilon$ is given. In Bayesian setting, it would be interesting to know what happens if we put a prior on it (and more importantly, what prior to put). \n\n(2) The results, implicitly, compute the interpolation errors (as $L^2$ error based on $p_0$). What about extrapolation errors (that is, to predict an \"outside\" the domain)? This is also important in most geostatistics problems. \n\n(3) For the numerical experiments, the authors only consider the synthetic examples, e.g. dragon, sphere. It is interesting to see if the intrinsic and extrinsic modeling can bring a big difference for some real data. One such example is provided in [14].", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for your thorough review of our work and for the very encouraging comments! Below we address key questions:\n\n**Further work:**\n\n(1) *Nugget and prior on $\\sigma_\\epsilon$*\n* Thank you for this question! In our work, the main reason we assume $\\sigma_\\epsilon$ is fixed is because in the Euclidean case other work has studied conditions on priors that ensure minimax optimality - at least up to ln n factor, as is done in [17] - **we expect these to be similar in both the Euclidean and geometric settings**. We therefore opted to not do this to avoid the paper becoming too long, but will amend the manuscript to add references so that readers interested in this case can find the relevant papers.\n\n(2) *Interpolation vs. extrapolation, $p_0$-norm vs empirical norm*\n* Thank you for the comment! We would argue that **our results do, in an appropriate technical sense, extrapolate** outside of the design points, since $p_0$ is an absolutely continuous measure - as opposed to, for instance, to the empirical measure of the data. Our proof technique uses the assumption that $p_0$ is lower-bounded, which implies that we can control convergence over the whole manifold, by moving between different distributions by changing the constants in the bound. Note that, if one considers regions where $p_0$ takes small values, the constant may degenerate, so this assumption does have limitations - but, we expect these properties are similar to the Euclidean case and not particularly specific to the manifold setting.\n\n(3) *Real data examples*\n* Thanks - this is a great point! We opted to focus on synthetic examples for simplicity, since this allows us to better control the moving parts in the experiment. However, we also wanted to note that in addition to [14] a similar performance difference was also observed in [12] in the context of medical data with a slightly-different prior. We will add a few remarks mentioning this to the experimental section to point reviewers toward these papers.\n\n**Questions:**\n\n(1) *Intrinsic vs. extrinsic terminology*\n* Thanks for this question! We called the processes this because the concepts “extrinsic” and “intrinsic” mirror the distinction between intrinsic and extrinsic properties/quantities **in the sense of differential geometry** (and mathematics more generally) where the former refers to concepts not needing any kind of embedding to be defined, while the latter refers to ones that need to be expressed through an embedding in a higher dimensional space/object. This is exactly the difference in how the two Matérn processes are constructed.\n\n(2) *Thanks for spotting the typo!*\n\n(3) *Notation for Hölder spaces*\n* This is a great point - thank you for spotting this! For $\\gamma = k + \\alpha$ with integer $k >= 0$ and $0 < \\alpha <= 1$, define $CH^\\gamma$ to be the space of k times differentiable functions with $k$th derivative being $\\alpha$-Hölder. We've **changed the notation to $CH$** to avoid confusion with ordinary smooth functions.\n\n(4) *Prior mean and identifiability*\n* Thank you for this observation! The mean of the prior processes is indeed kept fixed at 0, as is commonly done for proving contraction rates for Gaussian processes. This choice actually leads to optimal contraction rates, but could more generally be relaxed. Regarding identifiability, note that we are not trying to identify covariance parameters given a sample of a Matérn process in an infill asymptotic regime, but rather to show contraction of the posterior towards a fixed regression function. **Our nonparametric regression model's parameter is therefore identifiable**: the probability distribution is only indexed by the regression function, and our results imply the existence of a consistent estimator for it, by taking for instance the posterior mean of the process. We will add a few sentences and references on this particular point.\n\n(5) *Emphasis on $\\sigma_\\epsilon$*\n* This is a good idea - we'll add a sentence on this!\n\n(6) *Adaptive selection of $\\nu$*\n* This is an *excellent* point! We agree adaptivity to the smoothness of f_0 is a natural next step. For the intrinsic Matérn process adaptivity can be achieved by standard techniques that are not specific to the geometric setting. One way would be to follow the approach of Kirichenko & Van Zanten in “Estimating a smooth function on a large graph by Bayesian Laplacian regularization” where a multiplicative scaling parameter is introduced in the definition, and on which another prior is placed. Another solution could be to consider a truncated intrinsic Matérn process with a prior on the truncation, similar to Waaij & Van Zanten in “Full adaptation to smoothness using randomly truncated series priors with Gaussian coefficients and inverse gamma scaling”. For the extrinsic Matérn process, the problem is more complicated, but is tackled for the RBF kernel by Yang & Dunson in “Bayesian Manifold Regression” where they show that they can achieve adaptivity by placing a prior on the length scale of the prior process: **we believe a result of the same flavor could be shown in our context, although our proof technique is fundamentally different**. Studying adaptivity is, however, difficult and we believe should be the focus of a follow-up paper, but we are happy to add more comments on this in the revised version.\n\n(7) *Connection with other asymptotic rates*\n* For $\\beta = \\nu$, namely the less-realistic case where the right smoothness has been chosen in advance, we indeed recover the same contraction rate as the one you mention. This makes sense, as this is the minimax optimal rate of estimation of $f_0$ in this model. **It's a great idea to point this out with appropriate references** - thanks for bringing this to our attention!\n\n(8) *References*\n* Thank you very much for bringing these references to our attention! We will add them to help the reader frame our results with respect to identifiability and current research on Bayesian models."}, {"review_id": "0j0scMoDHl", "reviewer": "Reviewer_dGg7", "summary": "The paper concerns Gaussian processes on manifolds. The authors present theorems for contraction rates for Matern Gaussian processes defined intrinsically and extrinsically through and embedding in a higher-dimensional Euclidean space. The authors shows that the rates are asymptotically equal in the two settings. Additionally, they treat the case of finitely truncated series expansion of the kernels to get a similar rate. Finally, it is shown experimentally that it can still be beneficial to work with the intrinsic geometry in the small sample domain.", "questions": "no questions", "limitations": "yes", "rating": 8, "confidence": 3, "soundness": 4, "presentation": 4, "contribution": 4, "strengths": "- extremely well-written paper. I found the exposition very clear\n- novel theoretical results\n- interesting findings (in line with what the authors state, I would not have expected the manifold dimension to show up in the extrinsic case)\n- empirical study to cover the small sample case", "weaknesses": "- since the theorems are not in the main paper, one could perhaps consider if a longer format than NeurIPS (a journal paper) would be more suitable for the paper from a presentation point of view", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "Thank you very much for your review, and especially for the comment that you found our exposition \"very clear\" in spite of the technical nature of the subject! Below we comment on one of the points.\n\n---\n\n*\"since the theorems are not in the main paper, one could perhaps consider if a longer format than NeurIPS (a journal paper) would be more suitable for the paper from a presentation point of view\"*\n* While it certainly would have also been possible for us to write up our work in a longer format compared to NeurIPS such as JMLR, we believe that the NeurIPS format is also appropriate here, because the relatively short page limit means that a paper like ours needs to focus its main body on summarizing the main theoretical results, rather than on technical details in proofs. It is valuable to have papers like this, because they provides a picture of the state of affairs which are also readable by non-experts - we think this is particularly important for posterior contraction analysis, which tends to involve harder-to-parse mathematics than other areas. Therefore, we believe the NeurIPS format used by this work effectively complements other, much longer and more technical-detail-focused papers on related questions which can be found in the literature."}]

Author rebuttal:
We would like to thank the referees for their summaries and pertinent observations that will help us to improve the final version of the paper. Most of the comments were about the format, additional references, identifiability and insights regarding the definitions of the different processes; we have responded to each reviewer in detail.

Meta-review: The authors analyze Gaussian processes on Riemannian manifolds. This question is of significant interest in many areas of machine learning. They bound the posterior contraction rate (i.e. rate at which posterior converges to true data distribution) of intrinsic, truncated intrinsic, and extrinsic matern processes. (Interestingly the optimal rate is the same for intrinsic and extrinsic processes)

The reviewers all agree that the technical result is novel and interesting. Empirical results complement their results nicely, highlighting potential future directions for sharper bounds which will reveal the reveal the difference between intrinsic and extrinsic methods in certain settings.
Decision: Accept (spotlight)