Uncertainty Quantification of Spline Predictors on Compact Riemannian Manifolds

Uncertainty Quantification of Spline Predictors on Compact Riemannian Manifolds Charlie Sire ∗ , Mike Pereira Center for Statistics and Images, Mines P aris, PSL University , F ontainebleau, 77300, F rance Abstract T o predict smooth physical phenomena from observ ations, spline interpolation provides an inter - pretable framework by minimizing an energy functional associated with the Laplacian operator . This work proposes a methodology to construct a spline predictor on a compact Riemannian manifold, while quantifying the uncertainty inherent in the classical deterministic solution. Our approach lev erages the equiv alence between spline interpolation and univ ersal kriging with a spe- cific co variance kernel. By adopting a Gaussian random field frame work, we generate stochastic simulations that reflect prediction uncertainty . Howe ver , on compact manifolds, the cov ariance kernel depends on the generally unkno wn spectrum of the Laplace-Beltrami operator . T o address this, we introduce a finite element approximation based on a triangulation of the manifold. This leads to the use of intrinsic Gaussian Markov Random Fields (GMRF) and allo ws for the incor- poration of anisotropies through local modifications of the Riemannian metric. The method is validated using a temperature study on a sphere, where the operator’ s spectrum is known, and is further extended to a test case on a c ylindrical surface. K eywor ds: Uncertainty Quantification, Spline interpolation, Riemannian manifold, Intrinsic Gaussian Markov Random Fields, Uni versal Kriging ∗ Corresponding author Email addr ess: charlie.sire@minesparis.psl.eu (Charlie Sire) 1. Introduction Thin-plate spline interpolation [1] is one of the most widely used methods for predicting spa- tially distributed data. It consists in constructing an interpolating function that exactly matches the observed v alues at the measurement locations while being as smooth as possible according to a precise mathematical criterion. This notion of optimality , defined through the minimization of a roughness functional, is particularly well suited to settings with few observ ations in which the underlying phenomenon is expected to vary smoothly over space. Howe ver , thin-plate splines yield a single deterministic prediction over the domain and do not directly provide a measure of uncertainty . In many applications, uncertainty quantification is crucial, for instance to assess risk or compute predictiv e probabilities, which moti vates the need for a full predictiv e distribu- tion rather than a point estimate [2]. From this perspective, the spline predictor can be shown to be equiv alent to the uni versal kriging predictor, as first established by Matheron [3] and later exploited for data defined on R [4] and on R d [5, 6]. In this framew ork, the observations are viewed as a realization of a random field with an unknown mean function and a prescribed covariance structure, that are related to the eigen values and eigen vectors of the Laplacian operator − ∆ . Univ ersal kriging then relies on a statistical model to provide not only a prediction b ut also an associated predicti ve v ariance at e very point in the domain. T o fully characterize this predictive uncertainty , a Gaussian assumption is typically adopted for the prior random field, leading to Gaussian process re gression (GPR) [7], which enables computation of the posterior distrib ution of the random field, that is, the conditional distribution giv en the observed v alues at the measurement locations. While the predicti ve mean and variance obtained from kriging and GPR are identical, the latter frame work additionally allows for the simulation of posterior realizations of the random field [8]. These simulations provide a richer representation of uncertainty and facilitate the visualization and quantification of plausible scenarios under the assumed statistical model. Thus, exploiting Gaussian process regression with the appropriate covariance kernel makes it possible to incorporate uncertainty into the spline prediction problem. Howe ver , when the domain of interest is a non-planar surface, sev eral practical challenges arise. Studies on such non-Euclidean domains appear in a variety of conte xts, for e xample in fluid mechanics when analyzing flows around cylinders [9, 10], or in neuroimaging when collecting data on the brain [11]. In theory , the equiv alence between spline interpolation and universal kriging prediction still holds on compact Riemannian manifolds, as formally established in [12] using the frame work of Reproducing Kernel Hilbert Spaces (RKHS), building on the pioneering work of W ahba on the sphere [13]. On manifolds, the mean function and covariance kernel un- derlying this equiv alence are derived from the spectral decomposition of the Laplace-Beltrami operator , a natural generalization of the Euclidean Laplacian. In practice, ho wev er , the eigenv al- ues and eigenfunctions of the Laplace-Beltrami operator are generally unkno wn for surf aces that are neither planar nor spherical. This lack of spectral information pre vents the exact construction of spline predictors and their associated co variance structures for data defined on such manifolds. T o overcome this limitation, and as proposed in [14], one possible approach is to rely on a finite element approximation of random fields defined on compact manifolds, thereby working with a finite set of approximate eigen values and eigenfunctions. This type of approximation was first introduced by Lindgren [15] and later extended by Pereira [16] to incorporate local anisotropies into the co variance structure. This framew ork can then be adapted to construct spline predictors with preferred directions of dependence, which can significantly improve prediction 2 accuracy when the underlying phenomenon is anisotropic. Howe ver , a key di ffi culty arising in the finite-element approximation of spline predictors is that it requires working with a Gaussian Markov random field (GMRF , see [17]) whose cov ari- ance matrix is not in vertible, referred to as an intrinsic GMRF . This type of model, inv estigated in particular in [18], complicates both the formulation of the classical kriging equations and the generation of simulations. As shown in [19] and thoroughly studied in [20], such simulations can be obtained through linear conditioning of non-intrinsic GMRFs. This strategy enables an extension of the w ork of [14], which focused exclusiv ely on computing an approximate solution to the spline interpolation problem, without addressing the associated uncertainty quantification of the resulting predictions. Therefore, the contribution of this work is to propose a statistical model for the spline inter- polation problem, with the aim of providing both spline-based predictions at ev ery point of the in vestigated compact Riemannian manifold and a quantification of the associated uncertainty . A strategy for generating conditional simulations of the underlying random function is also pro- posed, enabling the visualization of predictive scenarios arising from spline interpolation prob- lems. The paper is organized as follows. Section 2 introduces the theoretical solution to the spline interpolation problem on compact Riemannian manifolds and establishes its equiv alence with kriging. Section 3 presents the finite element approximation of this solution. Section 4 then details the strategy for generating conditional simulations of the target intrinsic GMRF , while Section 5 addresses the introduction of anisotropy and the estimation of the unknown parame- ters. Finally , the application test cases and corresponding results are presented in Section 6, and Section 7 summarizes the contributions and outlines potential directions for future w ork. 2. Spline interpolation solution on compact manif olds In this manuscript, the inv estigated domain is a smooth compact and connected manifold M of dimension d , equipped with a Riemannian metric g , and − ∆ denotes the Laplace-Beltrami operator on ( M , g ). The integral of a function f ov er M is denoted by R M f d µ g , where d µ g is the canonical measure associated to ( M , g ) and ∂ M denotes its boundary . The spaces L 2 ( M ) and H 2 ( M ) denote the set of square-integrable functions and the Sobolev space of order 2, re- spectiv ely . The objective of this article is to obtain a prediction of a phenomenon based on n observations denoted by y = ( y i ) n i = 1 at points S = ( s i ) n i = 1 with s i ∈ M , 1 ≤ i ≤ n . 2.1. F ormulation of the spline pr ediction pr oblems The spline interpolation problem on M , denoted by P 0 , consists in finding the smoothest function, in a precise mathematical sense, that interpolates the observed data. More precisely , it consists in solving the following v ariational problem [21, 22]. Definition 1 (Spline interpolation problem on M ) . P 0 the spline interpolation pr oblem on M consists in finding a function u ∈ H 2 ( M ) such that            E ( u ) = min v ∈ H 2 ( M ) Z M | ∆ v | 2 d µ g u ( s i ) = y i , 1 ≤ i ≤ n . (1) 3 When the observed data are corrupted by noise, seeking an exact interpolating solution is no longer appropriate. In this case, the smoothing spline problem, denoted by P τ , is more suit- able, as it consists in minimizing the energy functional penalized by the prediction error at the observation points. Definition 2 (Smoothing spline problem on M ) . P τ the smoothing spline pr oblem on M consists in finding, for τ > 0 a function u ∈ H 2 ( M ) minimizing 1 n n X i = 1 ( y i − u ( s i ) ) 2 + τ 2 Z M | ∆ v | 2 d µ g (2) The analytical solutions of these problems are very similar , and are related to the spectrum of the operator − ∆ . 2.2. Analytical solution to the spline pr ediction pr oblems T o deriv e the solutions to the problems P 0 and P 1 , we rely on the Spectral Theorem [23], recalled in Appendix A, which provides the countable spectrum of − ∆ , denoted by ( λ k , ϕ k ) k ∈ N . Here, ( λ k ) k ∈ N ⊂ R + are the eigenv alues and ( ϕ k ) k ∈ N ⊂ C ∞ ( M ) the associated eigenfunctions, which form an orthonormal basis of L 2 ( M ). The eigen value λ 0 = 0 corresponds to constant functions [24], with ϕ 0 = 1 R M d µ g . As explained in [12], the theory of RKHS is needed to obtain the solution of problems of the type P τ for τ ≥ 0 . Some details were recalled in [14], but are not ev oked here. W e simply introduce the kernel that will be central in the solution, which is defined as K 1 ( s 1 , s 2 ) = X k ∈ N ⋆ 1 λ 2 k ϕ k ( s 1 ) ϕ k ( s 2 ) , ( s 1 , s 2 ) ∈ M 2 , (3) and the solution to the spline prediction problems is provided by Proposition 1. Proposition 1 (Spline prediction [25]) . The solution to the pr oblem P τ , which corr esponds to the interpolating spline when τ = 0 and to the smoothing spline when τ > 0 , is given by u ⋆ τ ( s ) = a τ ϕ 0 ( s ) + k 1 ( s ) ⊤ K − 1 1 ,τ  y − a τ ϕ obs 0  , (4) with                        K 1 ,τ = h K 1 ( s i , s j ) i 1 ≤ i , j ≤ n + τ 2 I n k 1 ( s ) = ( K 1 ( s 1 , s ) , . . . , K 1 ( s n , s ) ) ⊤ ∀ s ∈ M ϕ obs 0 = ( ϕ 0 ( s 1 ) , . . . , ϕ 0 ( s n ) ) a τ =  ( ϕ obs 0 ) ⊤ K − 1 1 ,τ ϕ obs 0  − 1 ( ϕ obs 0 ) ⊤ K − 1 1 ,τ y (5) The inv ertibility of the matrix K 1 ,τ is discussed in Appendix B. For τ = 0 , K 1 ,τ is positive definite if and only if an y observation v ector y can be interpolated by a function in L 2 ( M ), while the in versibility is guaranteed when τ > 0 . 4 2.3. Equivalence with kriging The solution u τ coincides with the univ ersal kriging predictor [26, 27], with mean a ϕ 0 , where a ∈ R is an unknown parameter, cov ariance kernel σ 2 K 1 , and observation noise variance σ 2 τ 2 when τ > 0. Here, σ 2 > 0 is the unkno wn variance. Note that since ϕ 0 is constant in this setting, the model reduces to ordinary kriging [28], which assumes a constant but unkno wn mean. Under these assumptions, and with predictive mean equal to u τ ( s 0 ) at a new point s 0 ∈ M , the associated kriging variance is gi ven by V UK ( s 0 ) = σ 2 h v ( s 0 ) + K 1 ( s 0 , s 0 ) − k 1 ( s 0 ) ⊤ K − 1 1 ,τ k 1 ( s 0 ) i (6) with v ( s 0 ) = h ϕ 0 ( s 0 ) − k 1 ( s 0 ) ⊤ K − 1 1 ,τ ϕ obs 0 i ⊤   ϕ obs 0  ⊤ K − 1 1 ,τ ϕ obs 0  − 1 h ϕ 0 ( s 0 ) − k 1 ( s 0 ) ⊤ K − 1 1 ,τ ϕ obs 0 i . The term σ 2 v ( s 0 ) reflects the uncertainty associated with the unkno wn coe ffi cient a appearing in the mean of the random field, while σ 2 h K 1 ( s 0 , s 0 ) − k 1 ( s 0 ) ⊤ K − 1 1 ,τ k 1 ( s 0 ) i corresponds to the pre- dictiv e variance of the simple kriging predictor , that is, when the mean is assumed to be kno wn. This analogy with univ ersal kriging naturally pav es the way for uncertainty quantification in spline prediction. Howe ver , to obtain the full predicti ve distribution rather than just the mean and variance, the idea is to introduce Gaussian assumptions and to w ork with Gaussian random fields (GRFs) Z and E . Since the parameter a governing the mean is unknown, the Bayesian kriging framew ork [29] is adopted by modeling it as a random variable A . The observations y are then interpreted as noisy realizations of the GRF Z such that, Y i = Z ( s i ) + στ E i , 1 ≤ i ≤ n where • ( Z | A = a ) is a GRF with mean a ϕ 0 and cov ariance kernel σ 2 K 1 , for all a ∈ R • A ∼ N ( µ, α ) , which is the a priori distribution of the unknown coe ffi cient of the mean • E = ( E i ) n i = 1 is a centered Gaussian vector with Co v( E ) = I n , independent of A and Z . This corresponds to the a priori information on the model, where Z again represents the underlying phenomenon of interest and στ E accounts for an additional measurement noise when τ > 0. Under the Gaussian assumptions, the posterior distribution of the random field Z can be computed explicitly , namely ( Z | Y obs = y ) , (7) where Y obs denotes the random vector of observations [ Y 1 , . . . , Y n ] ⊤ . This conditional random field is itself a GRF , and its distrib ution is therefore entirely characterized by its mean and covari- ance kernel. As explained in [30] and recalled here in Proposition 2, its behavior when α → ∞ is of particular interest. Proposition 2 (Con vergence of the conditional GRF .) . When α , the variance of A, tends to infinity , corresponding to a non-informative prior , the conditional field ( Z | Y obs = y ) con verg es to a GRF whose mean and covariance coincide with the universal kriging predictor . They ar e given by          lim α →∞ E ( Z ( s ) | Y obs = y ) = u ⋆ τ ( s ) lim α →∞ Cov  Z ( s ) , Z ( s ′ ) | Y obs = y  = σ 2 h c ( s , s ′ ) + K 1 ( s , s ′ ) − k 1 ( s ) ⊤ K − 1 1 ,τ k 1 ( s ′ ) i (8) 5 with c ( s , s ′ ) = h ϕ 0 ( s ) − k 1 ( s ) ⊤ K − 1 1 ,τ ϕ obs 0 i ⊤   ϕ obs 0  ⊤ K − 1 1 ,τ ϕ obs 0  − 1 h ϕ 0 ( s ′ ) − k 1 ( s ′ ) ⊤ K − 1 1 ,τ ϕ obs 0 i . The proof is pro vided in Appendix C. Note that the a priori mean µ of the random v ariable A does not a ff ect the results when α → ∞ . Therefore, without loss of generality , we set µ = 0 for the remainder of the analysis. The objecti ve of the study is then to in vestigate this conditional GRF , with a particular focus on generating simulations that reflect both the spline predictor and its associated uncertainty . Howe ver , the covariance kernel K 1 is linked to the spectrum of the Laplace-Beltrami operator − ∆ . For a general compact and connected manifold, these eigenv alues and eigenfunctions are not av ailable in closed form, making standard Gaussian process regression computations intractable. T o overcome this di ffi culty , and as already proposed in [14], a finite- element approximation of the conditional GRF can be considered. 3. Finite element approximation of Gaussian Random Fields Here, the objectiv e is to construct a discretization of Gaussian random fields on the Rie- mannian manifold M by exploiting the finite element approximation framework, as introduced in [15]. W e discretize the manifold M by introducing a polyhedral approximation obtained by a triangulation T with m nodes c 1 , . . . , c m ∈ M . In particular , we assume (as it is standard with surface finite element methods, cf. [31]) that the resulting triangulated surface is “close enough” to the manifold M so that we can assume there is a one-to-one mapping between the points of M , and those of its discretized counterpart. As a result, we use with an abuse of notation the same letter M to denote the triangulated surface. Associated with this triangulation, we define a set of compactly supported basis functions { ψ j } m j = 1 on M . Each basis function ψ j is piecewise linear over T, equals 1 at the node c j , and vanishes at all other nodes, corresponding to a first-order finite element discretization. This construction naturally yields the following matrices: • the mass matrix M ∈ R m × m , with entries [ M ] i j = ⟨ ψ i , ψ j ⟩ L 2 ( M ) , • the sti ff ness matrix F ∈ R m × m , with entries [ F ] i j = ⟨∇ ψ i , ∇ ψ j ⟩ L 2 ( M ) , • the projection matrix P n ∈ R n × m mapping the triangulation nodes to the observation points, with entries [ P n ] i j = ψ j ( s i ). The mass matrix M is symmetric positiv e definite, and the sti ff ness matrix F is symmetric positiv e semi-definite [16]. Consequently , one can introduce a matrix square root √ M ∈ R m × m satisfying √ M √ M ⊤ = M , which may be obtained, for example, via a Cholesky factoriza- tion [32]. Owing to the local support of the basis functions { ψ j } m j = 1 , both M and F are sparse matrices. 3.1. F inite element statistical model Let us propose a finite-element approximation of the statistical model introduced in Sec- tion 2.3. First, it should be reminded that a finite-element approximation relies solely on the values at the triangulation nodes and is then e xtended over the entire domain by linear combina- tion. For instance, the random field Z modeling the phenomenon will be approximated by Z ( m ) ( s ) = m X j = 1 Z j ψ j ( s ) = P ( s ) ⊤ Z , s ∈ M (9) 6 where for 1 ≤ j ≤ m , Z j = Z ( m ) ( c j ) corresponds to the v alue of the finite-element approximation Z ( m ) at the mesh node c j , Z = ( Z 1 , . . . , Z m ) ⊤ , and P ( s ) = ( ψ 1 ( s ) , . . . , ψ m ( s ) ) ⊤ . W e therefore suppose that from now on, the observations y from our initial spline problems in Section 2.1 are noisy realizations of the random field Z ( m ) , i.e. y i is a realization of Y ( m ) i = Z ( m ) ( s i ) + στ E i , 1 ≤ i ≤ n (10) where E = ( E i ) n i = 1 is a centered Gaussian vector with Co v( E ) = I n , independent of Z . Finally , following the same approach as in Section 2.3, we assume that the distribution of Z , which entirely charecterizes the distribution of Z ( m ) arises from the following two-stage hierar- chical model: •  Z ( m ) ( s ) | A = a  =  P ( s ) ⊤ Z | A = a  is a centered Gaussian vector • A ∼ N (0 , α ) ( µ = 0 without loss of generality) The objecti ve is then to define the distribution of the random vector ( Z | A = a ) , so that  Z ( m ) | A = a  is the finite-element approximation of ( Z | A = a ) . As a reminder , ∀ a ∈ R , ( Z | A = a ) is a GRF with mean a ϕ 0 and covariance kernel σ 2 K 1 . The result is presented and justified in [33, 34], and is briefly recalled here. Let us denote − ∆ m the Galerkin approximation of − ∆ on the triangulation T, defined explic- itly in Appendix D, and  λ ( m ) k  m − 1 k = 0 and  ϕ ( m ) k  m − 1 k = 0 denote its eigen values and eigen vectors. Then, to characterize the distrib ution of ( Z | A = a ) , we need to introduce the follo wing matrix S defined as S =  √ M  − 1 F  √ M  −⊤ , (11) which is real, symmetric, and positiv e semi-definite with same eigenv alues as − ∆ ( m ) . [34]. Then, using the diagonalization S = V Diag  λ ( m ) 0 , . . . , λ ( m ) m − 1  V ⊤ , we hav e the following result. Proposition 3 (Distrib ution of the random vector ( Z | A = a ) .) . F or all a ∈ R , the random vector ( Z | A = a ) is a Gaussian vector with • mean a ϕ 0 • covariance matrix σ 2 Σ with Σ =  √ M  −⊤ f ( S )  √ M  − 1 , wher e ϕ 0 =  ϕ ( m ) 0 ( c j )  m j = 1 = 1 m ∥ √ M ⊤ 1 m ∥ 2 (see Appendix E), f is the function such that f ( λ ) = 1 λ 2 1 λ> 0 and f ( S ) is the matrix function defined fr om the eigendecomposition of S as f ( S ) = V Diag  f ( λ ( m ) 0 ) , . . . , f ( λ ( m ) m − 1 )  V ⊤ . Using proposition 3 and Equation 9, it provides the distrib ution of the finite-element random field  Z ( m ) ( s ) | A = a  . Corollary 1 (Distrib ution of the finite-element approximation of the GRF .) . F or all a ∈ R ,  Z ( m ) ( s ) | A = a  is a Gaussian random field with mean a P ( s ) ⊤ ϕ 0 and covariance kernel σ 2 K ( m ) 1 with K ( m ) 1 ( s , s ′ ) = P ( s ) ⊤ Σ P ( s ′ ) , ( s , s ′ ) ∈ M 2 . 7 In practice, and as commonly adopted in finite element implementations, the mass matrix M is replaced by its mass lumping approximation [35]. This approximation consists of a diagonal matrix whose entries are giv en by [ M ] j j = ⟨ ψ j , 1 ⟩ for 1 ≤ j ≤ m . For notational simplicity , we retain the same symbol M for the mass lumping matrix. This simplification facilitates access to  √ M  − 1 and ensures that S is sparse. Now that we have characterized the distribution of the discretized field Z ( m ) used to model the observations Y ( m ) obs = h Y ( m ) 1 , . . . , Y ( m ) n i ⊤ introduced in Equation 10, we can move on to the characterization of the posterior distribution  Z ( m ) | Y ( m ) obs = y  when the variance α of A tends to infinity , which will be used to make predictions at new locations. 3.2. Pr edictive distribution of the finite element model Follo wing once again the approach in Section 2.3, we wish to use the posterior distribution  Z ( m ) | Y ( m ) obs = y  in the case V ( A ) = α → ∞ as a predicti ve distribution to recover the observed field at any location. Noting that the vector of observ ations Y ( m ) obs can be written as Y ( m ) obs = P n Z + στ E , (12) with P n the projection matrix defined in Section 3, we can adapt Proposition 2 to our finite- element statistical model to obtain the following result. Proposition 4 (Con ver gence of the finite element conditional GRF) . Let α = V ( A ) in the hier- ar chical model defining Z ( m ) . Then,  Z ( m ) | Y ( m ) obs = y  is a GRF such that (i) F or any s ∈ M , the mean E  Z ( m ) ( s ) | Y ( m ) obs = y  satisfies lim α →∞ E  Z ( m ) ( s ) | Y ( m ) obs = y  = a ( m ) τ P ( s ) ⊤ ϕ 0 + P ( s ) ⊤ Σ P ⊤ n  K ( m ) τ  − 1 ( y − a ( m ) τ P n ϕ 0 ) wher e a ( m ) τ =  ϕ ⊤ 0 P ⊤ n  K ( m ) τ  − 1 P n ϕ 0  − 1 ϕ ⊤ 0 P ⊤ n  K ( m ) τ  − 1 y and K ( m ) τ = P n Σ ( P n ) ⊤ + τ 2 I n (ii) F or any s , s ′ ∈ M , the covariance kernel Co v  Z ( s ) , Z ( s ′ ) | Y ( m ) obs = y  satisfies lim α →∞ Cov  Z ( s ) , Z ( s ′ ) | Y ( m ) obs = y  = σ 2 h c ( m ) ( s , s ′ ) + P ( s ) ⊤ Σ P ( s ′ ) − P ( s ) ⊤ Σ P ⊤ n  K ( m ) τ  − 1 P n Σ P ( s ′ ) i . wher e c ( m ) ( s , s ′ ) = b ( s )  ϕ ⊤ 0 P ⊤ n  K ( m ) τ  − 1 P n ϕ 0  − 1 b ( s ′ ) and b ( s ) = P ( s ) ⊤ ϕ 0 − P ( s ) ⊤ Σ P ⊤ n  K ( m ) τ  − 1 P n ϕ 0 . Again, the limiting field defined by these moments r emains a GRF . The proof is identical to that of Proposition 2 giv en in Appendix C, replacing ϕ 0 ( s ) with ϕ ( m ) 0 ( s ) = P ( s ) ⊤ ϕ 0 and K 1 ( s , s ′ ) with K ( m ) 1 ( s , s ′ ) = P ( s ) ⊤ Σ P ( s ′ ). The inv erisibility of the matrix K ( m ) τ is discussed in Appendix F. Using these results, the objectiv e is now to simulate, at the triangulation points, the GRF  Z ( m ) | Y ( m ) obs = y  of Proposition 4, which corresponds to the vector 8  Z | Y ( m ) obs = y  as α → ∞ . This will be the purpose of the next section. As an extension, in Ap- pendix G, we provide a strategy to generate simulations of  Z ( m ) | Y ( m ) obs = y , A = a ( m ) τ  , that is a GRF with the same mean as the limit presented in Proposition 2, then matching the finite element approximation of the spline prediction, but whose covariance kernel does not contain the term c ( m ) ( s , s ′ ) that refers to the uncertainty on the a priori mean, as in the simple kriging framew ork. 4. Predicti ve posterior simulations 4.1. Overview of the appr oach In this section, we define Z post ∼ N ( m post ( y ) , Σ post ), the limiting predictive distribution of Proposition 4, i.e. we take • m post ( y ) = lim α →∞ E  Z | Y ( m ) obs = y  as the posterior expectation of  Z | Y ( m ) obs = y  in the limit V ( A ) = α → ∞ . • Σ post = lim α →∞ Cov  Z | Y ( m ) obs = y  as the corresponding cov ariance kernel, which in particular does not depend on y (see Proposition 4). Our objective is to generate simulations of this Gaussian vector . T wo approaches may natu- rally be considered as possible starting points. First, within this uni versal kriging frame work, simulations of Z post can be deri ved from prior simulations of ( Z | A = a ) , that has a cov ariance matrix σ 2 Σ (as established in Proposition 3). Howe ver , computing Σ requires the diagonalization of S (defined in Equation 11). Because S is an m × m matrix, where m denotes the number of triangulation nodes, this computation becomes computationally intractable for large m . T o address this, the precision matrix, defined as the in verse of the cov ariance matrix Σ , is exploited in practice. The idea is that Gaussian random vectors defined on a triangulation mesh can be viewed as GMRFs, which implies that their precision matrix is sparse. Howe ver , direct manipulation of the precision matrix is not feasible in our setting. Since rank( Σ ) = rank  f ( S )  = m − 1, with f ( λ ( m ) 0 ) = 0 and f ( λ ( m ) k ) > 0 for 1 ≤ k ≤ m − 1, the random vector ( Z | A = a ) is an intrinsic GMRF , rendering classical computational methods in volving the precision matrix inapplicable. Alternativ ely , one may approximate the α → ∞ limit by simulating  Z | Y ( m ) obs = y  for a su ffi ciently large α . These simulations leverage the in vertibility of ˜ Σ α , the cov ariance matrix of Z , defined as ˜ Σ α = σ 2 Σ + α ϕ 0 ϕ ⊤ 0 . The computations leading to this cov ariance matrix are identical to those detailed in Appendix C.1. Its in verse, the precision matrix of Z , is defined as: ˜ Q α = ˜ Σ − 1 α = σ − 2 √ MS 2 √ M ⊤ + 1 α  M ϕ 0   M ϕ 0  ⊤ . (13) The structure of ˜ Q α is well-suited for sparse computations. Specifically , one can compute the Cholesky factorization of the sparse matrix √ MS 2 √ M ⊤ and subsequently obtain the Cholesky factorization of ˜ Q α via a rank-one update. Ho wever , this strategy leads to numerical di ffi culties that depend on the value of α and relies on the approximation of a limit, which prev ents exact computation. The approach proposed here inv olves lev eraging the matrix ˜ Q α for a fixed α to deriv e exact simulations of Z post , thereby av oiding the need for a limiting approximation. This process is decomposed into the following steps: 9 1. Generate simulations of the intrinsic GMRF ( Z | A = a ) using ˜ Q α (Section 4.2). 2. Compute the posterior expectation m post ( y ) using ˜ Q α , without limit estimation (Section 4.3) 3. Construct exact simulations of Z post by combining the results from the previous two steps (Section 4.4) The ov erall procedure is summarized in Section 4.5. 4.2. Simulating an intrinsic GMRF Here, we detail the procedure to obtain simulations, for a giv en a ∈ R , of the intrinsic GMRF Z a : = ( Z | A = a ) ∼ N ( a ϕ 0 , σ 2 Σ ) , using the GMRF Z ∼ N ( 0 m , ˜ Q − 1 α ). The work of [19] can be adapted here and leads to the following result, which is v alid for all α > 0 . Proposition 5 (Simulation of intrinsic GMRF) . Considering Z ∼ N ( 0 m , ˜ Σ α ) , then  a ϕ 0 + Z | h  M ϕ 0  ⊤ Z = 0 i = a ϕ 0 + Z − ϕ 0  M ϕ 0  ⊤ Z ∼ N ( a ϕ 0 , σ 2 Σ ) . Pr oof. The vector Z − ϕ 0  M ϕ 0  ⊤ Z =  I m − ϕ 0  M ϕ 0  ⊤  Z is a centered Gaussian vector , since it is a linear transformation of the centered Gaussian vector Z . Its cov ariance is  I m − ϕ 0  M ϕ 0  ⊤  Cov( Z )  I m − ϕ 0  M ϕ 0  ⊤  ⊤ =  I m − ϕ 0  M ϕ 0  ⊤   σ 2 Σ + α ϕ 0 ϕ ⊤ 0   I m − M ϕ 0 ϕ ⊤ 0  . As shown in Appendix F, Σ M ϕ 0 = 0 m . Then,  I m − ϕ 0  M ϕ 0  ⊤   σ 2 Σ + α ϕ 0 ϕ ⊤ 0   I m − M ϕ 0 ϕ ⊤ 0  = σ 2 Σ + α ϕ 0 ϕ ⊤ 0 − 2 α ϕ 0  M ϕ 0  ⊤ ϕ 0 ϕ ⊤ 0 + α ϕ 0  M ϕ 0  ⊤ ϕ 0 ϕ ⊤ 0 M ϕ 0 = σ 2 Σ + α ϕ 0 ϕ ⊤ 0 − 2 α ϕ 0 ϕ ⊤ 0 + α ϕ 0 ϕ ⊤ 0 using ϕ ⊤ 0 M ϕ 0 = 1 see Appendix E . = σ 2 Σ Since the simulation of Z ∼ N ( 0 m , ˜ Σ α ) can be carried out by solving a linear system in volving L α the Cholesky factor of ˜ Q α , and a vector of independent standard Gaussian components, this yields the following algorithm for simulating Z a . Algorithm 1 Simulation of Z a ∼ N ( a ϕ 0 , σ 2 Σ ) Require: M , F , ϕ 0 , σ , α, a and a vector ε m of m independent standard Gaussian components 1: S ← √ M − 1 F √ M −⊤ 2: Q ← √ M S 2 √ M ⊤ 3: Compute L α the Cholesky factor ˜ Q α using sparsity of σ − 2 Q and a rank-one update 4: Solve the linear system L ⊤ α z = ε m 5: retur n z a = a ϕ 0 + z − ϕ 0  M ϕ 0  ⊤ z 10 4.3. Computing the posterior expectation The objecti ve here is to compute m post ( y ). The work of [14] establishes an analytical rela- tionship between m ( α ) post ( y ) : = E  Z | Y ( m ) obs = y  = ˜ Q − 1 α P ⊤ n h P n ˜ Q − 1 α P ⊤ n + σ 2 τ 2 I n i − 1 y , and m post ( y ) = lim α →∞ m ( α ) post ( y ) . Proposition 6 (Computation of the posterior expectation.) . The posterior expectation m post ( y ) can be computed using m ( α ) post ( y ) and m ( α ) post  P n ϕ 0  , using the r elation m post ( y ) = m ( α ) post ( y ) + ϕ 0 − h h m ( α ) post  P n ϕ 0  i  M ϕ 0  ⊤ m ( α ) post  P n ϕ 0  − ϕ 0 + h h m ( α ) post  P n ϕ 0  i !  M ϕ 0  ⊤ m ( α ) post ( y ) , wher e h [ m ] = m − ϕ 0 ( M ϕ 0 ) ⊤ m 1 − ( M ϕ 0 ) ⊤ m , ∀ m ∈ R m . Using Proposition 6 which is valid for all α > 0, the posterior expectation m post ( y ) can be computed for any y ∈ R n , , pro vided that m ( α ) post ( y ) is tractable for any y ∈ R n . For the computation of m ( α ) post ( y ), two scenario are distinguished. 1. The first scenario corresponds to the study of interpolating splines, i.e., setting the noise τ = 0. T o facilitate the computations in this situation, we require that the observation points S = ( s i ) n i = 1 are included among the triangulation nodes c 1 , . . . , c m . More precisely , let I = { i 1 , . . . , i n } denote the set of indices such that s k = c i k for 1 ≤ k ≤ n . This requirement is not restrictiv e in practice, as the triangulation can typically be designed to include the observation points, while simplifying the computations. Indeed, with this hypothesis, we hav e ( P n ) k , j =        1 if j = i k , 0 otherwise , for 1 ≤ k ≤ n and 1 ≤ j ≤ m , leading to P n ˜ Q − 1 α P ⊤ n = h ˜ Q − 1 α i I , I , by denoting h ˜ Q α i I 1 , I 2 the submatrix of ˜ Q α consisting of the rows indexed by I 1 and the columns indexed by I 2 for any inde x sets I 1 , I 2 . Then, we hav e the simplified results         m ( α ) post ( y )  ¯ I =  ˜ Q α  − 1 ¯ I , ¯ I  ˜ Q α  ¯ I , I y  m ( α ) post ( y )  I = y , and the vector m ( α ) post ( y ) can be computed with the following procedure. 11 Algorithm 2 Computation of m ( α ) post ( y ) in the first scenario Require: M , F , ϕ 0 , I , σ , α, y 1: S ← √ M − 1 F √ M −⊤ 2: Q ← √ M S 2 √ M ⊤ 3: ϕ 0 ← 1 m ∥ √ M ⊤ 1 m ∥ 2 4: Compute the Cholesky decomposition of h ˜ Q α i ¯ I , ¯ I = σ − 2 [ Q ] ¯ I , ¯ I + 1 α  M ϕ 0  ¯ I  M ϕ 0  ⊤ ¯ I using sparsity of σ − 2 [ Q ] ¯ I , ¯ I and a rank-one update 5: Solve h ˜ Q α i ¯ I , ¯ I ω = h ˜ Q α i ¯ I , I y using sparse linear systems 6: retur n       [ m ( α ) post ( y )] ¯ I [ m ( α ) post ( y )] I       = " ω y # 2. The second scenario corresponds to the computation of smoothing splines, i.e., setting the noise τ > 0. In this case, the assumption that the observation points belong to the triangulation nodes is no longer required to obtain computational simplifications. Indeed, one has ˜ Q − 1 α P ⊤ n  P n ˜ Q − 1 α P ⊤ n + σ 2 τ 2 I n  − 1 y =  σ 2 τ 2 ˜ Q α + P ⊤ n P n  − 1 P ⊤ n y by the W oodbury formula. Using the fact that P ⊤ n P n is sparse by construction of the basis functions ( ψ j ) m j = 1 , the decomposition of σ 2 τ 2 ˜ Q α + P ⊤ n P n can be e ffi ciently obtained by expressing this matrix as the sum of a sparse matrix and a rank-one update. Algorithm 3 then provides the computation of m ( α ) post ( y ) for τ > 0. Algorithm 3 Computation of m ( α ) post ( y ) in the second scenario Require: M , F , ϕ 0 , α, P n , y , τ > 0 , σ > 0 1: S ← √ M − 1 F √ M −⊤ 2: Q ← √ M S 2 √ M ⊤ 3: ϕ 0 ← 1 m ∥ √ M ⊤ 1 m ∥ 2 4: Compute the Cholesk y of σ 2 τ 2 ˜ Q α + P ⊤ n P n = τ 2 Q + P ⊤ n P n + σ 2 τ 2 α  M ϕ 0   M ϕ 0  ⊤ using sparsity of τ 2 Q + P ⊤ n P n and a rank-one update 5: Solve  σ 2 τ 2 ˜ Q α + P ⊤ n P n  ω = P ⊤ n y using sparse linear systems 6: retur n m ( α ) post ( y ) = ω 4.4. Generating the conditional simulations Our approach exploits the prior simulations of ( Z a = Z | A = a ) (Section 4.2) and the calcu- lation of the posterior expectation (Section 4.3) to obtain simulations of Z post . T o this end, we use the following result, which specifically addresses the simulation of residuals [36]. Proposition 7 (Simulation of the residuals.) . F or all a ∈ R , let Z a ∼ N ( a ϕ 0 , σ 2 Σ ) , and E ∼ N (0 , I n ) . Then, the vector of r esiduals Z r es = m post ( P n Z a + στ E ) − Z a is a center ed Gaussian vector with covariance matrix Σ post . 12 Pr oof. The objectiv e is to compute          E ( Z res ) = E  m post ( P n Z a + στ E )  − a ϕ 0 Cov( Z res ) = Cov h m post ( P n Z a + στ E ) i − 2Cov h m post ( P n Z a + στ E ) , Z a i + Cov [ Z a ] (14) using Cov ( P n Z a + στ E ) = σ 2 K ( m ) τ Cov [ Z a ] = σ 2 Σ . Step 1: Compute m post ( P n Z a + στ E ) . Using Proposition 4, hav e m post ( P n Z a + στ E ) = ˜ a ( m ) τ ϕ 0 + Σ P ⊤ n  K ( m ) τ  − 1 ( P n Z a + στ E − ˜ a ( m ) τ P n ϕ 0 ) with ˜ a ( m ) τ =  ϕ ⊤ 0 P ⊤ n  K ( m ) τ  − 1 P n ϕ 0  − 1 ϕ ⊤ 0 P ⊤ n  K ( m ) τ  − 1 ( P n Z a + στ E ) . Let us denote b = ϕ 0 − Σ P ⊤ n  K ( m ) τ  − 1 P n ϕ 0 and β = ϕ ⊤ 0 P ⊤ n  K ( m ) τ  − 1 P n ϕ 0 . Then, m post ( P n Z a + στ E ) = b ˜ a ( m ) τ + Σ P ⊤ n  K ( m ) τ  − 1 [ P n Z a + στ E ] = β − 1 b ϕ ⊤ 0 P ⊤ n  K ( m ) τ  − 1 [ P n Z a + στ E ] + Σ P ⊤ n  K ( m ) τ  − 1 [ P n Z a + στ E ] = h β − 1 b ϕ ⊤ 0 P ⊤ n + Σ P ⊤ n i  K ( m ) τ  − 1 [ P n Z a + στ E ] . (15) Step 2: Compute E h m post ( P n Z a + στ E ) i and Co v h m post ( P n Z a + στ E ) i . From step 1, we ha ve E h m post ( P n Z a + στ E ) i = a h β − 1 b ϕ ⊤ 0 P ⊤ n + Σ P ⊤ n i  K ( m ) τ  − 1 P n ϕ 0 = a  b + Σ P ⊤ n  K ( m ) τ  − 1 P n ϕ 0  = a ϕ 0 from the definition of b . For the co variance, we obtain Cov h m post ( P n Z a + στ E ) i = σ 2 h β − 1 b ϕ ⊤ 0 P ⊤ n + Σ P ⊤ n i  K ( m ) τ  − 1 h β − 1 b ϕ ⊤ 0 P ⊤ n + Σ P ⊤ n i ⊤ = σ 2  β − 1 b ⊤ b + h 2 β − 1 b ϕ ⊤ 0 P ⊤ n + Σ P ⊤ n i  K ( m ) τ  − 1 P n Σ  . Step 3: Compute Cov h m post ( P n Z a + στ E ) , Z a i . Using Equation 15, we hav e Cov h m post ( P n Z a + στ E ) , Z a i = h β − 1 b ϕ ⊤ 0 P ⊤ n + Σ P ⊤ n i  K ( m ) τ  − 1 Cov [ P n Z a + στ E , Z a ] = σ 2 h β − 1 b ϕ ⊤ 0 P ⊤ n + Σ P ⊤ n i  K ( m ) τ  − 1 P n Σ 13 Step 4: Gathering the results. It comes E ( Z res ) = 0, and using Cov( Z a ) = σ 2 Σ , it leads to Cov( Z res ) = σ 2  b ⊤ β − 1 b − Σ P ⊤ n  K ( m ) τ  − 1 P n Σ + Σ  = lim α →∞ Cov  Z | Y ( m ) obs = y  using Proposition 4 = Σ post Then, the strategy adopted here to obtain a simulation of Z post ∼ N ( m post ( y ) , Σ post ) is to simulate Z a ∼ N ( a ϕ 0 , σ Σ ) for a given a and then return Z post = m post ( y ) + m post ( P n Z a + στ E ) − Z a . (16) As the result does not depend on a , we work with a = 0, without loss of generality . It leads to the following o verall procedure, detailed in Section 4.5. 4.5. Overall simulation algorithm Follo wing the procedure outlined in the previous sections, Algorithm 4 presents the method for generating spline simulations, with or without observation noise. Algorithm 4 Generation of splines simulations Require: M , F , ϕ 0 , α, τ , σ , y , a vector ε m of m independent standard Gaussian components 1: Generate a simulation z 0 , centered with cov ariance matrix σ 2 Σ , using Algorithm 1 2: if τ = 0 then 3: Require : I 4: Compute m ( α ) post ( y ), m ( α ) post ( [ z 0 ] I ) and m ( α ) post  P n ϕ 0  using Algorithm 2 5: Compute m post ( y ) and m post ( [ z 0 ] I ) using Proposition 6 6: else 7: Require : P n and ε n a vector of n independent standard Gaussian components 8: Compute m ( α ) post ( y ), m ( α ) post ( P n z 0 + στ ε n ) and m ( α ) post  P n ϕ 0  using Algorithm 3 9: Compute m post ( y ) and m post ( P n z 0 + στ ε n ) using Proposition 6 10: end if 11: retur n z post = m post ( y ) + m post ( P n z 0 ) − z 0 Note that these results are theoretically valid for any α > 0. Ho wever , in practice, this parameter cannot be chosen too large or too small in order to av oid numerical instabilities. As detailed in [14] and adapted here to account for the additional parameter σ , it is important to ensure that σ √ α ∈ [ λ ( m ) 1 , λ ( m ) m − 1 ] , to control the condition number of the matrices in volved in the Cholesky decomposition. The power -iteration method can be used to estimate these eigen values [37]. 14 5. Anisotropy and parameter estimation 5.1. Anisotr opic splines These spline-based simulations are particularly valuable as they quantify uncertainty along- side the predicted field. Ho wev er, the cov ariance kernels introduced above, K 1 and its finite- element approximation K ( m ) 1 , are isotropic, implying that the modeled phenomenon must exhibit identical behavior in all spatial directions. A natural way to incorporate directional anisotropy , while preserving the exact modeling frame work presented here, is to modify the parametrization of the manifold ( M , g ) via its Riemannian metric g . This strategy is described in detail in [33] and is briefly summarized below . Considering that M is embedded in R d , for any point s ∈ M we consider a coordinate chart ( U , x ) such that s ∈ U , where U is an open subset of M and x : U → R d is a homeomorphism onto an open subset of R d . The Riemannian metric expressed in local coordinates is gi ven by g s ( u s , v s ) =  u x s  ⊤ G x ( s ) v x s , (17) where u x s and v x s denote the coordinate representations of the tangent vectors u s , v s ∈ T s M in the chart ( U , x ), and G x ( s ) is the matrix representation of the Riemannian metric at s in these local coordinates. The matrix G x ( s ) is symmetric positiv e definite and admits the decomposition G x ( s ) = R x ( s ) D x ( s ) 2 ( R x ( s ) ) ⊤ , which can be interpreted as a local deformation operator . For d ∈ { 2 , 3 } , the matrix R x ( s ) cor- responds to a rotation, while D x ( s ) is a diagonal matrix containing scaling factors. Under this interpretation, the Riemannian metric induces a local deformation of the manifold by first ap- plying a rotation and then an anisotropic scaling that emphasizes preferred directions. These rotations and scalings are parametrized and can be specified locally on each element of the tri- angulation. The corresponding parameters are collected in a vector β , which modifies the inner product on L 2 ( M ) and, consequently , the matrices M and S . 5.2. P arameter estimation W e then hav e unknown parameters β defining the anisotropies, and the variance σ 2 . A stan- dard approach is to rely on the universal kriging statistical model for the observations and to maximize the associated log-likelihood. It corresponds to the log-density of the random vector  Y ( m ) obs | A = a  , which is maximized with respect to  a , σ 2 , β  L ( a , σ , β ) = − 1 2  n log(2 π ) + log | σ 2 K ( m ) τ | + σ − 2  y − a P n ϕ 0  ⊤  K ( m ) τ  − 1  y − a P n ϕ 0   , where β impacts both K and ϕ 0 . The optimization leads to              a ⋆ = a ( m ) τ σ ⋆ = " 1 n  y − a P n ϕ 0  ⊤  K ( m ) τ  − 1  y − a P n ϕ 0  # 1 / 2 (18) and we end up maximizing in β the concentrated log-likelihood [38] L ( C ) ( β ) = − n log  σ ⋆  − 1 2 log | K ( m ) τ | . The computation of σ ⋆ and | K ( m ) τ | is detailed in [14]. Note that, in the smoothing spline setting, the noise variance τ > 0 can be estimated jointly with β . 15 6. Application test cases Here, we present two application test cases in volving the prediction of smooth phenomena from a limited number of observations. The first application concerns temperature prediction on the Earth, modeled as a spherical surface. The second considers an analytical function defined on a cylindrical surface. In both cases, anisotropies are estimated in local coordinate charts, namely spherical coordinates for the sphere and cylindrical coordinates for the cylinder . For both test cases, the locations of the observation points are fully controlled, so we consider the first scenario and focus on the interpolation problem. The results obtained in the second scenario are similar and are provided in the supplementary material [39]. 6.1. Pr ediction on Earth surface Here, we consider the daily analysed Sea Surface T emperature (SST) field for December 26 th , 2023, obtained from the Copernicus satellite SST dataset [40]. It corresponds to a gap-free, gridded estimate of sea surface temperature at approximately 20 cm depth, obtained through optimal interpolation combining observations from multiple satellite sensors. The SST field ev olves smoothly over the oceans and exhibits strong anisotropy , as latitude is expected to have a major influence on temperature values. W e collect the data on a latitude-longitude grid with a 1 . 5 ◦ × 1 . 5 ◦ resolution, which defines the triangulation mesh ( c j ) n j = 1 . Obviously , the data points are restricted to oceanic regions only . For the prediction study , we use only n = 20 observ ation points y = ( y ( s i ) ) n i = 1 , in order to illustrate the e ff ectiveness of splines for reconstructing a smooth phenomenon from a limited number of observations. The observation locations are selected using a Maximin Latin Hypercube Sampling design (LHS, [41]). In the first scenario, the observation points must belong to the mesh nodes ( c j ) n j = 1 . The final design is therefore obtained by selecting the mesh nodes closest to the points generated by the LHS. The remaining mesh nodes, for which the data values are known, are then used as a validation dataset. For the anisotropic model, two parameters were considered: the rotation angle in spherical coordinates and the ratio between the two scaling parameters in this coordinate chart. These parameters were estimated by maximum likelihood, as described in Section 5.2. Figure 1: Predictions of the SST ( ◦ C ) on Earth, illustrated in 2D using the spherical coordinates ( θ, ϕ ) . n = 20 observ ation points are shown as black dots. Left: true values. Middle: isotropic splines. Right : Anisotropic splines. Figure 1 displays the prediction mean obtained on the sphere ov er the oceanic regions using the isotropic and the anisotropic spline. It clearly highlights the influence of anisotropy along 16 (a) Distribution of the prediction errors and associated RMSE in the isotropic and the anisotropic case. (b) Distribution of the scores for the isotropic and the anisotropic prediction. Figure 2: Comparison of the prediction errors and scores for the isotropic and anisotropic SST predictions on the Earth. the longitudinal direction, which substantially improves the prediction, as illustrated by the dis- tribution of prediction errors in Figure 2a and the associated Root Mean Square Error (RMSE). Figure 3 displays the uncertainty associated with the prediction, obtained from 500 simulations, along the line of longitude closest to ϕ = 170 ◦ in the triangulation mesh, crossing the Pacific Ocean. Once again, the anisotropic splines better capture the behavior of the phenomenon. Howe ver , this figure also emphasizes the importance of properly assessing uncertainty . In the isotropic case, although the prediction mean is less accurate, the associated uncertainty remains consistent. For example, when computing the 95% prediction interv al over all points of the tri- angulation mesh, the true SST value lies within this interval for more than 98% of the points in the isotropic setting. Howe ver , such coverage is not a relev ant performance criterion, since an arbitrarily large prediction variance would ensure that all true v alues fall within a giv en confidence interval. As discussed in [42], a rele vant metric that promotes a proper balance between predictive accuracy and uncertainty quantification is obtained by considering, at a prediction point s , the score S ( s ) = − ˆ y ( s ) − y true ( s ) ˆ σ 2 ( s ) ! 2 − log  ˆ σ 2 ( s )  , (19) which is directly related to the log-likelihood of the Gaussian predictive distribution at s , with predictiv e mean ˆ y ( s ) and predictive variance ˆ σ 2 ( s ). This score is computed at all points of the triangulation mesh, and its distrib ution is sho wn in Figure 2b, with higher scores indicating better predictiv e performance. The results once again confirm the benefit of incorporating anisotropy in the spline model for this phenomenon. Note that the case of the sphere is particularly interesting, since the eigenfunctions of − ∆ are the spherical harmonics [43]. As detailed in [14], the kernel K 1 ( s 1 , s 2 ) = X k ∈ N ⋆ 1 λ 2 k ϕ k ( s 1 ) ϕ k ( s 2 ) can therefore be computed explicitly by truncating the series expansion. In practice, the sum is truncated at K = 40, following the recommendations of [21]. Classical kriging operations can then be performed using this kernel. The resulting predictions serve for comparison with 17 (a) Predictiv e distribution in the isotropic case. (b) Predictiv e distribution in the anisotropic case. Figure 3: Comparison between the uncertainty quantification in the isotropic and the anisotropic SST prediction on the Earth, along the longitude closest to 170 ◦ . the finite element approximation, in order to validate our numerical approach, and are provided in Appendix H. 6.2. Pr ediction on the surface of a cylinder Here, we consider a cylinder of radius 1 and height 20, since in many industrial applications the height of a cylinder is typically much larger than its radius (e.g., in the study of subsea pipelines [44]). W e define C as the surface under in vestigation: C = n ( x , y , z ) ∈ R 2 × [0 , 20] | x 2 + y 2 = 1 o . As a case study , we consider the three-dimensional Franke function [45], which is widely used for interpolating smooth functions [46, 47]. It is defined on [0 , 1] 3 by g ( x , y , z ) = 0 . 75 exp − ( a x (9 x − 2)) 2 + ( a y (9 y − 2)) 2 + ( a z (9 z − 2)) 2 4 ! + 0 . 75 exp − ( a x (9 x + 1)) 2 49 + ( a y (9 y + 1)) 2 10 + ( a z (9 z + 1)) 2 10 !! + 0 . 5 exp − ( a x (9 x − 7)) 2 + ( a y (9 y − 3)) 2 + ( a z (9 z − 5)) 2 4 ! − 0 . 2 exp  −  ( a x (9 x − 4)) 2 + ( a y (9 y − 7)) 2 + ( a z (9 z − 5)) 2  , where a x , a y , and a z are coe ffi cients that introduce anisotropy , set to a x = 0 . 4, a y = 0 . 4, and a z = 1. The phenomenon on C is then modeled by the function ˜ g defined as ∀ ( x , y , z ) ∈ C , ˜ g ( x , y , z ) = g x + 1 2 , y + 1 2 , z 20 ! , so that g is applied on [0 , 1] 3 . 18 (a) 3D plots of the predictions compared to the true function. The cylinder’ s height is normalized for visualization purposes. (b) 2D plots of the predictions compared to the true function. The x -axis represents the longitudinal direction, while the y -axis represents the polar angle in degrees. Figure 4: Predictions on the cylinder . Left: true function. Middle: isotropic spline. Right: anisotropic spline. n = 20 observation points are sho wn as black dots. A triangulation mesh is constructed using a regular grid in the cylindrical coordinate chart ( θ, z ), with a 5 ◦ resolution in θ and a resolution of 0 . 5 in z . A total of n = 10 observation points are selected as described in Section 6.1. The anisotropy parameters are estimated using the same approach as before, i.e. with the rotation angle and the ratio between the scaling parameters in the cylindrical coordinate chart. The resulting prediction mean is displayed in Figure 4, while the prediction errors are shown in Figure 5a. Once again, estimating anisotropy substantially improves predictive accuracy . The uncer- tainty quantification appears consistent with the true values, as illustrated in Figure 6, which presents the 95% prediction interval along the line θ = 270 ◦ . Finally , the distribution of the score defined in Equation 19, which balances predicti ve accuracy and uncertainty quantifica- tion, further confirms the improved performance obtained when incorporating anisotropy (see Figure 5b). 7. Summary and perspectives This article extends the work of [14] by proposing a framework for spline interpolation on compact Riemannian manifolds, together with a quantification of the associated uncertainty . The approach relies on the equiv alence between the classical solution to the interpolation problem and the uni versal kriging predictor . Howe ver , it is based on the spectrum of the Laplace-Beltrami operator − ∆ , which is generally unknown on compact Riemannian manifolds. 19 (a) Distribution of the prediction errors and associated RMSE in the isotropic and the anisotropic case. (b) Distribution of the scores for the isotropic and the anisotropic prediction. Figure 5: Comparison of the prediction errors and scores for the isotropic and anisotropic predictions on the cylinder . (a) Predictiv e distribution in the isotropic case. (b) Predictiv e distribution in the anisotropic case. Figure 6: Comparison between the uncertainty quantification in the isotropic and the anisotropic prediction on the c ylin- der , along the line θ = 270 ◦ . 20 T o overcome this di ffi culty , a finite-element approximation of the predictor is introduced, based on a triangulation of the manifold as proposed in [33]. This strategy avoids the need for the explicit spectrum of − ∆ and enables the incorporation of anisotropies in the covariance kernel through a parametrization of the Riemannian metric associated with the manifold. Howe ver , it leads to an intrinsic Gaussian Markov random field (GMRF) with a singular covariance matrix, which complicates both computations and simulations. Adapting the methodology of [20], we represent this intrinsic GMRF through linear condi- tioning of a non-intrinsic GMRF . The approach is illustrated on two test cases: the Earth, mod- eled as a spherical surface, and a cylinder . For both examples, the impact of anisotropy and the role of uncertainty quantification are emphasized. The predictiv e distribution is displayed o ver a selected subdomain to demonstrate that the predictiv e variance is consistent with the actual pre- diction error , even in situations where isotropic splines yield inaccurate point estimates. Finally , a scoring rule is e valuated to assess the overall predicti ve performance, taking into account both accuracy and uncertainty quantification. Although this work o ff ers a thorough in vestigation, sev eral perspectiv es for future research remain. First, in this work we only considered constant anisotropy parameters. Howe ver , for certain phenomena, one may need to estimate a spatially varying anisotropy field. Such an ex- tension would naturally lead to a very large number of parameters to infer through maximum log-likelihood estimation. Alternative strategies would therefore be worth exploring, for in- stance by modeling the anisotropy parameters themselves as a random field within a hierarchical framew ork. This would raise significant challenges, including the specification of a suitable prior statistical model for the anisotropy field and the implementation of a fully Bayesian approach. Another promising direction is to inv estigate the link between spline prediction and the Matérn stochastic partial di ff erential equation (SPDE), defined as follows [48]: ( κ 2 − ∆ ) α/ 2 Z = W , (20) where ( κ 2 − ∆ ) α/ 2 is a pseudo-di ff erential operator and W denotes Gaussian white noise. As extensi vely studied in [18] in R d , the solution to this SPDE is an intrinsic Gaussian field when κ = 0. On a compact manifold M , when α = 2, this solution is strongly related to our spline predictor as it admits a cov ariance kernel K 1 , and a then detailed in vestigation of this result would be of particular interest. More generally , this perspective opens the way to extending to Riemannian manifolds the work of [18] on intrinsic Whittle–Matérn fields, and to comparing our approach with theirs, which in volv es the computation of the Moore–Penrose pseudoinv erse. It may then be interesting to study the behavior of our spline predictor, and more generally of intrinsic fields, in e xtrapolation scenarios, i.e. when predictions are made far from the observ ation points, and to compare it with results obtained using classical Matérn solutions. Finally , in our work we used the score defined in Equation 19 to assess the predictive perfor- mance at each point, while accounting for the associated uncertainty . This metric provides a con- venient balance between accurac y and uncertainty quantification, and remains straightforward to compute when predictive simulations are available. Howe ver , it only relies on the marginal pre- dictiv e variance at each location and does not incorporate the full covariance structure of the predictiv e distribution. More sophisticated scoring rules could therefore be considered, although their implementation becomes challenging when dealing with an intrinsic GMRF and a large number of prediction points (e.g., the entire triangulation mesh). In such situations, estimating the predictiv e cov ariance matrix empirically would require a very large number of independent 21 simulations to obtain a su ffi ciently accurate estimator, while its analytical computation is delicate due to the manipulation of singular matrices. Acknowledgements This work was fully supported by the Chaire Geolearning funded by Andra, BNP-Paribas, CCR and SCOR Foundation. Declaration of Competing Interest The authors declare that they hav e no known competing financial interests or personal rela- tionships that could hav e appeared to influence the work reported in this paper . Data A vailability The data and code used to generate the results presented in this article are provided in the Supplementary Material [39]. References [1] G. W ahba, Y . W ang, Spline function: Overvie w , Wile y StatsRef: Statistics Reference Online (1990) 1–19. [2] T . J. Santner , B. J. W illiams, W . I. Notz, The design and analysis of computer experiments, volume 1, Springer , New Y ork, 2003. [3] G. Matheron, Splines and kriging: their formal equivalence, Down to Earth Statistics: Solutions looking for geological problems (1981) 77–95. [4] O. Dubrule, Comparing splines and kriging, Computers and Geosciences 10 (1984) 327– 338. [5] D. E. Myers, Interpolation with positive definite functions, Sciences de la T erre 28 (1988) 251–265. [6] D. E. Myers, Kriging, cokriging, radial basis functions and the role of positiv e definiteness, Computers & Mathematics with Applications 24 (1992) 139–148. [7] C. W illiams, C. Rasmussen, Gaussian Processes for Machine Learning, volume 2, MIT Press, Cambridge, MA, USA, 2006. [8] L. Hewing, E. Arcari, L. P . Fröhlich, M. N. Zeilinger , On simulation and trajectory pre- diction with gaussian process dynamics, in: Learning for Dynamics and Control, PMLR, 2020, pp. 424–434. [9] G. Fournier , S. Pellerin, L. T . Phuoc, Contrôle par rotation ou par aspiration de l’écoulement autour d’un cylindre calculé par simulation des grandes échelles, Comptes rendus. Mé- canique 333 (2005) 273–278. 22 [10] W . Gao, D. Nelias, Z. Liu, Y . L yu, Numerical inv estigation of flow around one finite circular cylinder with two free ends, Ocean Engineering 156 (2018) 373–380. [11] E. Lila, J. A. Aston, L. M. Sangalli, Smooth principal component analysis over two- dimensional manifolds with an application to neuroimaging (2016). [12] P . Kim, Splines on riemannian manifolds and a proof of a conjecture by wahba, preprint (http: // www . uoguelph. ca / pkim / research. html) (2001). [13] G. W ahba, Spline interpolation and smoothing on the sphere, SIAM Journal on Scientific and Statistical Computing 2 (1981) 5–16. [14] C. Sire, M. Pereira, T . Romary , Spline interpolation on compact riemannian manifolds, arXiv preprint arXi v:2510.11239 (2025). [15] F . Lindgren, H. Rue, J. Lindström, An explicit link between gaussian fields and gaussian markov random fields: the stochastic partial di ff erential equation approach, Journal of the Royal Statistical Society: Series B (Statistical Methodology) 73 (2011) 423–498. [16] M. Pereira, Generalized random fields on Riemannian manifolds : the- ory and practice, Theses, Univ ersité Paris sciences et lettres, 2019. URL: https://pastel.hal.science/tel-02499376 . [17] H. Rue, L. Held, Gaussian Marko v random fields: theory and applications, Chapman and Hall / CRC, Boca Raton, FL, USA, 2005. [18] D. Bolin, P . Braunsteins, S. Engelke, R. Huser, Intrinsic whittle–matern fields and sparse spatial extremes, arXiv preprint arXi v:2512.23395 (2025). [19] D. P . Simpson, I. W . T urner , A. N. Pettitt, Sampling from gaussian marko v random fields conditioned on linear constraints, The Proceedings of ANZIAM 48 (2006) C1041–C1053. [20] D. Bolin, J. W allin, E ffi cient methods for gaussian markov random fields under sparse linear constraints, Advances in Neural Information Processing Systems 34 (2021) 9882–9894. [21] W . Keller , A. Borko wski, Thin plate spline interpolation, Journal of Geodesy 93 (2019) 1251–1269. [22] T . Duchamp, W . Stuetzle, Spline smoothing on surfaces, Journal of Computational and Graphical Statistics 12 (2003) 354–381. [23] M. Craioveanu, M. Puta, T . M. Rassias, Spectral Properties of the Laplace-Beltrami Oper- ator and Applications, Springer Netherlands, Dordrecht, 2001, pp. 119–211. [24] H. Urakawa, Geometry of laplace-beltrami operator on a complete riemannian manifold, Progress in di ff erential geometry 22 (1993) 347–406. [25] G. W ahba, Spline models for observational data, SIAM, Philadelphia, P A, USA, 1990. [26] G. Matheron, Le krigeage univ ersel, v olume 1, École nationale supérieure des mines de Paris, P aris, 1969. 23 [27] T . Hengl, A practical guide to geostatistical mapping, volume 52, Hengl Amsterdam, Netherlands, 2009. [28] H. W ackernagel, Ordinary kriging, in: Multi variate geostatistics: an introduction with applications, Springer , New Y ork, 2003, pp. 79–88. [29] P . J. Diggle, P . J. Ribeiro Jr , Bayesian inference in gaussian model-based geostatistics, Geographical and en vironmental modelling 6 (2002) 129–146. [30] C. Helbert, D. Dupuy , L. Carraro, Assessment of uncertainty in computer experiments from uni versal to bayesian kriging, Applied Stochastic Models in Business and Industry 25 (2009) 99–113. [31] G. Dziuk, C. M. Elliott, Finite element methods for surface pdes, Acta Numerica 22 (2013) 289–396. [32] C. Benoit, Note sur une méthode de résolution des équations normales provenant de l’application de la méthode des moindres carrés à un système d’équations linéaires en nom- bre inférieur à celui des inconnues (procédé du commandant cholesky), Bulletin géodésique 2 (1924) 67–77. [33] M. Pereira, N. Desassis, D. Allard, Geostatistics for large datasets on riemannian manifolds: A matrix-free approach, Journal of Data Science 20 (2022) 512–532. doi:10.6339 / 22- JDS1075. [34] A. Lang, M. Pereira, Galerkin–chebyshev approximation of gaussian random fields on compact riemannian manifolds, BIT Numerical Mathematics 63 (2023) 51. [35] A. Quarteroni, Modellistica numerica per problemi di ff erenziali, volume 97, Springer , Cham, Switzerland, 2016. [36] J.-P . Chiles, P . Delfiner , Conditional Simulations, John Wile y and Sons, Ltd, Chichester , 2012, pp. 478–628. doi:https: // doi.org / 10.1002 / 9781118136188.ch7. [37] R. Mises, H. Pollaczek-Geiringer , Praktische verf ahren der gleichungsauflösung., ZAMM- Journal of Applied Mathematics and Mechanics / Zeitschrift für Angew andte Mathematik und Mechanik 9 (1929) 58–77. [38] D. R. Jones, A taxonomy of global optimization methods based on response surfaces, Journal of global optimization 21 (2001) 345–383. [39] C. Sire, M. Pereira, Supplementary material : Uncertainty quantifica- tion of spline predictors on compact riemannian manifolds, 2026. URL: https://doi.org/10.5281/zenodo.19186963 . [40] Copernicus Climate Change Service (C3S), Global sea surface temperature deriv ed from satellite observations, 2025. doi:10.24381 / cds.cf608234, accessed on DD-MMM-YYYY . [41] Q. Y . Kenn y , W . Li, A. Sudjianto, Algorithmic construction of optimal symmetric latin hypercube designs, Journal of statistical planning and inference 90 (2000) 145–159. 24 [42] T . Gneiting, A. E. Raftery , Strictly proper scoring rules, prediction, and estimation, Journal of the American statistical Association 102 (2007) 359–378. [43] F . Dai, Y . Xu, F . Dai, Y . Xu, Spherical harmonics, Approximation theory and harmonic analysis on spheres and balls (2013) 1–27. [44] D. Seth, B. Manna, J. T . Shahu, T . Fazeres-Ferradosa, F . T . Pinto, P . J. Rosa-Santos, Buck- ling mechanism of o ff shore pipelines: a state of the art, Journal of Marine Science and Engineering 9 (2021) 1074. [45] R. Frank e, A critical comparison of some methods for interpolation of scattered data, T ech- nical Report, 1979. [46] M. Challacombe, A n-body solver for free mesh interpolation, arXi v preprint arXiv:1511.01353 (2015). [47] R. Franke, K. Šalkauskas, Localization of multi variate interpolation and smoothing meth- ods, Journal of computational and applied mathematics 73 (1996) 79–94. [48] P . Whittle, On stationary processes in the plane, Biometrika (1954) 434–449. Appendix A. Spectral Theorem Proposition 8 (Spectrum of − ∆ ) . The Spectral Theor em pro vides the following pr operties [23]. (i) F or the closed, the Dirichlet, and the Neumann eigen value pr oblems, the spectrum of − ∆ is an infinite, countable sequence of r eal eigen values 0 ≤ λ 0 ≤ λ 1 ≤ · · · ≤ λ k ≤ . . . , wher e each eigen value λ k appears with finite multiplicity m k . (ii) Ther e exists an orthonormal basis ( ϕ k ) k ∈ N of L 2 ( M ) such that, for all k ∈ N , ϕ k ∈ C ∞ ( M ) is an eigenfunction associated with λ k . Consequently , for any u ∈ L 2 ( M ) , u = X k ∈ N u k ϕ k , u k = ⟨ u , ϕ k ⟩ L 2 ( M ) = Z M u ϕ k d µ g . (iii) F or the closed and the Neumann eigen value pr oblems, the eigen value λ 0 = 0 corr esponds to the constant functions [24], which implies that m 0 = 1 and ϕ 0 = 1 R M d µ g . (iv) F or the Dirichlet eigen value problem, all eig envalues ar e strictly positive. In this study , we consider λ 0 = 0, as in the closed and Neumann eigenv alue problems, noting that the case λ 0 > 0 (Dirichlet problem) is computationally simpler for the operations that will be performed. 25 Appendix B. In versibility of K 1 ,τ The in versibility of K 1 ,τ is guaranteed when τ > 0 . In the following, we focus on the case τ = 0 , referring to the interpolating problem. For k ∈ N ⋆ , let us denote w k = ( ϕ k ( s 1 ) , . . . , ϕ k ( s n ) ) ⊤ . ∀ x ∈ R n , x ⊤ K 1 , 0 x = n X i = 1 n X j = 1 x i x j X k > 0 1 λ 2 k ϕ k ( s i ) ϕ k ( s j ) = X k > 0 1 λ 2 k n X i = 1 n X j = 1 x i x j ϕ k ( s i ) ϕ k ( s j ) = X k > 0 1 λ 2 k        n X i = 1 x i ϕ k ( s i )        2 Then, x ⊤ K 1 , 0 x = 0 ⇔ ∀ k > 0 , n X i = 1 x i ϕ k ( s i ) = 0 ⇔ ∀ k > 0 , ⟨ x , w k ⟩ R n = 0 It comes K 1 , 0 is positiv e definite ⇔ span { w k } k > 0 = R n ⇔ ∀ x ∈ R n , ∃ ( a k ) k > 0 , x = X k > 0 a k w k ⇔ ∀ x ∈ R n , ∃ ( a k ) k > 0 , ∀ 1 ≤ i ≤ n , x i = X k > 0 a k ϕ k ( s i ) , ⇔ ∀ x ∈ R n , ∃ ( a k ) k > 0 , the function u = X k > 0 a k ϕ k interpolates the vector x As ( ϕ k ) k > 0 is a basis of L 2 ( M ) , K 1 , 0 is in vertible if and only if the an interpolating function exists in L 2 ( M ) for any v ector of observations x . Note that this condition cannot be satisfied when an observation point lies on the boundary ∂ M under Dirichlet boundary conditions. Appendix C. Proof of Proposition 2 Our random vector of observ ations Y obs = ( Y i ) n i = 1 verifies Y i = Z ( s i ) + στ E i , 1 ≤ i ≤ n , with • ( Z | A = a ) a GRF with mean a ϕ 0 and cov ariance kernel σ 2 K 1 , for all a ∈ R 26 • A ∼ N ( µ, α ) • E = ( E i ) n i = 1 a centered Gaussian vector with Co v( E ) = I n , independent of A and Z . The objective is to compute the mean and cov ariance kernel of ( Z ( s ) | Y obs = y ) when α tends to + ∞ . Appendix C.1. Distribution of Z ( s ) Here, we will show that Z ( s ) is a GRF with • mean µϕ 0 • cov ariance kernel αϕ 0 ( s ) ϕ 0 ( s ′ ) + Σ 2 K 1 ( s , s ′ ) Note that this step is not mandatory , but is included here since a very similar dev elopment can be used to show that Z ∼ N ( 0 m , ˜ Σ α ) in Section 4.1. Let us consider x 1 , . . . , x p , p points of M , and study the distribution of Z test = Z ( x 1 ) , . . . , Z ( x p ) . W e denote ϕ test 0 = ( ϕ 0 ( x i ) ) p i = 1 , and K test 1 = h K 1 ( x i , x j ) i 1 ≤ i , j ≤ p . Then density of Z test is p Z test ( z ) = Z R p p Z test | A = a ( z ) p A ( a ) d a = Z R 1 (2 π ) ( p + 1) / 2 | ( σ 2 K test 1 ) | 1 / 2 exp − 1 2 ( z − a ϕ test 0 ) ⊤ ( σ 2 K test 1 ) − 1 ( z − a ϕ test 0 ) ! 1 √ α exp − ( a − µ ) 2 2 α ! d a = 1 √ α (2 π ) ( p + 1) / 2 | ( σ 2 K test 1 ) | 1 / 2 exp − 1 2 z ⊤ ( σ 2 K test 1 ) − 1 z + µ 2 α !! Z R exp − 1 2 γ a 2 + ζ a ! d a with γ = ( ϕ test 0 ) ⊤ ( K test 1 ) − 1 ϕ test 0 + 1 α ! and ζ = z ⊤ ( σ 2 K test 1 ) − 1 ϕ test 0 + µ α . Then, p Z test ( z ) = 1 (2 π ) ( p + 1) / 2 | ( σ 2 K test 1 ) | 1 / 2 exp − 1 2 z ⊤ ( σ 2 K test 1 ) − 1 z + µ 2 α ! s 2 π γ exp ζ 2 2 γ ! ∝ exp − 1 2 z ⊤ " ( σ 2 K test 1 ) − 1 − 1 γ ( σ 2 K test 1 ) − 1 ϕ test 0 ( ϕ test 0 ) ⊤ ( K test 1 ) − 1 # z ! exp 1 γ z ⊤ ( σ 2 K test 1 ) − 1 ϕ test 0 µ α ! Let us denote K 1 ,α = σ 2 K test 1 + α ϕ test 0 ( ϕ test 0 ) ⊤ . From W oodbury formula, K − 1 1 ,α = ( σ 2 K test 1 ) − 1 − 1 γ ( σ 2 K test 1 ) − 1 ϕ test 0 ( ϕ test 0 ) ⊤ ( σ 2 K test 1 ) − 1 . And as we can remark that 1 γα = 1 − ( ϕ test 0 ) ⊤ ( σ 2 K test 1 ) − 1 ϕ test 0 γ , we hav e 27 1 γ z ⊤ ( σ 2 K test 1 ) − 1 ϕ test 0 µ α = z ⊤ ( σ 2 K test 1 ) − 1 − 1 γ ( σ 2 K test 1 ) − 1 ϕ test 0  ϕ test 0  ⊤ ( σ 2 K test 1 ) − 1 ! µ ϕ test 0 = z ⊤ K − 1 1 ,α µ ϕ test 0 . Finally , p Z test ( z ) ∝ exp − 1 2 z ⊤ K − 1 1 ,α z + z ⊤ K − 1 1 ,α µ ϕ test 0 ! . W e recognize the density of a Gaussian distribution, with mean vector µ ϕ test 0 and covari- ance matrix K α . That shows that Z is a GRF with function mean µϕ 0 and cov ariance kernel αϕ 0 ( s ) ϕ 0 ( s ′ ) + σ 2 K 1 ( s , s ′ ) . Appendix C.2. Distribution of ( Z ( s ) | Y obs = y ) From        Y obs = ( Z ( s i ) + στ E i ) n i = 1 Z is a GRF of mean µϕ 0 and cov ariance kernel σ 2 K 1 ( s , s ′ ) + αϕ 0 ( s ) ϕ 0 ( s ′ ) , we use the conditional Gaussian formulas and obtain that ( Z ( s ) | Y obs = y ) is a GRF with • mean m α ( s ) = µϕ 0 ( s ) + h σ 2 k 1 ( s ) + αϕ 0 ( s ) ϕ obs 0 i ⊤ h σ 2 K 1 ,τ + α ϕ obs 0 ( ϕ obs 0 ) ⊤ i − 1 h y − µ ϕ obs 0 i • cov ariance kernel k α ( s , s ′ ) = σ 2 K 1 ( s , s ′ ) + αϕ 0 ( s ) ϕ 0 ( s ′ ) − h σ 2 k 1 ( s ) + αϕ 0 ( s ) ϕ obs 0 i ⊤ h σ 2 K 1 ,τ + α ϕ obs 0 ( ϕ obs 0 ) ⊤ i − 1 h σ 2 k 1 ( s ′ ) + αϕ 0 ( s ′ ) ϕ obs 0 i Appendix C.3. Limit when α → ∞ Now , let us compute the limits lim α →∞ m α ( s ) and lim α →∞ k α ( s , s ′ ) . First, from W oodbury formula, h σ 2 K 1 ,τ + α ϕ 0 ( ϕ obs 0 ) ⊤ i − 1 = σ − 2 K − 1 1 ,τ        1 − σ − 2 α ϕ obs 0 ( ϕ obs 0 ) ⊤ K − 1 1 ,τ 1 + σ − 2 α ( ϕ obs 0 ) ⊤ K − 1 1 ,τ ϕ obs 0        Then, αϕ 0 ( s )( ϕ obs 0 ) ⊤ h σ − 2 K 1 ,τ + α ϕ 0 ( ϕ obs 0 ) ⊤ i − 1 = σ − 2 αϕ 0 ( s )( ϕ obs 0 ) ⊤ K − 1 1 ,τ − σ − 2 αϕ 0 ( s ) σ − 2 α ( ϕ obs 0 ) ⊤ K − 1 1 ,τ ϕ obs 0 1 + σ − 2 α ( ϕ obs 0 ) ⊤ K − 1 1 ,τ ϕ obs 0 ( ϕ obs 0 ) ⊤ K − 1 1 ,τ = σ − 2 αϕ 0 ( s )( ϕ obs 0 ) ⊤ K − 1 1 ,τ − σ − 2 αϕ 0 ( s )        1 − 1 1 + σ − 2 α ( ϕ obs 0 ) ⊤ K − 1 1 ,τ ϕ obs 0        ( ϕ obs 0 ) ⊤ K − 1 1 ,τ = σ − 2 ϕ 0 ( s )( ϕ obs 0 ) ⊤ K − 1 1 ,τ 1 α + σ − 2 ( ϕ obs 0 ) ⊤ K − 1 1 ,τ ϕ obs 0 28 and σ 2 k 1 ( s ) ⊤ h σ 2 K 1 ,τ + α ϕ 0 ( ϕ obs 0 ) ⊤ i − 1 = k 1 ( s ) ⊤ K − 1 1 ,τ − σ − 2 k 1 ( s ) ⊤ K − 1 1 ,τ ϕ obs 0 ( ϕ obs 0 ) ⊤ K − 1 1 ,τ 1 α + σ − 2 ( ϕ obs 0 ) ⊤ K − 1 1 ,τ ϕ obs 0 Con verg ence of m α ( s ) . Then, with 1 α → 0 , it gives lim α →∞ m α ( s ) = µϕ 0 ( s ) +        k 1 ( s ) ⊤ K − 1 1 ,τ + h ϕ 0 ( s ) − k 1 ( s ) ⊤ K − 1 1 ,τ ϕ obs 0 i ( ϕ obs 0 ) ⊤ K − 1 1 ,τ ( ϕ obs 0 ) ⊤ K − 1 1 ,τ ϕ obs 0        h y − µ ϕ obs 0 i = k 1 ( s ) ⊤ K − 1 1 ,τ y + a ( m ) τ h ϕ 0 ( s ) − k 1 ( s ) ⊤ K − 1 1 ,τ ϕ obs 0 i + µ h ϕ 0 ( s ) − k 1 ( s ) ⊤ K − 1 1 ,τ ϕ 0 − h ϕ 0 ( s ) − k 1 ( s ) ⊤ K − 1 1 ,τ ϕ obs 0 ii = u ⋆ τ ( s ) . Con verg ence of k α ( s , s ′ ) . h σ 2 k 1 ( s ) + αϕ 0 ( s ) ϕ obs 0 i ⊤ h σ 2 K 1 ,τ + α ϕ obs 0 ( ϕ obs 0 ) ⊤ i − 1 h σ 2 k 1 ( s ′ ) + αϕ 0 ( s ′ ) ϕ obs 0 i = k 1 ( s ) ⊤ K − 1 1 ,τ h σ 2 k 1 ( s ′ ) + αϕ 0 ( s ′ ) ϕ obs 0 i + h ϕ 0 ( s ) − k 1 ( s ) ⊤ K − 1 1 ,τ ϕ obs 0 i σ − 2 ( ϕ obs 0 ) ⊤ K − 1 1 ,τ 1 α + σ − 2 ( ϕ obs 0 ) ⊤ K − 1 1 ,τ ϕ obs 0 h σ 2 k 1 ( s ′ ) + αϕ 0 ( s ′ ) ϕ obs 0 i = σ 2 k 1 ( s ) ⊤ K − 1 1 ,τ k 1 ( s ′ ) + h ϕ 0 ( s ) − k 1 ( s ) ⊤ K − 1 1 ,τ ϕ obs 0 i ( ϕ obs 0 ) ⊤ K − 1 1 ,τ 1 α + σ − 2 ( ϕ obs 0 ) ⊤ K − 1 1 ,τ ϕ obs 0 k 1 ( s ′ ) + α k 1 ( s ) ⊤ K − 1 1 ,τ ϕ 0 ( s ′ ) ϕ obs 0 + α h ϕ 0 ( s ) − k 1 ( s ) ⊤ K − 1 1 ,τ ϕ obs 0 i σ − 2 α ( ϕ obs 0 ) ⊤ K − 1 1 ,τ ϕ obs 0 1 + σ − 2 α ( ϕ obs 0 ) ⊤ K − 1 1 ,τ ϕ obs 0 ϕ 0 ( s ′ ) = σ 2 k 1 ( s ) ⊤ K − 1 1 ,τ k 1 ( s ′ ) + h ϕ 0 ( s ) − k 1 ( s ) ⊤ K − 1 1 ,τ ϕ obs 0 i ( ϕ obs 0 ) ⊤ K − 1 1 ,τ 1 α + σ − 2 ( ϕ obs 0 ) ⊤ K − 1 1 ,τ ϕ obs 0 k 1 ( s ′ ) + α k 1 ( s ) ⊤ K − 1 1 ,τ ϕ 0 ( s ′ ) ϕ obs 0 + α h ϕ 0 ( s ) − k 1 ( s ) ⊤ K − 1 1 ,τ ϕ obs 0 i        1 − 1 1 + σ − 2 α ( ϕ obs 0 ) ⊤ K − 1 1 ,τ ϕ obs 0        ϕ 0 ( s ′ ) = σ 2 k 1 ( s ) ⊤ K − 1 1 ,τ k 1 ( s ′ ) + h ϕ 0 ( s ) − k 1 ( s ) ⊤ K − 1 1 ,τ ϕ obs 0 i ( ϕ obs 0 ) ⊤ K − 1 1 ,τ 1 α + σ − 2 ( ϕ obs 0 ) ⊤ K − 1 1 ,τ ϕ obs 0 k 1 ( s ′ ) + αϕ 0 ( s ) ϕ 0 ( s ′ ) − h ϕ 0 ( s ) − k 1 ( s ) ⊤ K − 1 1 ,τ ϕ obs 0 i ϕ 0 ( s ′ ) 1 α + σ − 2 ( ϕ obs 0 ) ⊤ K − 1 1 ,τ ϕ obs 0 And then it comes lim α →∞ k α ( s , s ′ ) = σ 2 K 1 ( s , s ′ ) − σ 2 k 1 ( s ) ⊤ K − 1 1 ,τ k 1 ( s ′ ) + σ 2 h ϕ 0 ( s ) − k 1 ( s ) ⊤ K − 1 1 ,τ ϕ obs 0 i h ( ϕ obs 0 ) ⊤ K − 1 1 ,τ ϕ obs 0 i − 1 h ϕ 0 ( s ′ ) − ( ϕ obs 0 ) ⊤ K − 1 1 ,τ k 1 ( s ′ ) i . 29 Con verg ence to a GRF . F or any finite set of p test points { x 1 , . . . , x p } ⊂ M , let us denote Z test = ( Z ( x 1 ) , . . . , Z ( x p )) ⊤ the corresponding random v ector . Since ( Z | Y obs = y ) is a Gaussian Random Field for any α > 0, the conditional vector ( Z test | Y obs = y ) is Gaussian with a mean vector and a cov ariance matrix that con verge, as α → ∞ , to m lim and K lim respectiv ely . It follows from Lévy’ s continuity theorem that ( Z test | Y obs = y ) L − − − − → α →∞ N ( m lim , K lim ) , (C.1) which ensures that all finite-dimensional distributions of the limiting process are Gaussian, thereby proving that the limit of ( Z | Y obs = y ) remains a GRF . Appendix D . Definition of − ∆ m This section defines the Galerkin approximation − ∆ m of the Laplace-Beltrami opperation − ∆ . Let us denote V m the linear span of ψ 1 , . . . , ψ m the piece wise linear functions associated with the triangulation T . − ∆ m is the endomorphism mapping any f ∈ V m to − ∆ m f , with ∀ u ∈ V m , ⟨− ∆ m f , u ⟩ L 2 ( M ) = ⟨− ∆ f , u ⟩ L 2 ( M ) . Appendix E. Computation of ϕ 0 By definition, the function ϕ ( m ) 0 has unit norm in L 2 ( M ), that is, * m X k = 1 [ ϕ 0 ] k ψ k , m X k = 1 [ ϕ 0 ] k ψ k + L 2 ( M ) = 1 , which implies that ϕ ⊤ 0 M ϕ 0 =  √ M ⊤ ϕ 0  ⊤ √ M ⊤ ϕ 0 = 1 . Since ϕ 0 is constant, it follows that ϕ 0 = 1 m ∥ √ M ⊤ 1 m ∥ 2 . (E.1) Appendix F . In versibility of K ( m ) τ Appendix F .1. F irst scenario As a reminder, Σ is positiv e semidefinite of rank m − 1. Now let us identify a v ector in K er( Σ ) . By definition of the basis function  ψ j  m j = 1 , we hav e ∀ s ∈ M , P m j = 1 ψ j ( s ) = 1 , and then ∀ s ∈ M , m X j = 1 ∇ ψ j ( s ) = 0 . As ϕ 0 is constant, S  √ M  ⊤ ϕ 0 =  √ M  − 1 F  √ M  −⊤  √ M  ⊤ ϕ 0 =  √ M  − 1 F ϕ 0 = 0 . Then,  √ M  ⊤ ϕ 0 ∈ K er( S ) ⊂ Ker( f ( S )) . 30 It comes, Σ M ϕ 0 =  √ M  −⊤ f ( S )  √ M  − 1 √ M  √ M  ⊤ ϕ 0 = 0 . And finally ker( Σ ) = span( M ϕ 0 ) . As a reminder , I = { i 1 , . . . , i n } denotes the set of indices such that s k = c i k for 1 ≤ k ≤ n . Then, ∀ y ∈ R n \ { 0 } , h P ⊤ n y i i =        y k if i = i k for some 1 ≤ k ≤ n , 0 otherwise , and P ⊤ n y < ker( Σ ) . It shows that P n Σ P ⊤ n is positiv e definite. Appendix F .2. Second scenario For the second scenario , the cov ariance matrix P n Σ P ⊤ n + τ 2 I n is inv ertible if and only if − τ 2 is not an eigenv alue of P n Σ P ⊤ n . As Σ is positi ve semidefinite, P n Σ P ⊤ n is positive semidefinite, and − τ 2 is not an eigen value of P n Σ P ⊤ n . Appendix G. Simulations with known mean Here, we provide a strategy to generate simulations of  Z ( m ) | Y ( m ) obs = y , A = a ( m ) τ  . From Proposition 3, we hav e  Z ( m ) | A = a ( m ) τ  is a Gaussian Process with mean a ( m ) τ ϕ 0 ( s ) and covari- ance kernel σ 2 P ( s ) ⊤ Σ P ( s ′ ) . Then, using conditionnal Gaussian process formulas,  Z ( m ) | Y ( m ) obs = y , A = a ( m ) τ  is a GRF with • Mean a ( m ) τ P ( s ) ⊤ ϕ 0 + P ( s ) ⊤ Σ P ⊤ n  K ( m ) τ  − 1 ( y − a ( m ) τ P n ϕ 0 ) • Cov ariance kernel σ 2  P ( s ) ⊤ Σ P ( s ′ ) − P ( s ) ⊤ Σ P ⊤ n  K ( m ) τ  − 1 P n Σ P ( s ′ )  . This GRF has the same mean as  Z ( m ) | Y ( m ) obs = y  in the limit as α → ∞ „ while the cov ari- ance kernels match except for the term c ( m ) ( s , s ′ ), which accounted for the uncertainty in the mean parameter A (see Proposition 4). Therefore, it remains relev ant to generate simulations of  Z ( m ) | Y ( m ) obs = y , A = a ( m ) τ  : the predictiv e mean coincides with the target value, and the associ- ated uncertainty corresponds to the simple kriging framew ork. Let us introduce the following proposition to obtain the target simulations. Proposition 9 (Simulations using the simple kriging framework) . Let Z 0 ∼ N ( 0 m , σ 2 Σ ) and ˜ Z 0 = Z 0 + ϕ 0 ( M ϕ 0 ) ⊤ m ( α ) post ( y ) ( M ϕ 0 ) ⊤ m ( α ) post ( P n ϕ 0 ) . Then, the vector Z SK = ˜ Z 0 + h I m − ϕ 0 ( M ϕ 0 ) ⊤ + h h m ( α ) post ( P n ϕ 0 ) i ( M ϕ 0 ) ⊤ i m ( α ) post  y − P n ˜ Z 0 − στ E  has same distribution as  Z | Y ( m ) obs = y , A = a ( m ) τ  . 31 Pr oof. Let us denote            m SK ( y ) = E  Z | Y ( m ) obs = y , A = a ( m ) τ  = m post ( y ) Σ SK = Cov  Z | Y ( m ) obs = y , A = a ( m ) τ  = σ 2  Σ − Σ P ⊤ n  K ( m ) τ  − 1 P n Σ  . As detailed in [14],                ( M ϕ 0 ) ⊤ m ( α ) post ( x ) = 1 1 α + σ − 2 ( P n ϕ 0 ) ⊤  K ( m ) τ  − 1 P n ϕ 0 σ − 2 ( P ϕ 0 ) ⊤  K ( m ) τ  − 1 x , ∀ x ∈ R m h h m ( α ) post ( P n ϕ 0 ) i = Σ P ⊤ n  K ( m ) τ  − 1 P n ϕ 0 . Then, ˜ Z 0 = Z 0 + a ( m ) τ ϕ 0 . (G.1) From Proposition 6, Z SK = ˜ Z 0 + m post  y − P n ˜ Z 0 − στ E  + h h m ( α ) post  P n ϕ 0  i − ϕ 0  M ϕ 0  ⊤ m ( α ) post  P n ϕ 0  ( M ϕ 0 ) ⊤ m ( α ) post  y − P n ˜ Z 0 − στ E  = ˜ Z 0 + m post  y − P n ˜ Z 0 − στ E  +  h h m ( α ) post  P n ϕ 0  i − ϕ 0  ( P ϕ 0 ) ⊤  K ( m ) τ  − 1  y − P n ˜ Z 0 − στ E  ( P n ϕ 0 ) ⊤  K ( m ) τ  − 1 P n ϕ 0 = ˜ Z 0 + m post  y − P n ˜ Z 0 − στ E  −  ϕ 0 − Σ P ⊤ n  K ( m ) τ  − 1 P n ϕ 0  ( P ϕ 0 ) ⊤  K ( m ) τ  − 1  y − P n ˜ Z 0 − στ E  ( P n ϕ 0 ) ⊤  K ( m ) τ  − 1 P n ϕ 0 = ˜ Z 0 + Σ P ⊤ n  K ( m ) τ  − 1  y − P n ˜ Z 0 − στ E  It comes E ( Z SK ) = m SK ( y ) and Cov( Z SK ) = Cov  ˜ Z 0  − 2Cov  ˜ Z 0 , P n ˜ Z 0 + στ E   K ( m ) τ  − 1 P ⊤ n Σ + Σ P ⊤ n  K ( m ) τ  − 1 Cov  P n ˜ Z 0 + στ E   K ( m ) τ  − 1 P ⊤ n Σ = σ 2 Σ − 2 σ 2 Σ P ⊤ n  K ( m ) τ  − 1 P ⊤ n Σ + σ 2 Σ P ⊤ n  K ( m ) τ  − 1 P ⊤ n Σ = Σ SK Using this result, Algorithm 5 details the procedure to generate simulations of  Z | Y ( m ) obs = y , A = a ( m ) τ  . 32 Algorithm 5 Generation of simple kriging simulations Require: M , F , α, τ , σ , y , a vector ε m of m independent standard Gaussian components 1: Generate a simulation z 0 , centered with cov ariance matrix σ 2 Σ , using Algorithm 1 2: if τ = 0 then 3: Require : I 4: Compute m ( α ) post ( y ), m ( α ) post ( [ z 0 ] I ) and m ( α ) post  P n ϕ 0  using Algorithm 2 5: Compute z SK Proposition 9 6: else 7: Require : P n and ε n a vector of n independent standard Gaussian components 8: Compute m ( α ) post ( y ), m ( α ) post ( P n z 0 + στ ε n ) and m ( α ) post  P n ϕ 0  using Algorithm 3 9: Compute z SK Proposition 9 10: end if 11: retur n z SK Appendix H. Spherical harmonics on the Earth Here we show that the results obtained using the kernel K 1 on the sphere, based on the well-known spherical harmonics, are really close to what we obtain with our finite-element ap- proximation. Figure H.7 provides the SST prediction obtained at the triangulation nodes on the Earth, while Figure H.8 shows the uncertainty associated with the prediction along the longitude closest to 170 ◦ . Figure H.7: Predictions of the SST ( ◦ C ) on Earth, illustrated in 2D using the spherical coordinates ( θ, ϕ ) . n = 20 observation points are shown as black dots. Left: Finite-element approximation of the isotropic spline. Right : Spline using spherical harmonics. 33 (a) Predictiv e distribution with the finite-element isotropic spline. (b) Predictiv e distribution with the kernel based on the spherical harmonics. Figure H.8: Comparison between the uncertainty quantification in SST prediction on the Earth using the finite-element approximation or the kernel based on the spherical harmonics, along the longitude closest to 170 ◦ . 34