title: Unlabeled Data Help in Graph-Based Semi-Supervised Learning: A Bayesian Nonparametrics Perspective
authors: Daniel Sanz-Alonso Ruiyi Yang
text: 
  Theory of three-magnon interaction in a vortex-state magnetic nanodot
      Andrei Slavin
      December 30, 2023
  =====================================================================
  
  
  
  
  
  
  In this paper we analyze the graph-based approach to semi-supervised learning under a manifold assumption. We adopt a Bayesian perspective and demonstrate that, for a suitable choice of prior constructed with sufficiently many unlabeled data, the posterior contracts around the truth at a rate that is minimax optimal  up to a logarithmic factor. Our theory covers both regression and classification. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  § INTRODUCTION
  
  Semi-supervised learning (SSL) refers to machine learning techniques that combine during training a small amount of labeled data with a large amount of unlabeled data. SSL has received increasing attention over the past decades because in many recent applications labeled data are expensive to collect while unlabeled data are abundant. Examples include analysis of body-worn videos and video surveillance <cit.>, text categorization and translation <cit.>, image classification <cit.> and protein structure prediction <cit.>.  The aim of this paper is to investigate whether unlabeled data can enhance learning performance. The answer to such question necessarily depends on the model assumptions and the methodology employed. Our main contribution is to show that under a standard manifold assumption, unlabeled data are helpful when using graph-based methods in a Bayesian setting. We do so by establishing that the optimal posterior contraction rate is achieved (up to a logarithmic factor) provided that the size of the unlabeled dataset grows sufficiently fast with the size of the labeled dataset. 
  
  The SSL problem of interest can be informally described as follows. Given labeled data { (X_1,Y_1),…,(X_n,Y_n) } and unlabeled data {X_n+1,…,X_N_n}, the goal is to predict Y from X. More precisely, we are interested in the  SSL problem of inferring the regression function f_0(x) :=𝔼 (Y|X=x) at the given features _N_n := {X_1, …, X_N_n} in either of these two settings:
  
      
    * Regression: Y=f_0(X)+η, where η∼𝒩(0,σ^2) with σ^2 known.
      
    * Binary classification: ℙ(Y=1|X)=f_0(X). 
  
  
  We analyze graph-based methods applied to this inference task under the manifold assumption  that X takes values on a hidden manifold .  This assumption is natural when the features _N_n live in a high-dimensional ambient space but admit a low-dimensional representation <cit.>.  For instance,  <cit.> demonstrate the low intrinsic dimension of standard datasets for image classification. Graph-based methods are well-suited under the manifold assumption as they allow to uncover the geometry of the hidden manifold and promote smoothness of the inferred f_0 along it. Indeed,  graph-based methods are among the most powerful classification techniques for SSL problems where similar features are expected to belong to the same class <cit.>. The central idea behind traditional graph-based methods is to  infer f_0 by minimizing an objective function comprising at least two terms:  (i) a regularization term constructed with a graph-Laplacian of the features 𝒳_N_n, which leverages the ability of unsupervised graph-based techniques to extract geometric information from ;  and (ii) a data-fidelity term that incorporates the labeled data. We adopt an analogous Bayesian perspective where: (i) the prior distribution (· | _N_n ) will be defined using a graph-Laplacian of the features _N_n  (hence the notation “given _N_n”)  to extract geometric information from ; and (ii) the likelihood function promotes matching the labeled data. Assuming the data are independent so that the likelihood factorizes, the posterior takes the form 
  
      (B  |  _N_n,_n)  ∝∫_B ∏_i=1^n L_Y_i|X_i(f)  d(f  | _N_n ),       B∈ℬ,
  
  where _n := { Y_1, …, Y_n} and ℬ is the Borel σ-algebra on ℝ^N_n (here we are identifying functions over 𝒳_N_n with ℝ^N_n).
  The conditional likelihood of Y_i|X_i is given by
  
      L_Y_i|X_i(f) =1/√(2πσ^2)exp(-|Y_i-f(X_i)|^2/2σ^2)
  
  in the regression setting and 
  
      L_Y_i|X_i(f) =f(X_i)^Y_i[1-f(X_i)]^1-Y_i
  
  in the classification setting. 
  
  
  We shall analyze our Bayesian approach in a frequentist perspective by assuming the data is generated from a fixed truth f_0 and studying the contraction rate of the posterior around f_0, defined as the sequence of real numbers ε_n such that 
  
      𝔼_X𝔼_f_0Π_n(f: d_n(f,f_0)≥ M_nε_n | 𝒳_N_n,𝒴_n) 0,
  
  for any sequence M_n→∞ and some suitable semi-metric d_n. Here the double expectation 𝔼_X𝔼_f_0 is taken first with respect to the conditional data distribution of 𝒴_n|X_1,…,X_n specified by f_0 and then with respect to the marginal of 𝒳_N_n. The idea of posterior contraction rate was formally introduced in <cit.> and has since then become a popular criterion for analyzing Bayesian methods. In particular, contraction of the posterior with rate ε_n gives a point estimator defined as 
  
      f_n=garg max [ Π_n(f:d_n(f,g)≤ M_nε_n | 𝒳_N_n,𝒴_n) ]
  
  that converges to f_0 with the same rate ε_n, which together with the minimax theory for statistical estimation quantifies the performance of our Bayesian SSL procedure.  
  
  
  Our main result, Theorem <ref>, shows that the posterior (<ref>) contracts around f_0  at the minimax optimal rate (up to logarithmic factors), provided that the prior distribution (· | _N_n ) is carefully designed and that N_n grows at a certain polynomial rate with n. More broadly, our theory suggests that graph-based methods require abundant unlabeled data in order to effectively extract geometric information to regularize SSL problems. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Related Work
  
  The question of whether unlabeled data enhance SSL performance has been widely studied.  Positive and negative conclusions have been reached depending on the assumptions made on the relationship between the target function and the marginal data distribution, and also on the methodology employed <cit.>. Two standard model assumptions  are the clustering assumption <cit.>, which states that the target function is smooth over high density regions,  and the manifold assumption that we have introduced. The latter has been extensively used in the statistical learning literature, and specifically in SSL, e.g. <cit.>. 
  
   Several methodologies for SSL have been developed based on generative modeling, support vector machines, semi-definite programming, graph-based methods, etc. (see the overview in <cit.>). Our focus is on  graph-based approaches that combine label information with geometric information extracted from the unlabeled data employing graphical unsupervised techniques <cit.>.  This heuristic has motivated the use of graph-based regularizations in a wide number of applications, but a rigorous analysis of the mechanisms by which unlabeled data enhance the performance of graph-based SSL methods is only starting to emerge. The recent papers <cit.> studied posterior consistency (the asymptotics of the posterior probability Π_n(f:d_n(f,f_0)≥ε | data) for a fixed small number ε)  for a fixed sample size in the small noise limit, whereas we consider the large n limit and further establish posterior contraction rates. Rates of convergence for optimization rather than Bayesian formulations of graph-based SSL have been recently established in  <cit.>.   In a Bayesian framework, <cit.> show the continuum limit of posterior distributions as the size of the unlabeled dataset grows, without increasing the size of the labeled dataset. These papers did not investigate whether the posteriors contract around the truth, and they did not demonstrate the value of unlabeled data in boosting learning performance. The work <cit.> studies fully-supervised function estimation on large graphs assuming that the truth changes with the size of the graph. In contrast, we investigate posterior contraction in a SSL setting with a fixed truth f_0 defined on the underlying manifold .
  
  Our analysis uses tools from Bayesian nonparametrics and spectral analysis of graph-Laplacians. We will provide the necessary background on posterior contraction in Section <ref>, and  we refer  to <cit.>[Chapters 6, 8 and 11] for a comprehensive introduction to this subject. While numerous results on spectral convergence of graph-Laplacians can be found in the literature, our analysis of posterior contraction requires bounds in L^∞-metric with rates, recently developed in <cit.>.
  
  
  
  
  
  
  
  
   §.§ Main Contributions and Scope
  
  
  Our main result, Theorem <ref>, is to our knowledge the first to establish posterior contraction rates for graph-based SSL. In doing so, we provide novel understanding on the relative value of labeled and unlabeled data, and we set forth a rigorous quantitative framework in which to analyze this question. We point out, however, two caveats. First, our theory is non-adaptive in the sense that the Bayesian methodology we analyze only achieves optimal contraction rates when a priori information on the smoothness of the truth is available. Our analysis is motivated by existing graph-based techniques, and the development and analysis of adaptive graph-based methods for SSL is an important research direction stemming from our work, but beyond the scope of this paper. Second, our results build on existing spectral convergence rates of graph-Laplacians that may be suboptimal; as a consequence,  the required sample size of the unlabeled data  that we establish may not be sharp. Due to the plug-in nature of our analysis, improved spectral convergence rates will automatically translate into a sharper bound on the sample complexity. Despite these caveats, our theory provides evidence that letting the unlabeled dataset  grow polynomially with the labeled dataset, as suggested by Theorem <ref>, is a fundamental requirement for standard graph-based methods to achieve optimal contraction. 
  
  A central part of our proof is devoted to analyzing the convergence of a discretely indexed Gauss-Markov random field in an unstructured data cloud to a Matérn-type Gaussian field on . This is formalized in Theorem <ref>, a result that we believe to be of independent interest. Finally, our work contributes to the Bayesian nonparametrics literature on manifolds <cit.>. 
  
  
  
  
    
  Outline
  The rest of this paper is organized as follows. Section <ref> introduces the construction of the graph-based prior and states our main result. Section <ref> provides the necessary background on posterior contraction and outlines our analysis. 
  Section <ref> proves our main result, and we close in Section <ref>.
  
  
  
    
  Notation We denote by ℒ(Z) the law of the random variable Z. For a,b two real numbers, we denote a∧ b= min{a,b} and a∨ b= max{a,b}. The symbol ≲ will denote less than or equal to up to a universal constant. For two real sequences {a_n} and {b_n}, we denote (i) a_n≪ b_n if lim_n (b_n/a_n)=0; (ii) a_n=O(b_n) if lim  sup_n (b_n/a_n)≤ C for some positive constant C; and  (iii) a_n≍ b_n if  c_1≤lim inf_n (a_n/b_n) ≤lim sup_n (a_n/b_n) ≤ c_2 for some positive constants c_1,c_2.  Finally we let γ denote the Lebesgue measure on ℝ.
  
  
  
  § PRIOR CONSTRUCTION AND MAIN RESULT
  
  In this section we introduce the graph-based prior (·  | _N_n) and we state the main result of this paper. Before doing so, we formalize our setting. 
  We assume to be given labeled data (X_1,Y_1),…,(X_n,Y_n)i.i.d.∼Ł(X,Y) and unlabeled data X_n+1,…,X_N_ni.i.d.∼Ł(X), where {X_i}_i=n+1^N_n are independent from {X_i}_i=1^n. 
  Recall that the goal is to estimate the regression function
   
      f_0(x) =𝔼(Y|X=x)
  
  at the given features _N_n = {X_1, …, X_N_n}. 
  We suppose that μ : = Ł(X) is supported on an m-dimensional smooth, connected, compact manifold ℳ without boundary embedded in ℝ^d, with the absolute value of the sectional curvature bounded and with Riemannian metric inherited from the embedding.  For technical reasons, we further assume that m≥ 2 and that  is a homogeneous manifold (the group of isometries acts transitively on , e.g. spheres and tori) normalized so that its volume is 1. We assume for simplicity that μ is the uniform distribution on .  As discussed in Section <ref>,  our results can be generalized  to nonuniform marginal density. 
  
  
  
  
  
   §.§ Graph-based Prior
  
  We now describe the construction of the graph-based prior  (·  | _N_n ) on f_0 restricted to the given features _N_n. The priors we consider have the general form  
  
      (·  | _N_n)=ℒ(Φ( ) | _N_n ),
  
  where  is a Gauss-Markov random field in ℝ^N_n whose covariance will be defined in terms of  a graph-Laplacian <cit.> of _N_n.  For the regression problem, Φ is taken to be the identity. For the classification problem, where f_0 takes values in (0,1),  Φ:ℝ→ (0,1) is a link function, which we assume throughout to be invertible and twice differentiable with Φ^'/(Φ(1-Φ)) uniformly bounded and ∫ (Φ^'')^2/Φ^'  dγ <∞. The logistic function, for instance, satisfies all these standard requirements.
  
  In the remainder of this subsection we introduce and motivate our construction of the Gauss-Markov random field . The starting point is to define a similarity matrix H∈ℝ^N_n× N_n whose entries H_ij≥ 0 represent the closeness between features X_i and X_j. For reasons discussed in Subsection <ref> below, we set
  
      H_ij :=2(m+2)/N_nν_mζ_N_n^m+21{|X_i-X_j|<ζ_N_n},
  
  where |·| is the Euclidean distance in ℝ^d, ν_m is the volume of the m-dimensional unit ball, and ζ_N_n is the connectivity of the graph, a user-specified parameter. Recall that a graph-Laplacian is a positive semi-definite matrix obtained by suitably transforming a similarity matrix. For concreteness, we work with the unnormalized graph-Laplacian matrix Δ_N_n:=D-H, where D is the diagonal matrix with entries D_ii=∑_j=1^N_nH_ij. For a vector v∈ℝ^N_n,  naturally identified with a function on the data cloud _N_n, we have 
  
      v^TΔ_N_nv=1/2∑_i=1^N_nH_ij|v(X_i)-v(X_j)|^2.
  
  We see that indeed Δ_N_n is positive semi-definite and that functions v that change slowly with respect to the similarity H have a small value of v^TΔ_N_nv. This naturally suggests considering the Gaussian distribution 𝒩(0,(I_N_n+Δ_N_n)^-s) since the likelihood for the samples will then have a term similar to (<ref>). Based on this idea we shall define the Gauss-Markov random field   as the principal components of it. Precisely, let {(λ_i^(N_n),ψ_i^(N_n))}_i=1^N_n be the ordered eigenpairs of Δ_N_n and define 
  
      =∑_i=1^k_N_n[1+λ_i^(N_n)]^-s/2ξ_i ψ_i^(N_n),     ξ_ii.i.d.∼𝒩(0,1),
  
  where k_N_n≪ N_n is to be determined. 
  Notice that W_n as in (<ref>) is a truncation of the Karhunen-Loève expansion of samples from  𝒩(0,(I_N_n+Δ_N_n)^-s) and therefore the terminology “principal components”.  
  
  The law of  depends  on three parameters: the graph connectivity ζ_N_n used to define the similarity matrix H, the smoothness parameter s, and the principal components truncation parameter k_N_n. The connectivity ζ_N_n determines the sparsity of the precision matrix of W_n, and it  should be taken to decrease with N_n to better resolve the geometry of  as more unlabeled data are available. The smoothness s controls the level of regularization.  Larger s leads to  faster decay of the coefficients in (<ref>) and smoother samples. The truncation parameter k_N_n≪ N_n allows us to keep only the components of the graph-Laplacian that contain useful geometric information on . Suitable choices and scalings of these three parameters, as well as further insights on their interpretation, will be given  as we develop our theory. Importantly, the construction does not require knowledge of the underlying manifold , but only of its dimension. These data-driven Gauss-Markov fields have been used within various intrinsic approaches to Bayesian learning, see e.g. <cit.>. 
   
  
  
  
  
  
  
  
  
  
    §.§.§ Interpretation as a Discretely Indexed Matérn Field
  
  The Gauss-Markov field  can be interpreted as a discretely indexed Matérn Gaussian field  <cit.>. Consider the Gaussian measure (0,  (I-Δ)^-s) ), where -Δ denotes the Laplace-Beltrami operator on  and the fractional order operator is defined spectrally. Then,  draws from (0,  (I-Δ)^-s) )   admit a representation 
  
      W^=∑_i=1^∞(1+λ_i)^-s/2ξ_i ψ_i,    ξ_ii.i.d.∼𝒩(0,1),
  
  where (λ_i,ψ_i)'s are the ordered eigenpairs of the Laplace-Beltrami operator. Note the analogy with (<ref>). The field W^ is a generalization of the Matérn model to compact manifolds  <cit.>. 
  An important step in our analysis of posterior contraction will be to show the convergence of  towards W^, in a sense to be made precise in Subsection <ref>, provided that the connectivity, smoothness and truncation parameters of  are suitably chosen.
   The similarity matrix defined in (<ref>) enables this convergence result, but other choices (e.g. K-nearest neighbors) are possible. Key to showing convergence of  to W^  is the spectral convergence of Δ_N_n to -Δ. The truncation parameter k_N_n is motivated by the fact that only the first eigenpairs of Δ_N_n accurately approximate those of -Δ; see Proposition  <ref> below for a precise statement.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
    §.§.§ Reconciling Optimization and Bayesian Perspectives
  
  To further motivate the definition of , we relate our prior construction to the optimization literature. To streamline the discussion, let us focus on the regression problem where the link function Φ in (<ref>) is the identity map. Classical graph-based optimization recovers f_0 at _N_n using f^TΔ_N_nf as a regularizer and an appropriate data-fidelity term to match the labeled data,  for instance,
  
      f̂_0:=f∈ℝ^N_narg min ∑_i=1^n |Y_i-f(X_i)|^2 +λ f^T  Δ_N_nf,
  
  where λ controls the level of regularization. The solution to (<ref>) is conceptually equivalent to the maximum a posteriori estimator in the Bayesian approach when the prior is chosen to be 𝒩(0, Δ_N_n^-1), where Δ_N_n^-1 denotes the pseudo-inverse of Δ_N_n. 
  The paper <cit.>  shows that (<ref>) is not well-posed  when m≥ 2, in the sense that as the number of unlabeled data points increases the solution degenerates to a noninformative function. The authors suggest that this issue can be alleviated by defining the regularization term as λ f^TΔ_N_n^s f with s > m/2, so that higher order “derivatives” of f are controlled.  
  Similar behavior has been observed with p-Laplacian regularizations <cit.>. The parameter s in (<ref>) plays the exact same role, and we shall see that suitably scaling s with m is needed to warrant consistent learning in the limit of large unlabeled datasets.
   
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Main Result
  
  Now we are ready to state our main result. Let δ>0 be arbitrary. Let p_m = 3/4 if m=2 and p_m=1/m otherwise. Let α_s=6m+6/2s-3m+1 if s≤9/2m+5/2 and α_s=1 otherwise. 
  
  Suppose Φ^-1(f_0)∈ B_∞,∞^β. Let  be the prior defined by (<ref>) and (<ref>) with s>3/2m-1/2. Consider the following scaling for ζ_N_n,k_N_n and N_n. 
  
      
    * m≤ 4: ζ_N_n≍ N_n^-1/m+4+δ(log N_n)^p_m/2,k_N_n≍ N_n^m/(m+4+δ)(2s-3m+1)α_s(log N_n)^-mp_m/(4s-6m+2)α_s with
      
      N_n≍ n^(m+4+δ)α_s (log n)^(m+4+δ)p_m/2.
  
      
    * m≥ 5: ζ_N_n≍ N_n^-1/2m(log N)^p_m/2,k_N_n≍ N_n^1/(4s-6m+2)α_s(log N_n)^-mp_m/(4s-6m+2)α_s with
      
      N_n≍ n^2mα_s(log n)^mp_m.
  
  
  Then, for _n a multiple of n^-(s-m/2)∧β/2s(log n)^(s-m/2)∧β/4s-2m and all M_n→∞,
  
      𝔼_X𝔼_f_0(f∈ L^∞(μ_N_n): f-f_0_n≥ M_n_n | _N_n,_n )0,
  
  where μ_N_n is the empirical measure of _N_n and  f-f_0^2_n:=1/n∑_i=1^n |f(X_i)-f_0(X_i)|^2. Here 𝔼_f_0 denotes expectation with respect to the conditional distribution of 𝒴_n|X_1,…,X_n determined by f_0 and 𝔼_X denotes expectation with respect to the marginal distribution of 𝒳_N_n.  
  
  
  
  
  Theorem <ref> presents the posterior contraction rates with respect to the priors constructed in Section <ref> under suitable choices of the parameters. Several remarks are in order.
  
  
  
  
  The truth f_0 is assumed to belong to a Besov-type space B_∞,∞^β on manifolds, defined in Subsection <ref> in analogy to the usual Besov space B_∞,∞^β(ℝ^m) on Euclidean space. We recall that B_∞,∞^β(ℝ^m) coincides with the  Hölder space 𝒞^β(ℝ^m) (the space of functions with ⌊β⌋ continuous derivatives whose ⌊β⌋-th derivative is β-⌊β⌋-Hölder continuous) when β∉ℕ and contains 𝒞^β(ℝ^m) when β∈ℕ (see e.g. <cit.>[Theorem 1.3.4]). Therefore,  B_∞,∞^β can be interpreted as an analog of the Hölder space with parameter β on manifolds, and in particular represents a space of β-regular functions. 
  
  The key implication of Theorem <ref> is that when s-m/2=β, we recover the rate n^-β/(2β+m)(log n)^1/4, which is minimax optimal up to a logarithmic factor for β-regular functions. The assumption s>3/2m-1/2 then requires β>m-1/2, so optimal contraction rates can only be attained if f_0 is not too rough. Theorem <ref> can be extended to hold for all s>m if the eigenfunctions of the Laplace-Beltrami operator on  are uniformly bounded, which holds for the family of flat manifolds <cit.> that include for instance the tori. As mentioned in Subsection <ref>, the above choice of s  requires knowing the regularity of f_0 and is not adaptive. 
  
  Another key ingredient of the result is the scaling for N_n as in (<ref>) and (<ref>), which are both larger than a multiple of n^2m since α_s≥ 1 in all cases. In other words, the required sample size of the unlabeled data should grow polynomially with respect to the sample size of the  labeled  data  in order to achieve the near optimal contraction rates described above, thereby justifying the claim that unlabeled data help.  
  
  We further remark that Theorem <ref> only gives an upper bound for the required sample size N_n, whose proof (see Section <ref>) in fact has a plug-in nature. Suppose the sequence of Gaussian-Markov fields  in (<ref>) converges in some semimetric d_n towards W^ in (<ref>) with rate 𝔼_ξ d_n(,W^)=R(N_n) for some function R. The required sample size is then obtained by matching R(N_n) with n^-1. In particular, any improvement of the rate R(N_n) shown in Subsection <ref> will lead to a reduction of the required sample size. However, the spectral convergence rate in Proposition <ref> and hence R(N_n) should not be faster than the resolution of the point cloud, which is shown to scale like N_n^-1/m in Proposition <ref>. Therefore by matching N_n^-1/m with n^-1 we expect a polynomial dependence of N_n on n to be necessary to achieve optimal contraction rate, further demonstrating the need of unlabeled data in our specific setting. 
  
  
  
  
  
  
  § POSTERIOR CONTRACTION: BACKGROUND AND SET-UP
  
  In this section we present necessary background on posterior contraction rates.  In Subsection <ref> we review the main results on posterior contraction that our theory builds on, and in Subsection <ref>  we explain how these general results are used in our proof. 
  
  It will be important to note that our prior is constructed with the X_i's, and we shall view our observations as the Y_i's only, as the notation in (<ref>) suggests. In particular, the Y_i's are independent but non-identically distributed (i.n.i.d.) and hence we will apply general results from <cit.>. This also explains the double expectation in (<ref>), where the randomness of the X_i's is treated separately. 
  
  
  
   §.§ General Principles
  
  Here we review general posterior contraction theory for the problem of estimating f_0 at the continuum level from i.n.i.d data, following <cit.>. 
  Suppose the data {Y_i}_i=1^n are generated according to a density P_f_0^(n)=∏_i=1^n p_f_0,i for some ground truth parameter f_0, where p_f_0,i is the individual density for each observation. Let Π_n be a sequence of priors over f_0 that is supported on some parameter space ℱ equipped with a semimetric d_n. Theorem 4 of <cit.> states that 
  
      𝔼_f_0Π_n(f∈ℱ: d_n(f,f_0)≥ M_n_n | _n) 0,
  
  provided that there exists sets ℱ_n⊂ℱ and positive real numbers _n, K so that 
  
      log N(_n, ℱ_n, d_n)    ≤ n_n^2, 
      Π_n(ℱ_n^c)    ≤  e^-(K+4)n_n^2,
      Π_n(B_n^*(f_0, _n))   ≥ e^-Kn_n^2,
  
  where N(_n,ℱ_n,d_n) is the minimum number of d_n-balls of radius _n needed to cover ℱ_n. Here
  
      B_n^*(f_0,_n):={f∈ℱ: 1/n∑_i=1^n (p_f,i,p_f_0,i)≤_n^2,   1/n∑_i=1^n v_(p_f,i,p_f_0,i)≤ C_n^2},
  
  where C is a universal constant, (p,q)=∫ p log (p/q)dγ is the Kullback-Leibler divergence and v_(p,q)=∫ p|log(p/q)-(p,q)|^2dγ.  Conditions (<ref>) and (<ref>) roughly state that there are approximating sieves ℱ_n which capture most of the prior probability while not being too large.
   Condition (<ref>) requires sufficient prior mass near the truth f_0 and together with (<ref>) further indicates that the priors should be “uniformly spread”. 
    In fact the three conditions are stronger than those in <cit.>[Theorem 4] but will suffice in our case for Gaussian process priors to be discussed shortly below. We shall refer to <cit.>[Section 2] and <cit.>[Section 2] for further discussion on the interpretation and relaxation of these conditions.
  
  The general theory can be used to analyze a wide range of statistical models but does not give a recipe for constructing the sieves ℱ_n and _n. However, when the priors are Gaussian there exists a simple relation between _n and the priors. Suppose in addition that Π_n=ℒ(w_n) are Gaussian measures on some Banach space (𝔹,·_𝔹) that converge to a fixed Gaussian measure Π=ℒ(w) on the same Banach space with 10𝔼w_n-w^2_𝔹≤ n^-1. Theorem 2.2 in <cit.> then states that the  contraction rate ε_n can be obtained by studying the concentration function of the limit prior Π, defined as 
  
      ϕ_f_0() := h∈ℍ:h-f_0_𝔹<inf h^2_ℍ - logℙ (w_𝔹<),
  
  where (ℍ,·_ℍ) is the  reproducing kernel Hilbert space (RKHS) of Π (see e.g. <cit.> for more details on Gaussian measures and the associated RKHSs). More precisely, if _n satisfies ϕ_f_0(_n)≤ n_n^2, then for the same _n and some ℱ_n, the three conditions (<ref>), (<ref>) and (<ref>) are satisfied (possibly for different proportion constants in front of n_n^2) with d_n and B_n^*(f_0,_n) replaced by ·_𝔹 and {f:f-f_0_𝔹<_n}, respectively.  
  
  
  
   
  
  
   §.§ Application to Our Setting
   
  In this subsection we discuss how we utilize the general theory outlined above, and provide a road map for the proof of Theorem <ref>. Notice that in (<ref>) the sequence of posteriors are supported on the discrete space L^∞(μ_N_n), whose size changes with n. To alleviate this issue we will reduce the analysis to a sequence of posteriors Π_n^(· | _N_n,_n) supported on the same continuum space L^∞(μ) in Subsection <ref>. The corresponding sequence of priors will turn out to have the form  Π_n^( ·  | _N_n) = ℒ(Φ(ℐW_n)), with W_n defined in (<ref>) and ℐ a suitable interpolation map so that ℐW_n is a Gaussian process on  that approximates W^. Our analysis will then consist of the following steps. We will establish the three conditions (<ref>), (<ref>), (<ref>) for W^ in Subsection <ref>, followed by a convergence rate analysis of ℐW_n towards W^ in Subsection <ref>, so that the three conditions are inherited by ℐW_n. The assumptions on the link function Φ will then allow us to conclude similar conditions for Φ(ℐW_n). 
  
  The discrepancy measure d_n can be taken as the empirical L^2-norm ·_n in the regression case (see e.g. <cit.>[Section 8.3.2]) and the root average square Hellinger distance in the classification setting (<cit.>[Section 3]), defined as 
  
      d^2_n,H(f,f^')=1/n∑_i=1^n ∫(√(p_f,i)-√(p_f^',i))^2 dγ.
  
  In the latter case, one can show that d_n,H is equivalent to ·_n. Indeed, since the densities are uniformly bounded, f-f^'_n is bounded above by a multiple of (<ref>). Furthermore (<ref>) is upper bounded by a multiple of f-f^'_n by <cit.>[Lemma 2.8(ii)] given our assumption on Φ. Hence this explains the choice of ·_n in (<ref>). 
  
  
  A natural choice for the Banach space in our setting is then  𝔹=L^∞(μ) since (i) d_n=·_n≤·_L^∞(μ); and (ii) B_n^*(f_0,_n) ⊃{f∈ℱ: f-f_0_L^∞(μ)≤C̃_n} for some universal constant C̃. The first point is clear and an upper bound on the metric entropy (<ref>) in ·_L^∞(μ) will automatically yield an upper bound in ·_n.  The second point follows from the fact that one can upper bound the two quantities in (<ref>) by f-f_0_n and hence f-f_0_L^∞(μ) in both the regression and classification setting. Therefore the prior mass condition (<ref>) for W^ in L^∞(μ) balls is sufficient to give the corresponding condition for B^*(f_0,_n). This motivates us to consider the continuum Gaussian field defined in (<ref>) as an element in L^∞(μ) in Subsection <ref> and the L^∞(μ) convergence rate in Subsection <ref>. 
  
  
  
  
  
  
  
  
  
  
  
  § PROOF OF THE MAIN RESULT
   
  In this section we prove Theorem <ref>. An important part of the proof is to formalize the convergence of the Gauss-Markov random field  in (<ref>), defined in the data cloud _N_n, to the Matérn field W^ in (<ref>), defined in .  To that end, we introduce an interpolation of  into 
  
      :=∑_i=1^k_N_n[1+λ_i^(N_n)]^-s/2ξ_i ψ_i^(N_n)∘ T_N_n,     ξ_ii.i.d.∼𝒩(0,1),
  
  where T_N_n:→{X_1,…,X_N_n} are transport maps to be specified in Proposition <ref> below. In Subsection <ref> we show that it suffices to establish posterior contraction with respect to the prior 
  
      ( · | _N_n) :=ℒ(Φ()  | _N_n).
  
   In Subsection <ref> we show concentration properties of the limiting Gaussian field W^. In Subsection <ref> we establish the L^∞(μ)-convergence rate of  towards W^.
  Combining the facts that  is close to W^ with the contraction properties of W^, we deduce that Ł( | _N_n ) satisfies the three conditions (<ref>), (<ref>), (<ref>). We complete the proof in Subsection <ref> by  lifting these conditions to Ł(Φ()  | _N_n ) in both the regression and classification problems. 
   
   
  
  
   §.§ Reduction via Interpolation
  
  
  Here we show that in order to establish Theorem <ref> it suffices to prove posterior contraction with respect to the continuum prior ( · | _N_n) defined in (<ref>). Let ℐ:=ℐ_N_n:L^∞(μ_N_n)→ L^∞(μ) be the interpolation map defined by ℐu:=u∘ T_N_n. Specifically, we claim that, in order to establish Theorem <ref>, it suffices to show that
  
      𝔼_X𝔼_f_0(f∈ℐ(L^∞(μ_N_n)):f-f_0_n≥ M_n_n | _N_n , _n)  0.
  
  The only property of the maps T_N_n in (<ref>) that we will use in this subsection is that T_N_n(X_i) = X_i.  
  In order to establish the claim, let A_n={f∈ L^∞(μ_N_n):f-f_0_n≥ M_n_n}. Since A_n⊂ℐ^-1(ℐ(A_n)) we have 
  
      (A_n | _N_n,_n) ≤(ℐ^-1(ℐ(A_n)) | _N_n,_n)=ℐ_♯[(· | _N_n,_n)](ℐ(A_n)),
  
  where ℐ_♯ denotes push-forward through the map ℐ. (Recall that the push-forward of a measure ν through the map ℐ is defined by the relationship (ℐ_♯ν)(B):=ν(ℐ^-1(B)).)  Therefore, it follows from Lemma <ref> below that 
  
      (A_n | _N_n,_n)  ≤( ℐ (A_n) | _N_n, _n),
  
  and hence it suffices to bound the right-hand side of (<ref>). 
  Since T_N_n(X_i)=X_i, we have
  
      ℐ(A_n)={f∘ T_N_n: f-f_0_n≥ M_n_n }   ={f∘ T_N_n: f∘ T_N_n-f_0_n≥ M_n_n }
         ={f∈ℐ(L^∞(μ_N_n)): f-f_0_n≥ M_n_n},
  
  and the claim is established. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  It holds that  (· | _N_n,_n)=ℐ_♯[ (· | _N_n, _n)].
  
  
  Step 1:  First we show that ( · | _N_n)  = ℐ_♯ [ ( ·  | _N_n)].
  Since ℐ is linear, we have that = ℐ. Then observe that ℐ(Φ())=Φ(ℐ()). Indeed, for x∈ T_N_n^-1({X_i}), we have
  
      ℐ(Φ())(x)=Φ()(X_i)=Φ((X_i))=Φ(ℐ()(x))=Φ(ℐ())(x).
  
  Therefore, we have that 
  
      ( · | _N_n) = ℒ(Φ() | _N_n) =ℒ(ℐ(Φ()) | _N_n).
  
  Finally, note that ℐ_♯ [( ·  | _N_n) ]=  ℒ(ℐ(Φ()) | _N_n), since
  
      ℐ_♯ (B | _N_n) = (ℐ^-1(B) | _N_n) = ℙ(Φ()∈ℐ^-1(B) | _N_n)= ℙ(ℐ(Φ())∈ B | _N_n).
  
  Step 2: Now we show that (· | _N_n,_n)=ℐ_♯[ (· | _N_n, _n)].
  By definition of pushforward measure, it suffices to show that, for any measurable B,
  
      (B | _N_n,_n)=(ℐ^-1(B) | _N_n,_n).
  
  Using Step 1, the left-hand side equals 
  
      (B | _N_n,_n) 
             =
          ∫_B ∏_i=1^n L_Y_i|X_i(f)  dℐ_♯[(f  | _N_n)]/∫_ℐ(L^∞(μ_N_n))∏_i=1^n L_Y_i|X_i(f)  dℐ_♯[(f | _N_n)],
  
  where 
  L_Y_i|X_i(f) is given in (<ref>) and (<ref>). Note that pointwise values of f are well-defined since ℐ_♯ is supported on ℐ(L^∞(μ_N_n)).
  By the change-of-variable formula for pushforward measures, 
  
      (<ref>) = ∫_ℐ^-1(B)∏_i=1^n L_Y_i|X_i∘ℐ(f_n)  d(f_n | _N_n) /∫_L^∞(μ_N_n)∏_i=1^n L_Y_i|X_i∘ℐ(f_n) d(f_n | _N_n),
  
  which equals (<ref>) with B replaced by ℐ^-1(B), by noticing that L_Y_i|X_i∘ℐ(f_n) is exactly the same as in (<ref>) and (<ref>). The result follows. 
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
   §.§ Regularity and Contraction Properties of the Limiting Field
  
  In this subection we study  the limit Gaussian measure π=𝒩(0,(I-Δ)^-s)=ℒ(W^). Recall that the samples admit Karhunen-Loève expansion
  
      W^=∑_i=1^∞(1+λ_i)^-s/2ξ_i ψ_i,    ξ_ii.i.d.∼𝒩(0,1),
  
  where (λ_i,ψ_i)'s are the eigenpairs of the Laplace-Beltrami operator -Δ on . Notice that a larger s leads to a faster decay of the coefficients and hence more regular sample paths. From  Weyl's law that λ_i≍ i^2/m (see e.g. <cit.>[Theorem 72]), setting s>m/2 makes π a well-defined measure on L^2(μ). To ensure almost sure continuity of the samples, so that point evaluations are well-defined and π is supported over L^∞(μ), we need the stronger assumption that s>m-1/2. 
  
  If s>m-1/2, then samples of π are almost surely continuous. 
   
  
  By <cit.>[Corollary 4.5], it suffices to show that
   
      𝔼_ξ |W^(x)-W^(y)|^2 ≤ Cd_(x,y)^η,
  
  for some η>0. We have
  
      𝔼_ξ|W^(x)-W^(y)|^2     =𝔼_ξ(∑_i=1^∞(1+λ_i)^-s/2ξ_i(ψ_i(x)-ψ_i(y)))^2 
         =∑_i=1^∞(1+λ_i)^-s|ψ_i(x)-ψ_i(y)|^2,
  
  since 𝔼_ξξ_iξ_j=0 for i≠ j. By Proposition <ref>,
  
      |ψ_i(x)-ψ_i(y)|    ≤ Cλ_i^m-1/4,
      
          |ψ_i(x)-ψ_i(y)|    ≤∇ψ_i _L^∞(μ) d_(x,y) ≤ Cλ_i^m+1/2 d_(x,y).
  
  Using Weyl's law, we have  
  
      𝔼_ξ|W^(x)-W^(y)|^2    ≲∑_i=1^∞ (1+λ_i)^-smin{λ_i^m+1d_(x,y)^2,λ_i^m-1/2}
         ≲∑_i=1^∞i^-2s/mmin{i^2m+2/md_(x,y)^2, i^m-1/m}
         ≲ d_(x,y)^2∫_1^K z^-2s+2m+2/m dz +∫_K^∞ z^-2s+m-1/mdz.
  
  Setting K= d_(x,y)^-2m/m+3, we have 
  
      𝔼_ξ|W^(x)-W^(y)|^2    ≲  d_(x,y)^4s-4m+2/m+3.
  
  The result  follows since  s>m-1/2.
  
  
  
  The argument above relies on the worst case L^∞(μ) bound on the eigenfunctions of -Δ given in Proposition <ref>. If the eigenfunctions are uniformly bounded (which holds for the family of flat manifolds <cit.> that includes, for instance, the torus), then continuity can be guaranteed for s>m/2. 
  
  
  From now on, we shall consider π as a measure over L^∞(μ) with continuous sample paths, where pointwise evaluation is well-defined. Using the series representation (<ref>), the reproducing kernel Hilbert space ℍ associated with π has the following characterization 
  
      ℍ={h=∑_i=1^∞ h_i  ψ_i: h_ℍ^2:=∑_i=1^∞  h_i^2   (1+λ_i)^s <∞}.
  
  From the general theory in <cit.>, the concentration properties of W can be characterized by the concentration function 
  
      ϕ_w_0() := h∈ℍ:h-w_0_L^∞(μ)<inf h^2_ℍ - logℙ (W_L^∞(μ)<),
  
  where w_0 belongs to the support of π.  For our purpose, w_0 will be set as Φ^-1(f_0). The three conditions (<ref>), (<ref>) and (<ref>) hold for the _n that satisfies ϕ_w_0(_n)≤ n_n^2. 
  
  Before stating our main result in this subsection, we follow <cit.> to define a Besov space which will be needed to characterize the regularity of w_0. Let Ψ be an even function in the Schwartz space 𝒮(ℝ) with 
  
      0≤Ψ≤ 1,  Ψ≡ 1  on   [-1/2,1/2],   supp(Ψ)⊂[-1,1].
  
  Define the Besov space 
  
      B_∞,∞^β :={w: w_B_∞,∞^β:= j∈ℕsup  2^β jΨ_j(√(Δ))w(·)-w(·)_L^∞(μ) <∞},
  
  where Ψ_j(·)=Ψ(2^-j·) and, for w=∑_i=1^∞w_i ψ_i, 
  
      Ψ_j(√(Δ))w:=∑_i=1^∞Ψ_j(√(λ_i))w_iψ_i.
  
  It was shown in <cit.>[Proposition 6.2] that the definition is independent of the choice of Ψ. 
  
  
  
  Suppose w_0∈ B_∞,∞^β. Consider the prior π=𝒩(0,(I-Δ)^-s) with  s>β∧ m-1/2.  Then there exists sets B_n⊂ L^∞(μ) so that
  
      log N(3_n,B_n,·_L^∞(μ))   ≤  6Cn_n^2, 
      π(B_n^c)   ≤ e^-Cn_n^2, 
      π(w-w_0_L^∞(μ)<2_n)   ≥ e^-n_n^2,
  
  where _n is a multiple of n^-(s-m/2)∧β/2s(log n)^(s-m/2)∧β/4s-2m  and C>1 is a constant satisfying e^-Cn_n^2<1/2. 
  
  
  As noted above, it suffices to show that the above _n satisfies ϕ_w_0(_n) ≤ n_n^2, where ϕ_w_0 is defined in (<ref>). To bound the small ball probability, it suffices by general results from <cit.> to bound the metric entropy of ℍ_1, the unit ball in ℍ. Notice that (<ref>) has the form of a Sobolev ball of regularity s. Classical results for Sobolev spaces on ℝ^m such as <cit.>[Theorem 3.3.2] have shown that the entropy in the L^∞ metric is bounded by a multiple of ^-m/s. The manifold case has been treated in <cit.>[Theorem 3.12]  assuming homogeneity, which includes an additional logarithmic factor, i.e., 
  
      log N(, ℍ_1,·_L^∞(μ))≲^-m/s(log1/)^1/2.
  
  The original theorem was about entropy numbers but can be translated to the above statement on metric entropy and its proof suggested that when p=2  the theorem holds with s>m/2. In particular their Sobolev class is defined as
  
      I_s U_2:={h=h_1+∑_i=2^∞λ_i^-s/2h_iψ_i: ∑_i=1^∞h_i^2≤ 1},
  
  which is compatible with our ℍ. (We remark that although <cit.> defined their Sobolev class with a further mean-zero condition, their proof actually applied to the general case.) 
  Now by (<ref>) and <cit.>[Theorem 1.2] we have
  
      -logℙ(w_L^∞(μ)<) ≲^-2m/2s-m(log1/)^s/2s-m.
  
  
  
  For the decentering function, let C_0=w_0_B^β_∞,∞ and consider h=Ψ_J(√(Δ))w_0 with J large enough so that C_02^-β J≤. Since w_0∈ B_∞,∞^β, we have
  
      Ψ_j(√(Δ))w_0-w_0_L^∞(μ)≤ C_02^-β j,      j∈ℕ.
  
  In particular h-w_0_L^∞(μ)≤ C_0 2^-β J≤. Suppose w_0 has the representation w_0=∑_i=1^∞w_iψ_i. We then have h=∑Ψ_J(√(λ_i)) w_i ψ_i∈ℍ since Ψ_J(√(λ_i))=0 for √(λ_i)>2^J. Moreover, since Ψ_J≤ 1,
  
      h_ℍ^2 ≤∑_√(λ_i)≤ 2^J w_i^2 (1+λ_i)^s≤∑_√(λ_i)≤ 12^sw_i^2+∑_j=1^J ∑_2^j-1<√(λ_i)≤ 2^j w_i^2(1+λ_i)^s.
  
  By (<ref>), we have 
  
      ∑_2^j-1< √(λ_i)≤ 2^j w_i^2(1+λ_i)^s ≲ 2^2js∑_2^j-1< √(λ_i) w_i^2    ≤2^2jsΨ_j-1(√(Δ))w_0-w_0_2^2
         ≲ 2^2jsΨ_j-1(√(Δ))w_0-w_0_∞^2≲ 2^2(s-β)j.
  
  Therefore recalling that C_02^-β J≤,
  
      (<ref>)≲∑_√(λ_i)≤ 12^sw_i^2 +∑_j=1^J 2^2(s-β)j≲ 2^2(s-β)J≲^-2(s-β)/β,
  
  since the first sum remains bounded as J→∞. Combining the above with (<ref>), we get 
  
      ϕ_w_0() ≲^-2m/2s-m(log1/)^s/2s-m + ^-2(s-β)/β,
  
  which implies 
  
      ϕ_w_0()/^2≲^-2s/s-m/2(log1/)^s/2s-m + ^-2s/β.
  
  Therefore by choosing 
  
      _n =Cn^-(s-m/2∧β)/2s (log n)^(s-m/2)∧β/4s-2m
  
  for a large enough constant C, the condition ϕ_w_0(_n)≤ n_n^2 is satisfied. The three assertions then follow from <cit.>[Theorem 2.1].
  
  
  
  
  
  
   §.§ Convergence of Gaussian Fields in L^∞(μ)
   
  In this subsection we establish L^∞(μ) convergence  of  to W^.  We shall denote N:=N_n throughout the rest of this subsection to simplify our notation.
  Recall that 
  
      =∑_i=1^k_N[1+λ_i^(N)]^-s/2ξ_i ψ_i^(N)∘ T_N,     ξ_ii.i.d.∼𝒩(0,1),  
      
           W^   =∑_i=1^∞(1+λ_i)^-s/2ξ_i ψ_i,    ξ_ii.i.d.∼𝒩(0,1),
  
  In order to show the convergence, the transportation maps T_N will be assumed to be close to the identity in the sense made precise in the following proposition from <cit.>[Theorem 2].
   
  For γ>1, there exists a transportation map T_N:→{X_1,…,X_N}  satisfying T_N_n(X_i)=X_i  so that, with probability 1-O(N^-γ),
  
      ρ_N:=x∈ℳsup  d_ (x,T_N(x))≲(log N)^p_m/N^1/m,
  
  where p_m = 3/4 if m=2 and p_m=1/m otherwise. 
  
  The transport map T_N_n is a measure-preserving transformation in the sense that μ(T_N_n^-1(U))=μ_N_n(U) for all U⊂ measurable.  The additional requirement T_N(X_i)=X_i that was absent in <cit.>[Theorem 2] is valid here  since modifying T_N on a set of μ-measure 0 does not affect the measure preserving property and (<ref>) still holds in this case. The scaling in (<ref>) can be thought as the resolution of the point cloud,   and is important in suitably defining the choice of connectivity ζ_N, as we shall see. 
  
  
  
  
  
  
  
  
  
  
  
  The main result of this subsection is the following. 
  
  Let δ>0 be arbitrary.
  
      
    * 3/2m-1/2<s≤9/2m+5/2. Set
      
      ζ_N≍ N^-1/m+4+δ(log N)^p_m/2,     k_N≍ N^1/(m+4+δ)(6+6/m)(log N)^-mp_m/12m+12,         m≤ 4,
      ζ_N≍ N^-1/2m(log N)^p_m/2,    k_N≍ N^1/12m+12(log N)^-mp_m/12m+12,        m≥ 5.
  
  Then, with probability tending to 1,
  
      𝔼_ξ-W^^2_L^∞(μ)≲
          N^-2s-3m+1/(m+4+δ)(6m+6)(log N)^p_m(2s-3m+1)/12m+12,        m≤ 4,
      
          N^-2s-3m+1/m(12m+12) (log N)^p_m(2s-3m+1)/12m+12,       m≥ 5.
  
  
  
    * s>9/2m+5/2. Set
  
      ζ_N≍ N^-1/m+4+δ(log N)^p_m/2,    k_N≍ N^m/(m+4+δ)(2s-3m+1)(log N)^-mp_m/4s-6m+2,         m≤ 4, 
      ζ_N≍ N^-1/2m(log N)^p_m/2,    k_N≍ N^1/4s-6m+2(log N)^-mp_m/4s-6m+2,         m≥ 5.
  
  Then, with probability tending to 1,
  
      𝔼_ξ-W^^2_L^∞(μ)≲
          N^-1/m+4+δ(log N)^p_m/2,       m≤ 4,
      
          N^-1/2m(log N)^p_m/2,       m≥ 5.
  
  
  
  
  We remark that the statement “with probability tending to 1” refers to the randomness coming from the X_i's. 
  By solving for N so that the rate matches n^-1, we get the following. 
  
  Let δ>0 be arbitrary. Let
  
      N≍ n^α_1(log n)^α_2,
  
  where 
  
      α_1=(m+4+δ)(6m+6)/2s-3m+1,  α_2=p_m(m+4+δ)/2,       if   m≤ 4   ,  3/2 m-1/2<s≤9/2m+5/2,
      α_1=m(12m+12)/2s-3m+1,  α_2=mp_m,       if   m ≥ 5    ,  3/2 m-1/2<s≤9/2m+5/2,
      α_1=m+4+δ,  α_2=p_m(m+4+δ)/2,        if   m≤ 4    ,   s> 9/2m+5/2,
      α_1=2m,  α_2=mp_m,        if   m≥ 5    ,   s> 9/2m+5/2.
  
  Let ζ_N and k_N have the same scaling as in Theorem <ref>. Then, with probability tending to 1,
  
      𝔼_ξ-W^^2_L^∞(μ)≲ n^-1.
  
  
  
  
  For flat manifolds the results of Theorem <ref> and Corollary <ref> can be shown to hold for s>m with corresponding modifications in the scaling of the parameters.
  
  
  The key to show the above results is to derive convergence rates for 
  
      |λ^(N)_i-λ_i|,    and    ψ_i^(N)∘ T_N-ψ_i_L^∞(μ).
  
  We shall build our analysis on the existing results from the literature.
  Recall that ζ_N is the connectivity of the graph and the resolution ρ_N  defined in (<ref>). Assuming we are in the event that (<ref>) holds, the following results from <cit.>[Theorem 4.6 & 4.7] bound the spectral approximations. 
   
  Suppose  ρ_N≪ζ_N and ζ_N√(λ_k_N)≪ 1 for N large. Then there exists orthonormalized eigenfunctions {ψ_i^(N)}_i=1^N and {ψ_i}_i=1^∞ so that, for i=1,…,k_N,
  
      |λ_i^(N)-λ_i|    ≲λ_i(ρ_N/ζ_N+ζ_N√(λ_i)),
      ψ_i^(N)∘ T_N-ψ_i_L^2(μ)^2    ≲  i^3  (ρ_N/ζ_N+ζ_N√(λ_i)).
  
  
  
  
  Proposition <ref> bounds the eigenfunction approximation in L^2(μ) norm and we need to lift such result to the L^∞(μ_N) norm. We remark that L^∞(μ_N) convergence rates have already been established in the literature (see e.g. <cit.>). However these results do not contain an explicit proportion constant in terms of the index i as in (<ref>). Instead of going through their details to find the explicit constants, which would make our presentation much more involved, we directly show L^∞(μ_N) bounds based on Proposition <ref>, by following the same idea as in <cit.>. Since we build our results from (<ref>), which is not the sharpest bound in the literature, our results in Theorem <ref> and Corollary <ref> suffer the same drawback. Nevertheless, the goal of this paper is to demonstrate the idea that unlabeled data helps and hence the finding sharpest scaling of N_n is less essential.  
  
  
  
  
  Below we record four results from <cit.>[Theorem 3.3], <cit.>[Theorem 1.2] and <cit.>[Equation (2.10)],  <cit.>[Corollary 2.5] that will be needed. 
  
   
  Let f∈𝒞^3(). Then with probability  1-2nexp(-cNζ_N^m+4),
  
      Δ_Nf(x)-Δ f(x)_L^∞(μ_N)≤ C(1+f_𝒞^3())ζ_N.
  
  
  
  
  
   
  Let ψ be an L^2(μ)-normalized eigenfunction of -Δ associated with eigenvalue λ≠ 0. Then, for k∈ℕ,
  
      ψ_L^∞(μ)   ≤ C λ^m-1/4, 
      ∇^k ψ_L^∞(μ)   ≤ Cλ^k+m-1/2.
  
  
  
  
   
  Suppose ζ_N≤c/Λ+1 where c is a sufficiently small constant depending only on . Then with probability 1-Cζ_N^-6mexp(-cNζ_N^m+4)-2Nexp(-cN(Λ+1)^-m) 
  
      u_L^∞(μ_N)≤ C(Λ+1)^m+1u_L^1(μ_N)
  
  for all u∈ L^2(μ_N) with λ_u<Λ, where 
  
      λ_u :=Δ_N u_L^∞(μ_N)/u_L^∞(μ_N).
  
  
  Notice that the high probability condition is satisfied only if 
  
      ζ_N^-6mexp(-cNζ_N^m+4)-2Nexp(-cN(Λ+1)^-m)→ 0
  
  and this requires some care in setting the scaling of ζ_N and k_N. 
  In particular a sufficient condition for ζ_N is ζ_N ≳ N^-1/m+4+δ for an arbitrarily small δ>0.  On the other hand the condition for Λ reduces to the scaling of k_N. 
  Indeed if u=ψ_i^(N), then λ_u=λ^(N)_i and therefore Λ can be chosen as a multiple of λ_i given the spectral approximation in Proposition <ref>. We will show that the same choice of Λ suffices when u=ψ_i^(N)-ψ_i. Since we are interested in bounding the first k_N eigenfunctions, we can set Λ to be multiple of λ_K_N. Therefore in order to make 2Nexp(-cN(Λ+1)^-m)≪ 1, it is sufficient to have λ_k_N≲ N^1-δ/m, i.e., if k_N ≲ N^1-δ/2 by Weyl's law. Therefore we should keep the following in mind when we set the scaling later
  
      N^-1/m+4+δ≲ζ_N,     k_N ≲ N^1-δ/m,    ζ_N k_N^2/m≲ 1,
  
  where the third requirement corresponds to the condition that ζ_N (Λ+1)≤ c.
  
  Now we are ready to show the L^∞(μ_N) bound of eigenfunction approximations.  
  
  Under the same conditions and the intersection of the high probability events as in Propositions <ref> and <ref> we have, for i=1,…,k_N,
  
      ψ_i^(N)∘ T_N-ψ_i_L^∞(μ)≲λ_i^m+1i^3/2√(ρ_N/ζ_N+ζ_N√(λ_i)).
  
  
  
  
  We will first modify the L^2(μ)-bound (<ref>) to an L^2(μ_N)-bound using the regularity of eigenfunctions of -Δ, which is then lifted to an L^∞(μ_N)-bound using Proposition <ref> and <ref>. Finally we use regularity of the eigenfuctions again to transfer the L^∞(μ_N)-bound to an L^∞(μ)-bound. To start, we notice that the transport T_N induces a partition of  by the sets {U_i=T_N^-1({X_i})}_i=1^N. Furthermore, the measure preserving property gives μ(U_i)=1/N and Proposition <ref> implies that U_i⊂ B_(X_i,ρ_N). We then have
  
      ψ_i^(N)-ψ_i^2_L^2(μ_N)=1/N∑_i=1^N      |ψ_i^(N)(X_i)-ψ_i(X_i)|^2=∑_i=1^N ∫_U_i|ψ_i^(N)(X_i)-ψ_i(X_i)|^2 dμ(x)  
         ≤∑_i=1^N ∫_U_i 2|ψ_i^(N) (X_i) -ψ_i(x)|^2+2|ψ_i(x)-ψ_i(X_i)|^2 dμ(x),
  
  where 
  
      ∫_U_i |ψ_i(x)-ψ_i(X_i)|^2 dμ(x) ≤1/N∇ψ_i^2_∞ ρ_N^2.
  
  Hence by Proposition <ref> and Weyl's law
  
      ψ_i^(N)-ψ_i^2_L^2(μ_N)≤ 2ψ_i^(N)∘ T_N-ψ_i^2_L^2(μ) +Ci^m+1/mρ_N^2.
  
  Since i^m/m+1ρ_N^2 is of higher order than (<ref>), we conclude that 
  
      ψ_i^(N)-ψ_i^2_L^2(μ_N)≲  i^3  (ρ_N/ζ_N+ζ_N√(λ_i)).
  
  Now let g=ψ_i^(N)-ψ̃_i with ψ̃_̃ĩ=ψ_i|_{X_1,…,X_N}. We have by Proposition <ref>
  
      Δ_N g_L^∞(μ_N)   ≤Δ_N ψ_i^(N) - Δψ_i_L^∞(μ_N) +Δψ_i - Δ_N ψ̃_i_L^∞(μ_N)
         ≤λ_i^(N)ψ_i^(N)-λ_i ψ_i_L^∞(μ_N)+C(1+ψ_i_𝒞^3())ζ_N
         ≤λ_i^(N)ψ_i^(N)-ψ_i_L^∞(μ_N) + |λ_i^(N)-λ_i|ψ_i_L^∞(μ_N)+ C(1+ψ_i_𝒞^3())ζ_N
         ≤λ_i^(N)g_L^∞(μ_N)+Cλ_i^m+5/2(ρ_N/ζ_N+ζ_N),
  
  where we have used (<ref>) and Proposition <ref> in the last step. Therefore recalling the definition in (<ref>), we have 
  
      λ_g ≤λ_i^(N) +Cg^-1_L^∞(μ_N)λ_i^m+5/2(ρ_N/ζ_N+ζ_N).
  
  If λ_g≥λ_i^(N)+1, then we have 
  
      g_L^∞(μ_N)≤ Cλ_i^m+5/2(ρ_N/ζ_N+ζ_N).
  
  Otherwise if λ_g≤λ_i^(N)+1, then Proposition <ref> and (<ref>) implies 
  
      g_L^∞(μ_N)≤ C[λ_i^(N)+1]^m+1g_L^1(μ_N)≤ Cλ_i^m+1g_L^2(μ_N)≤ Cλ_i^m+1i^3/2√(ρ_N/ζ_N+ζ_N√(λ_i)),
  
  where we have used the fact that λ_i^(N)≤ Cλ_i and g_L^1(μ_N)≤g_L^2(μ_N). 
  Comparing  (<ref>) with (<ref>), we see that the error in (<ref>) dominates. To finish, we need again to lift the L^∞(μ_N)-bound to a  L^∞(μ)-bound using regularity of the ψ_i's. Notice that  
  
      ψ_i^(N)∘ T_N-ψ_i_L^∞(μ)   =1≤ j≤ Nmax x∈ U_jsup  |ψ_i^(N)∘ T_N(x)-ψ_i(x)| 
         ≤1≤ j ≤ Nmax x∈ U_jsup (|ψ_i^(N)(X_j)-ψ_i(X_j)|+|ψ_i(X_j)-ψ_i(x)|)
         ≤ψ_i^(N)-ψ_i_L^∞(μ_N)+∇ψ_i_L^∞(μ)ρ_N
         ≲λ_i^m+1i^3/2√(ρ_N/ζ_N+ζ_N√(λ_i)),
  
  where in the last step ∇ψ_i_L^∞(μ)ρ_N is a higher order term that we drop.  
  
  
  
  
  Now we are finally ready to show the L^∞(μ) convergence of  towards W^.
  
  Recall that
  
      =∑_i=1^k_N[1+λ_i^(N)]^-s/2ξ_i ψ_i^(N)∘ T_N, 
      
           W^   =∑_i=1^∞ (1+λ_i)^-s/2ξ_iψ_i.
  
  Consider two intermediate quantities: 
  
      W_n^   =∑_i=1^k_N(1+λ_i)^-s/2ξ_iψ_i^(N)∘ T_N, 
      W_n^   =∑_i=1^k_N(1+λ_i)^-s/2ξ_iψ_i.
  
  By the triangle inequality, 
  
      𝔼_ξ-W^^2_L^∞(μ)   ≤𝔼_ξ(-W_n^_L^∞(μ)+W_n^-W_n^_L^∞(μ)+W_n^-W^_L^∞(μ))^2
         ≤ 3(𝔼_ξ-W_n^^2_L^∞(μ) + 𝔼_ξW_n^-W_n^^2_L^∞(μ)+ 
          𝔼_ξW_n^-W^^2_L^∞(μ)),
  
  so it suffices to bound each term. 
  First, by Proposition <ref> and Weyl's law,
  
      𝔼_ξW_n^-W^_L^∞(μ)^2
             ≤𝔼_ξ[∑_i=k_N+1^∞∑_j=k_N+1^∞(1+λ_i)^-s/2(1+λ_j)^-s/2|ξ_i||ξ_j| ψ_i_L^∞(μ)ψ_j_L^∞(μ)]
         ≲∑_i=k_N+1^∞∑_j=k_N+1^∞(1+λ_i)^-s/2(1+λ_j)^-s/2λ_i^m-1/4λ_j^m-1/4
         =[∑_i=k_N+1^∞(1+λ_i)^-s/2λ_i^m-1/4]^2
          
          ≲[∫ _k_N^∞ x^-2s+m-1/2mdx]^2 
          ≲   k_N^-2s+3m-1/m.
   
  Similarly,
  
      𝔼_ξ-W_n^^2_L^∞(μ)   ≲[∑_i=1^k_N|[1+λ_i^(N)]^-s/2-[1+λ_i]^-s/2|ψ_i^(N)∘ T_N_L^∞(μ)]^2.
  
  By Lipschitz continuity of x^-s/2 on [1,∞) and (<ref>), for i=1,…,k_N,
  
      |[1+λ_i^(N)]^-s/2-[1+λ_i]^-s/2|
             ≤[(1+λ_i^(N))∧ (1+λ_i)]^-s/2-1|λ_i^(N)-λ_i| 
         ≲[(1+λ_i^(N))∧ (1+λ_i)]^-s/2-1λ_i (ρ_N/ζ_N+ζ_N√(λ_i))
         ≲λ_i^-s/2(ρ_N/ζ_N+ζ_N√(λ_i)).
  
  By (<ref>), if we choose k_N so that 
  
      k_N^3(ρ_N/ζ_N+ζ_N√(λ_k_N))≪ 1,
  
  then we have, for i=1,…,k_N,
  
      ψ_i^(N)∘ T_N_L^∞(μ)≤ψ_i^(N)∘ T_N-ψ_i_L^∞(μ)+ψ_i_L^∞(μ)≲λ_i^m-1/4.
  
  Therefore 
  
      𝔼_ξ-W_n^^2_L^∞(μ)≲[∑_i=1^k_Nλ_i^-s/2(ρ_N/ζ_N+ζ_N√(λ_i))λ_i^m-1/4]^2
          ≲(ρ_N/ζ_N+ζ_N√(λ_k_N))^2,
  
  where we have used that λ_i^-2s+m-1/4 is summable for s>3/2m-1/2. 
  Lastly, by Lemma <ref>,
  
      𝔼_ξW_n^-W_n^^2_L^∞(μ)   ≲[∑_i=1^k_N(1+λ_i)^-s/2ψ_i^(N)∘ T_N-ψ_i_L^∞(μ)]^2
         ≲[∑_i=1^k_N (1+λ_i)^-s/2λ_i^m+1i^3/2√(ρ_N/ζ_N+ζ_N√(λ_i))]^2
         ≲[1∨ k_N^-2s+9m+5/m](ρ_N/ζ_N+ζ_N).
  
  Combining (<ref>), (<ref>), (<ref>), we have
  
      𝔼_ξ-W^^2_L^∞(μ)   ≲  k_N^-2s+3m-1/m  + [1∨ k_N^-2s+9m+5/m](ρ_N/ζ_N+ζ_N).
  
  Now we will set ζ_N and k_N based on the dimension m and the smoothness parameter s. 
  
  Case 1: s≤9/2m+5/2. Let δ>0 be arbitrary and set
  
      ζ_N≍ N^-1/m+4+δ(log N)^p_m/2,     k_N≍ N^1/(m+4+δ)(6+6/m)(log N)^-mp_m/12m+12,        m≤ 4,
      ζ_N≍ N^-1/2m(log N)^p_m/2,    k_N≍ N^1/12m+12(log N)^-mp_m/12m+12,        m≥ 5.
  
  One can check that the conditions in (<ref>) and (<ref>) are satisfied and
  we get 
  
      𝔼_ξ-W^^2_L^∞(μ)≲
          N^-2s-3m+1/(m+4+δ)(6m+6)(log N)^p_m(2s-3m+1)/12m+12 ,       m≤ 4,
      
          N^-2s-3m+1/m(12m+12) (log N)^p_m(2s-3m+1)/12m+12,        m≥ 5.
  
  
  Case 2: s>9/2m+5/2.
  Now (<ref>) simplifies to
  
      𝔼_ξ-W^^2_L^∞(μ)≲  k_N^-2s+3m-1/m+(ρ_N/ζ_N+ζ_N).
  
  Therefore by setting
  
      ζ_N≍ N^-1/m+4+δ(log N)^p_m/2,    k_N≍ N^m/(m+4+δ)(2s-3m+1)(log N)^-mp_m/4s-6m+2,        m≤ 4, 
      ζ_N≍ N^-1/2m(log N)^p_m/2,    k_N≍ N^1/4s-6m+2(log N)^-mp_m/4s-6m+2,        m≥ 5,
  
  we have 
  
      𝔼_ξ-W^^2_L^∞(μ)≲
          N^-1/m+4+δ(log N)^p_m/2,        m≤ 4,
      
          N^-1/2m(log N)^p_m/2,        m≥ 5.
  
  Again one can check that the conditions in (<ref>) and (<ref>) are satisfied. This finishes the proof.
  
  
  
  
  
  
  
  
  
   §.§ Putting Everything Together
  
  Now we are ready to prove Theorem <ref>. 
  
  By Subsection <ref>, it suffices to show that 
  
      𝔼_X𝔼_f_0(f∈ℐ(L^∞(μ_N_n)):f-f_0_n≥ M_n_n | _N_n, _n)  0,
  
  with the said _n. 
  Let A_n be the high probability event in Corollary <ref>. Denote F_n(_N_n,_n):= (f∈ℐ(L^∞(μ_N_n)):f-f_0_n≥ M_n_n  | _n). We have 
  
      𝔼_X𝔼_f_0 F_n(_N_n,_n) = 𝔼_X[𝔼_f_0F_n(_N_n,_n)]1_A_n + 𝔼_X[𝔼_f_0F_n(_N_n,_n)]1_A_n^c.
  
  Since F_n(_N_n,_n)≤ 1, the second term is upper bounded by ℙ_X(A_n^c)→ 0. It then suffices to show 𝔼_f_0F_n(_N_n, _n)→ 0 in the event of A_n. By Corollary <ref> we can have 10𝔼_ξ-W^^2_L^∞(μ)≤ n^-1 if we set a large enough proportion constant for N_n. This together with Theorem <ref> and <cit.>[Theorem 2.2] implies that there exists sets  B_n⊂ℐ( L^∞(μ_N_n)) (by taking the intersection of the sets B_n in Theorem <ref> and ℐ( L^∞(μ_N_n)))  so that, for the same _n in the theorem statement,
  
      log N(6_n,B_n,·_L^∞(μ))   ≤  24Cn_n^2, 
      (B_n^c)   ≤ e^-4Cn_n^2,
      (w-Φ^-1(f_0)_L^∞(μ)<4_n)   ≥ e^-4n_n^2,
  
  where =ℒ(). 
  
  Case 1: Regression. In the regression case, since Φ is the identity, we have = and the above three conditions are true with  replaced by . Furthermore, since ·_n is upper bounded by ·_L^∞(μ), the above conditions remain true for the empirical norm. This together with the general results in <cit.>[Section 8.3.2] proves the assertion. 
  
  Case 2: Classification. In the classification setting, by <cit.>[Lemma 2.8] the average Kullback-Leibler divergence and variance in (<ref>) between the densities after applying the link function Φ are upper bounded by a multiple of the empirical norm. In particular if we set ℬ_n=Φ(B_n), then the conditions in (<ref>), (<ref>) and (<ref>) hold, for a possibly different proportion constant for _n. Moreover, as discussed before, the root average square Hellinger distance d_n in (<ref>) is equivalent to the empirical norm. Hence the result follows by <cit.>[Theorem 4].
  
  
  
  
  § DISCUSSION
  
  In this paper we have analyzed graph-based SSL using recent results on spectral convergence of graph Laplacians and standard Bayesian nonparametrics techniques.  We show that, for a suitable choice of prior constructed with sufficiently many unlabeled data, the posterior contracts around the truth at a rate that is minimax optimal up to logarithmic factor. Our theory applies to both regression and classification. 
  
  We have assumed throughout that the X_i's are uniformly distributed on . Our results can be generalized to nonuniform positive density q with respect to the volume form. In such a case, the continuum field is the Gaussian measure 𝒩(0, (I-Δ_q))^-s where  -Δ_q:=-1/qdiv (q^2 ∇ ) is a weighted Laplacian-Beltrami operator, for which spectral convergence results can be found in <cit.>. Since we do not have an explicit dependence of the proportion constant on the index i as in Proposition <ref>, we have chosen to present our result in the uniform case. However, similar conclusions should be expected to hold in the general case.
  
  The Bayesian methodology we analyzed is inspired by popular existing graph-based optimization methods for SSL <cit.>. In order to achieve optimal contraction, the prior smoothness parameter needs to match the regularity of the truth, which is rarely available in applications. An important research direction stemming from our work is the development and analysis of adaptive Bayesian SSL methodologies that can achieve optimal contraction without a priori smoothness information.  We expect that the existing results on adaptive estimation on manifolds <cit.> and on large graphs <cit.> will be an important stepping stone in this direction.
  
  
  
  
  
  
  
  § ACKNOWLEDGMENT
  
  Both authors are thankful for the support of NSF and NGA through the grant DMS-2027056. 
  The work of DSA was also partially supported by the NSF Grant DMS-1912818/1912802. 
  abbrvnat 
  
  
  