---
abstract: 'While spectral embedding is a widely applied dimension reduction technique in various fields, so far it is still challenging to make it scalable and robust to handle “big data”.  Motivated by the need of handling such data, we propose a novel spectral embedding algorithm, which we coined [*Robust and Scalable Embedding via Landmark Diffusion*]{} (ROSELAND).  In short, we measure the affinity between two points via a set of landmarks, which is composed of a small number of points, and “diffuse” on the dataset via the landmark set to achieve a spectral embedding.  The embedding is not only scalable and robust, but also preserves the geometric properties under the manifold setup.  The Roseland can be viewed as a generalization of the commonly applied spectral embedding algorithm, the [*diffusion map*]{} (DM), in the sense that it shares various properties of the DM.  In addition to providing a theoretical justification of the Roseland under the manifold setup, including handling the U-statistics like quantities and providing a spectral convergence rate, we show various numerical simulations and compare the Roseland with other existing algorithms. ** Keywords**: graph Laplacian, diffusion maps, Nyström, landmark, scalability, robustness, Roseland'
address:
- 'Department of Mathematics, Duke University, Durham, NC, USA'
- 'Department of Mathematics and Department of Statistical Science, Duke University, Durham, NC, USA; Mathematics Division, National Center for Theoretical Sciences, Taipei, Taiwan'
author:
- Chao Shen
- 'Hau-Tieng Wu'
bibliography:
- 'Roseland-20191230.bib'
title: 'Scalability and robustness of spectral embedding: landmark diffusion is all you need'
---

Introduction {#Sect:Introduction}
============

Unsupervised learning is arguably the holy grail in the field of artificial intelligence, and it is arguably that the more data we have, the better we can learn.  So far there have been many unsupervised learning algorithms proposed and it is still an active studying field.  In general, those algorithms share a common ground.  The learner designer constructs an optimization framework that captures the intended properties of the learning process, and then designs an algorithm to solve the optimization problem.  Based on the nature of an algorithm, it can be roughly classified into two classes – spectral or not.  Spectral algorithms include ISOMAP [@tenenbaum2000global], locally linear embedding (LLE) [@roweis2000nonlinear], Hessian LLE [@Donoho_Grimes:2003], eigenmap [@belkin2003laplacian], diffusion map (DM) [@Coifman2006], vector DM (VDM) [@singer2012vector], to name but a few.  Those algorithms have been widely applied to various scientific fields and various theoretical foundations have been established to support those algorithms in the past decades.  Under the manifold setup, we have had a rich knowledge about the geometric and asymptotic behavior of those algorithms.  For example, the DM and VDM are both based on the diffusion process [@Coifman2006; @singer2012vector], and asymptotically they converge to the Laplace-Beltrami operator or connection Laplacian so that the spectral geometry theory can be applied; the LLE algorithm is not diffusion-based and the underlying kernel is asymmetric and depends on the geometry of the dataset [@Wu_Wu:2017].  However, there are still various open problems remain toward a better unsupervised learning framework.  One critical challenge is how to make an algorithm scalable, which is always a numerical challenge and is critical in this “big data era”.  Another critical challenge is how to handle the inevitable noise, particularly when the noise is large and high dimensional.  How to handle these two challenges together is the focus of this paper. 

Take the DM algorithm as an example.  The DM is based on the eigendecomposition of the graph Laplacian (GL) matrix, which is constructed by determining the affinity among every pair of points in the database.  The algorithm has been shown to perform well and “robustly” when the database is “tiny”, like in the order of $10^3\sim 10^4$. However, when the database gets larger, like in the order of $10^6$ or above, the algorithm needs a modification.  Specifically, if the GL is dense, a full eigendecomposition is not feasible, and the usually applied trick like the k-nearest neighbor (kNN) scheme is not robust to noise.  Specifically, when the dataset is noisy and the neighboring information is [*not*]{} provided, obtaining a reliable kNN information is challenging, and the geometric information cannot be captured.  Another practical solution is directly subsampling the dataset, and then recovering the information of interest by the [*Nyström extension*]{} [@czaja2017overview; @belabbas2009landmark; @williams2001using].  This approach is also called the [*Nyström low-rank approximation*]{} [@chang2013asymptotic], the [*kernel extension method*]{} [@fowlkes2004spectral], or in general the [*interpolative decomposition*]{} [@martinsson2011randomized].  This approach has various theoretical backups, for example [@chang2013asymptotic], and has been widely applied.  While it works well for some missions, this approach is limited by the information loss during the subsampling process; for example, the geometric information is lost when the dataset is sampled from a manifold.  Moreover, to the best of our knowledge, how does it perform when combined with the above-mentioned spectral based unsupervised learning algorithms is not yet well explored, not to mention how it impacts the algorithm under the manifold setup, or when the data is noisy and nonuniform sampled.  Yet another approach is speeding up the matrix decomposition by taking randomization into account.  For example, we can construct a thin matrix by taking a random subset of columns of the GL matrix and speed up the algorithm by taking the singular value decomposition (SVD) into account.  For example, [@martinsson2011randomized] provides an efficient algorithm to approximate the SVD.  While this approach has been widely applied, to the best of our knowledge, we have limited knowledge about how it helps the spectral embedding algorithms, and how robust it is to the inevitable noise.  Unlike the above, in this paper we propose a novel algorithm that resolves these two challenges, robustness and scalability, as well as geometry recovery guarantee.  In short, we measure the affinity between any pair of two points via a [*landmark set*]{}, which is composed of a small number of points.  With this affinity, the embedding is carried out via the spectral decomposition.  As we will make clear soon, this algorithm has an interpretable geometric meaning, and is directly related to the diffusion process.  Hence, we coin the proposed algorithm the [*RObust and Scalable Embedding via LANdmark Diffusion*]{} (Roseland).  Moreover, we provide a series of theoretical supports for the Roseland, like the spectral convergence of eigenvalues/eigenvectors of the Roseland to the eigenvalues/eigenfunctions of the Laplace-Beltrami operator in the $L^\infty$ sense.  Based on the theoretical results, we propose a [*design-based*]{} landmark set sampling scheme to handle the inevitable non-uniform sampling in the real world data.  This is particularly useful when we want to recover the intrinsic Laplace-Beltrami operator of the manifold. 

Related work – scalability
--------------------------

To better position our contribution, we summarize various related work in this section.  The review paper by Czaja, et.al. [ @czaja2017overview] contains a comprehensive categorization of numerical acceleration techniques for nonlinear dimension reduction.  The acceleration algorithms are roughly classified into three categories.  In the first category, the dataset is compressed so that the relationship between pairs of points are well preserved.  For example, the random projection can be applied, or the dataset is converted to a well-designed basis under the compressed sensing framework.  This step can save us a little bit of time when computing pairwise distances used for neighborhood search.  In the second category, we may try to accelerate the kNN search step.  A brutal force method for computing the exact kNN graph requires $\Theta(n^{2})$. Many faster algorithms, deterministic or randomized, exact or approximate, have been developed in the past decades.  In the third category, we may accelerate the eigen-decomposition step.  For example, the kernel decomposition is approximated by classical iteration-based algorithms and the matrix decomposition can be evaluated by randomized algorithm, where rigorous analyses have been developed [@halko2011finding; @rokhlin2009randomized].  However, it is indicated in [@czaja2017overview] that the error bounds are usually pessimistic when compared with results of numerical experiments.  A summary with citations of methods in each category can be found in [@czaja2017overview]. 

We now review various algorithms that are directly related to our work.  The closest algorithm to the Roseland is the one introduced in [@haddad2014texture] to handle the texture separation problem.  We call this algorithm the [*HKC*]{}, which stands for initials of three authors in [@haddad2014texture].  The authors first convert an image into a collection of small patches, and choose a collection of specific patterns of interest as the [*reference set*]{}.  Then one can construct an affinity matrix associated with the set of patches, where the affinity between patches are based on the landmark set.  However, the normalization in the HKC is different from the Roseland, and this difference turns out to be significant.  Moreover, it is not clear how does the HKC performs under the manifold setup.  The HKC can be classified as the third category. 

Another directly related algorithm is the common Nyström extension [@belabbas2009landmark; @fowlkes2004spectral; @williams2001using].  We run eigen-decompostion on a small subset of the whole database, and then extend the eigenvectors to the whole dataset.  This Nyström extension approach can be classified as the third category.  There are various extensions or refinements of the Nyström extension method, for example, [@bermanis2013multiscale].  From the theoretical perspective, to the best of our knowledge, Chang, Bai, et.al [@chang2013asymptotic] is the only existing literature in this field.  The authors analyzed the asymptotic spectral error bounds between the ground truth spectrum of the kernel function, full kernel matrix and the Nyström low-rank approximation of the full kernel matrix. 

Yet another and fundamental approach is designing a better to speeding up the basic eigendecomposition or SVD themselves.  But this direction is out of the scope of our work.  There are some closely related but different algorithms in the field, for example, CUR decomposition [@mahoney2009cur], “UBV” decomposition [@cheng2005compression], or some studies focusing solely on accelerating the spectral clustering.  For the readers’ convenience, we summarize them in Appendix \[Appendix:MoreRelatedWork\].  To the best of our knowledge, none of the above-mentioned work, except [@chang2013asymptotic], provides theoretical analysis to answer questions like what is the [*asymptotical behavior*]{} of the algorithm?  However, even in [@chang2013asymptotic], it is not clear how much geometric information is lost.  In general, if we model the nonlinearity of the dataset by a manifold, we would like to know if we still have a convergence to the Laplace-Beltrami operator.  Without these understanding, we cannot answer questions like how to choose, or even design, landmark set so that the performance is guaranteed to some extent. 

Related work – robustness {#section: related work: robustness}
-------------------------

Compared with the scalability, there are less studies focusing on the robustness issue.  One intuitive idea is “denoising” the dataset before applying the spectral embedding algorithm.  However, it is in general a challenging problem since we usually do not know the structure of the dataset, and extracting the structure of the dataset is the main target.  Under the manifold setup, it has been proposed to denoise the dataset via the “reverse diffusion” scheme [@NIPS2006_2997].  While it might work for practical problems, the algorithm might not be scalable.  A commonly encountered issue is the stability of the kNN scheme.  When the neighboring information is not provided, usually we apply the kNN to construct the affinity graph on the dataset.  However, it is well known that finding neighbors via the kNN is noise-sensitive unless the pairwise distance is robust to noise.  Usually, unless the data point has extra structure so that a robust metric can be applied, for example, in the image analysis [@cheng2009learning], it is challenging to achieve a robust pairwise distance.  Some authors propose to take the tangent plane structure to determine neighbors [@wang2005adaptive]; however, it is well known that in the high dimensional setup, the tangent space estimation via the principle component analysis is biased [@johnstone2007high] and the benefit might be limited in the practical setup.  Another approach is taking the label into account to improve the stability of the kNN scheme [@ROHBAN20121363], but this approach is out of the scope of this work.  In short, while the kNN is also an acceleration tool (the second category of acceleration [@czaja2017overview]), it is only useful when the neighboring information can be accurately estimated.  When the edge information is known, in [@Steinerburger2016], the author proposes to design a self-consistency Markov chain before the spectral embedding by modifying the non-lazy random walk via diffusion.  In general, to the best of knowledge, the robustness of the GL-based algorithms was first studied in [@ElKaroui:2010a] under the random matrix framework.  The result was later extended to handle a large noise setup [@el2016graph], where the authors suggest to take a complete non-lazy random walk to stabilize the spectral embedding methods.  Obviously this approach is not scalable.  In that paper, it is also shown that the number of nearest neighbors should be chosen large enough to stabilize the kNN scheme. 

Our contribution
----------------

We propose the Roseland algorithm that simultaneously handle the above-mentioned two challenges – scalability and robustness.  The algorithm is intuitive and can be summarized in three steps.  First, we find a “small” subset of points from the whole dataset, either randomly or by design, or collect a separate clean point cloud of small size, which we call a [*landmark set*]{}.  Second, we construct an affinity matrix recording the affinities between points in the whole dataset and the landmark set, and normalize it properly.  This normalized affinity matrix is thin; that is, there are fewer columns then rows.  Third, evaluate the singular vectors and singular values of the normalized affinity matrix, and embed the dataset with the singular vectors and singular values. 

From the theoretical perspective, we provide a series of theoretical supports when the point cloud is sampled from a manifold, including the pointwise convergence and the spectral convergence, to show that the Roseland recovers the Laplace-Beltrami operator of the underlying manifold supporting the data set.  First, we observed a peculiar kernel behavior of the Roseland; specifically, the “effective kernel” associated with the Roseland is not fixed but adaptive to the chosen landmark set.  This is different from the ordinary kernel method, where the applied kernel is universal.  Second, we handle the dependence coming from the landmark set.  This dependence forbids us from using the standard large deviation theory to study the variance, and we count on the tools originated from the U-statistics to handle this challenge.  As a result, we provide not only the spectral convergence of the eigenvalues/eigenvectors of the Roseland to the eigenvalues/eigenfunctions of the Laplace-Beltrami operator, but also the convergence rate.  We argue that the convergence rate is controlled by the size of the landmark set.  Based on the analysis, we provide a natural diffusion interpretation of the role of selecting the landmark set in the Roseland, which allows us to link it to the widely applied algorithm, the diffusion map (DM).  A direct theoretical consequence is a guarantee that the diffusion distance (DD) approximates the underlying geodesic distance, and we can design a subsampling scheme to handle the non-uniform sampling, both of which are critical in practical applications.  From the application perspective, we provide a series of numerical simulations comparing the Roseland with other algorithms, including the Nyström extension and the HKC, under various situations.  We observe from simulations that the Roseland is more robust to noise than its competitors, which suggests an potential application of using the Roseland to “denoise” the original data. 

Organization of the paper
-------------------------

In Section \[algs\] we recall the DM algorithm and introduce the Roseland.  In Section \[setup\], we provide the mathematical set up of the problem.  In Sections \[Section:PointConvergence\] and \[Section:spectral convergence\], we state our main theoretical results.  In Section \[numerical\], we provide numerical results and analysis.  In Section \[discuss\], discussion and conclusion are provided. 

The proposed Roseland algorithm and relevant algorithms {#algs}
=======================================================

In this section, we assume we have a data set or point cloud $\mathcal{X}=\{x_{i}\}_{i=1}^{n}\subseteq\mathbb{R}^{q}$. Take a set $\mathcal{Y}=\{y_{k}\}_{k=1}^{m}$, which might or might not be a subset of $\mathcal{X}$. We call $\mathcal{Y}$ the [*landmark set*]{}.  Fix a non-negative kernel function $K:\mathbb{R}_+\to \mathbb{R}_+$ with proper decay and regularity; for example, a Gaussian function $K(t)=\frac{1}{\sqrt{2\pi}}e^{-t^2/2}$.

The proposed algorithm – Roseland
---------------------------------

We now introduce the Roseland.  Construct a [*landmark-set affinity matrix*]{} $W^{(\textup{r})}\in\mathbb{R}^{n\times m}$, which is defined as $$\label{landmark affinity matrix in Roseland}
W^{(\textup{r})}_{ik} = K_{\epsilon}(x_{i},y_{k}) := K\left(\frac{\|x_{i}-y_{k}\|_{\mathbb{R}^{q}}}{\sqrt{\epsilon}}\right)\,.$$ That is, the $(i,k)$-th entry of $W^{(\textup{r})}$ is the affinity between $x_i\in \mathcal{X}$ and $y_k\in \mathcal{Y}$. Construct a diagonal matrix $D^{(\textup{r})}$ as $$D^{(\textup{R})}_{ii}:=e_i^\top W^{(\textup{r})}(W^{(\textup{r})})^{\top}\boldsymbol{1}\,,$$ where $\boldsymbol{1}$ is a $n\times 1$ vector with all entries $1$, and $e_i$ is the unit vector with $1$ in the $i$-th entry.  With $W^{(\textup{r})}$ and $D^{(\textup{R})}$, we evaluate the SVD of $(D^{(\textup{R})})^{-1/2}W^{(\textup{r})}$: $$\label{DR-1/2Wr SVD in Roseland}
(D^{(\textup{R})})^{-1/2}W^{(\textup{r})}=U\Lambda V^\top,$$ where the singular values are $\sigma_1\geq \sigma_2\geq\ldots\geq\sigma_m\geq 0$. Set $\bar{U}:=(D^{(\textup{R})})^{-1}U$, and set $\bar{U}_m\in \mathbb{R}^{n\times m}$ to be a matrix consisting of the second to the $(m+1)$-th left singular vectors.  Also set $L_m:=\texttt{diag}(\sigma^2_2,\ldots,\sigma^2_{m+1})$. The Roseland embedding is defined by $$\label{embedding in Roseland}
\Phi^{(\textup{R})}_t:x_{i}\mapsto e_i^\top \bar{U}_{m}L_{m}^t\,,$$ where $t>0$ is the chosen diffusion time.  With the Roseland embedding, we have the associated Roseland diffusion distance (RDD) $$\label{DD in Roseland}
D^{(\textup{R})}_t(x_i,x_j):=\|\Phi^{(\textup{R})}_t(x_{i})-\Phi^{(\textup{R})}_t(x_j)\|_{\mathbb{R}^{m}}\,.$$ The Roseland algorithm is summarized in Algorithm \[refdmalg\]. 

![ Main idea of Roseland: to measure the similarity between $x_{1}$ to $x_{2}$, instead of diffuse from $x_{1}$ to $x_{2}$ directly, we take a detour and first diffuse $x_{1}$ to the landmarks $y_{1},y_{2},y_{3}$, and then diffuse from the landmarks to $x_{2}$.\[Figure:refDMillustration\]](ref_illu.pdf){width="75.00000%"}

Construct the affinity matrix $W^{(\textup{r})}$ w.r.t.  the landmark set.  Construct the degree matrix $D^{(\textup{R})}$. Run SVD $(D^{(\textup{R})})^{-1/2}W^{(\textup{r})} = U\Lambda V^{\top}$, where $U\in O(n)$ and $V\in O(m)$. Denote singular values as $\sigma_1\geq\sigma_2\geq\ldots\geq 0$. Set $\bar{U} = (D^{(\textup{R})})^{-1/2}U$. Let $\bar{U}d\in\mathbb{R}^{n\times q'}$ be the second to the $(q'+1)$-th columns of $\bar U$. Set $L_{q'}:=\texttt{diag}(\sigma_2^2,\sigma_3^2,\ldots,\sigma_{q'+1}^2)$. Embed $\mathcal{X}$ via $\Phi^{(\textup{R})}_t:x_{i}\mapsto e_i^\top \bar{U}_{q'}L_{q'}^t$ Output $\{\Phi^{(\textup{R})}_t(x_i)\}_{i=1}^n\subset \mathbb{R}^{q'}$ and/or the associated RDD. 

We now take a closer look at the proposed Roseland.  Given $W^{(\textup{r})}\in \mathbb{R}^{n\times m}$, construct a new matrix $$W^{(\textup{R})} := W^{(\textup{r})}(W^{(\textup{r})})^{\top}\in \mathbb{R}^{n\times n}\,,$$ which can be viewed as a new affinity matrix.  Indeed, since the kernel is chosen to be a non-negative function, $W^{(\textup{r})}$ is a matrix with non-negative entries, and so is $W^{(\textup{R})}$. Therefore, we can view $W^{(\textup{R})}$ as an affinity matrix defined on $\mathcal{X}$, where the affinity between $x_{i}$, $x_{j}$ via the landmark set $\mathcal{Y}$ is $$\label{definition WR in Roseland}
W^{(\textup{R})}_{ij}=\sum_{k=1}^{m}K_{\epsilon}(x_{i},y_{k})K_{\epsilon}(y_{k},x_{j}).$$ We call it the [*landmark-affinity matrix*]{}.  We mention that unlike the traditional affinity matrix, in general we cannot find a fixed kernel $\bar{K}$ and a bandwidth $\bar \epsilon$ so that $W^{(\textup{R})}_{ij}= \bar K_{\bar \epsilon}(x_i,x_j)$ for $i,j=1,\ldots,n$. Later, we will see how this “new” kernel function depends on $\mathcal{Y}$. Thus, by construction, $D^{(\textup{R})}$ is nothing but the degree matrix associated with the landmark-affinity matrix $W^{(\textup{R})}$. Clearly, $$\label{transition matrix definition of Roseland}
A^{(\textup{R})}:=(D^{(\textup{R})})^{-1}W^{(\textup{R})}$$ is a transition matrix on $\mathcal{X}$. In this sense, we have a Markov process, or diffusion, on $\mathcal{X}$, where if we want to diffuse from $x_i$ to $x_j$, we always go through $\mathcal{Y}$. Moreover, note that we have $$\begin{aligned}
&(D^{(\textup{R})})^{-1/2}W^{(\textup{r})}[(D^{(\textup{R})})^{-1/2}W^{(\textup{r})}]^\top
=(D^{(\textup{R})})^{-1/2}W^{(\textup{R})}(D^{(\textup{R})})^{-1/2}\,,\end{aligned}$$ where the right hand side is symmetric.  If $(D^{(\textup{R})})^{-1/2}W^{(\textup{R})}(D^{(\textup{R})})^{-1/2}$ is non-negative definite, the SVD of $(D^{(\textup{R})})^{-1/2}W^{(\textup{r})}$ recovers the eigen-structure of the non-negative definite matrix $(D^{(\textup{R})})^{-1/2}W^{(\textup{R})}(D^{(\textup{R})})^{-1/2}$.

Graph Laplacian and Diffusion Maps
----------------------------------

We now recall the GL, and see the similarity between the Roseland and a well-known algorithm, the diffusion map (DM).  Define the affinity matrix $W\in \mathbb{R}^{n\times n}$ by $$\label{Definition affinity matrix W}
W_{ij} := 
  K\left(\frac{\|x_{i}-x_{j}\|_{\mathbb{R}^{q}}}{\sqrt{\epsilon}}\right)\,$$ and define the corresponding [*degree matrix*]{} $D\in \mathbb{R}^{n\times n}$, which is a diagonal matrix defined as $$D_{ii}:=\sum_{j=1}^nW_{ij}.$$ In some cases, we may consider an $\alpha$-normalization [@Coifman2006] to refine the affinity, where $\alpha\in[0,1]$. The [*$\alpha$-normalized*]{} affinity matrix $W^{(\alpha)}\in \mathbb{R}^{n\times n}$ is defined as $$W_{ij}^{(\alpha)}:=\frac{W_{ij}}{D_{ii}^{\alpha}D_{ii}^{\alpha}}\,,$$ where $W_{ij}^{(\alpha)}$ is the $\alpha$-normalized affinity between $x_i$ and $x_j$, Note that $W^{(0)}=W$ defined in .  With the affinity matrix $W^{(\alpha)}$, we define an associated degree matrix $D^{(\alpha)}\in \mathbb{R}^{n\times n}$, which is defined as $$D_{ii}^{(\alpha)}:= \sum_{j=1}^nW_{ij}^{(\alpha)}\,.$$ The GL is defined as $L^{(\alpha)}:=I-A$, where $$\label{Eq A ordinary DM}
A:=(D^{(\alpha)})^{-1}W^{(\alpha)}$$ is the associated transition matrix.  As a result, we have a Markov process on the dataset $\mathcal{X}$. Note that the landmark-affinity matrix $W^{(\textup{R})}$ defined in , and the associated transition matrix $(D^{(\textup{R})})^{-1}W^{(\textup{R})}$ can be viewed as a different way of constructing a Markov process on the dataset $\mathcal{X}$. With the above preparation, we run the eigen-decomposition of the transition matrix In general, $A$ is not symmetric but is similar to the symmetric matrix $$\label{Eq bar A ordinary DM}
\bar A := (D^{(\alpha)})^{-1/2}W^{(\alpha)}(D^{(\alpha)})^{-1/2}\,,$$ which we have the eigen-decomposition $\bar A = U\Lambda U^{\top}$, where $U\in O(n)$ and $\Lambda=\text{diag}(\lambda_1,\lambda_2,\ldots,\lambda_n)$, where the eigenvalues are ordered by $\lambda_1\geq\lambda_2\geq\ldots\geq\lambda_n$. The right eigenvectors of $A$ can be recovered by $\tilde{U} = (D^{(\alpha)})^{-1/2}U$. Denote $\phi_i$ to be the $i$-th right eigenvector of $A$.  Note that the decomposition in is a parallel step of this eigendecomposition. 

With the spectral decomposition of the GL, we can carry out various missions, like the spectral clustering or the spectral embedding.  Among various algorithms, we show the DM, which is a famous algorithm that we will see to be closely related to the Roseland.  For a chosen [*diffusion time $t$*]{}, the DM embeds $\mathcal{X}$ via the map $$\Phi_t:x_{i}\mapsto (\lambda_2^t\phi_2(i),\ldots,\lambda_{q'+1}^t\phi_{q'+1}(i))\in \mathbb{R}^{q'},$$ where $q'$ is the dimension chosen by the user.  With the DM, the diffusion distance (DD) with the diffusion time $t>0$ is defined as $$D_t(x_i,x_j):=\|\Phi_t(x_{i})-\Phi_t(x_j)\|_{\mathbb{R}^{q'}}\,.$$ Clearly, the Roseland embedding and the RDD are closely related to the DM and the DD. 

In practice, to reduce the spatial and computational complexity of the DM, it is common to use the kNN scheme to determine the affinity matrix; that is, replacing by $$\label{Definition affinity matrix W kNN}
W_{ij} :=\left\{ \begin{array}{ll} 
  K\left(\frac{\|x_{i}-x_{j}\|_{\mathbb{R}^{q}}}{\sqrt{\epsilon}}\right) & \textup{if}\; x_{j}\;\textup{is the first $k$-th neighbor of}\;x_i \\
  0 & \textup{otherwise}\,,
\end{array}\right.$$ where $k$ is pre-determined by the user.  The kNN scheme is related to, but different from, choosing a compactly supported kernel.  For example, if $K(t)$ is $1$ when $t\in [0,1]$ and $0$ when $t>1$, the $(i,j)$-th entry of the affinity matrix is non-zero only if the distance between $x_i$ and $x_j$ is less than $\sqrt{\epsilon}$. In general, the affinity matrix constructed by the kNN scheme is not symmetric, while that constructed by the compactly supported kernel is symmetric.  We refer the reader with interest in the relationship between the kNN scheme and the compactly support kernel approach to, for example, [@Wu_Wu:2017 Section 5]. 

The HKC algorithm
-----------------

We now summarize the HKC algorithm [@haddad2014texture].  First, form the affinity matrix between the data set and the landmark set just like in the Roseland, $$W^{(\textup{HKC})}_{ik}= K_{\epsilon}(x_{i},y_{k})\,,$$ where $i=1,\ldots,n$ and $k=1,\ldots,m$; that is, $W^{(\textup{HKC})}=W^{(\textup{r})}$. Note that in [@haddad2014texture] the authors call the landmark set the reference set.  Then, unlike the Roseland, we compute a $n\times n$ diagonal matrix by $$D^{(\textup{HKC})}_{ii}=\sum_{j=1}^m W^{(\textup{HKC})}_{i,j}\,,$$ where $i=1,\ldots,n$. Then, convert $W^{(\textup{HKC})}$ to be row stochastic by: $$\begin{aligned}
\label{HKC affinity}
A^{(\textup{HKC})} = (D^{(\textup{HKC})})^{-1} W^{(\textup{HKC})}\in \mathbb{R}^{n\times m}\,.\end{aligned}$$ Clearly, $A^{(\textup{HKC})}_{ij}$ describes how we diffuse from $y_j$ to $x_i$. The goal is to embed the data via the eigenvectors $\psi_{j}$ of the matrix $\bar W^{(\textup{HKC})} = A^{(\textup{HKC})}(A^{(\textup{HKC})})^{\top}\in\mathbb{R}^{n\times n}$, which can be computed efficiently by $$\begin{aligned}
\label{HKC extend}
\psi^{(\textup{HKC})}_{j}=\frac{1}{(\lambda^{(\textup{HKC})}_{j})^{1/2}}A^{(\textup{HKC})}\phi^{(\textup{HKC})}_{j}\end{aligned}$$ where $\phi^{(\textup{HKC})}_{j}$ is the $j$-th eigenvectors of the matrix $\tilde W^{(\textup{HKC})}:=(A^{(\textup{HKC})})^{\top}A^{(\textup{HKC})}\in\mathbb{R}^{m\times m}$ associated with the eigenvalue $\lambda^{(\textup{HKC})}_{j}$. The denominator $(\lambda^{(\textup{HKC})}_{j})^{1/2}$ is to ensure that $\|\psi^{(\textup{HKC})}_{j}\|_{2}=\|\phi^{(\textup{HKC})}_{j}\|_{2}$.

In case of non-uniform sampling, we have the following slightly modified version of the algorithm [@haddad2014texture].  First, compute $A^{(\textup{HKC})}$ as in (\[HKC affinity\]) and $\tilde W^{(\textup{HKC})}=(A^{(\textup{HKC})})^{\top}A^{(\textup{HKC})}$. Next, obtain the degree matrix of $W^{(\textup{HKC})}$, denoted by $D^{(\textup{HKC})}_{1}$. Let $A^{(\textup{HKC})}_{1}:=A^{(\textup{HKC})}(D^{(\textup{HKC})}_{1})^{-1}$. Then, make $A^{(\textup{HKC})}_{1}$ row stochastic as above, denoted by $\tilde{A}^{(\textup{HKC})}_{1}.$ Finally, embed the data via the eigenvectors $\psi_{j}$ of the matrix $\bar W^{(\textup{HKC})}_1 = \tilde{A}^{(\textup{HKC})}_{1}(\tilde{A}^{(\textup{HKC})}_{1})^{\top},$ which can be computed efficiently similar to (\[HKC extend\]).  Note that this step is similar to the $1$-normalization step in the DM when the sampling is non-uniform. 

In summary, we see that the HKC algorithm is close to the Roseland, while the normalization is different.  Indeed, due to the normalization , $\bar W^{(\textup{HKC})}$ can only be viewed as a new affinity matrix on $\mathcal{X}$, and is different from $(D^{(\textup{R})})^{-1/2}W^{(\textup{R})}(D^{(\textup{R})})^{-1/2}$. As we will see later, this different normalization step plays a significant role. 

Nyström Extension
-----------------

The direct competitor of the proposed Roseland algorithm is the well known Nyström extension.  This idea has been used, for example, in [@belabbas2009landmark; @fowlkes2004spectral; @williams2001using; @belabbas2009spectral; @coifman2006geometric].  In this algorithm, we first find a subset of the given dataset, run the eigendecomposition on it, and then extend the eigenvectors to the whole dataset.  A direct application is extending the eigenvectors to a new arrival data point.  To be more precise, take a kernel matrix $M\in\mathbb{R}^{N\times N}$, which is assumed to be positive semi-definite, and choose a subset $I\subset \{1,\ldots,N\}$ of size $L\ll N$. To simplify the terminology, we also call $I$ the [*landmark set*]{}.  Rewrite $M$ as: $$M = 
\begin{bmatrix} 
M_{L} & B \\
B^{\top} & C 
\end{bmatrix}\label{segmentation of matrix M}$$ where $M_{L}$ is the submatrix formed by taking the rows and columns corresponding to the landmarks.  Since $L\ll N$, the eigen-decomposition $M_{L}=U_{L}\Lambda_L U^{\top}_{L}$, where $U_L\in O(L)$ and $\Lambda_L$ is diagonal and contains eigenvalues of $M_L$, is much faster than doing the eigen-decomposition of $M$. Furthermore, if we denote $$\widetilde{\Lambda}=\begin{bmatrix}\Lambda_L & 0\\ 0&0\end{bmatrix}\,\,\mbox{ and }\,\,\widetilde{U}=
\begin{bmatrix} 
U_{L} \\
B^{\top}U_{L}\Lambda_L^{-1}
\end{bmatrix}\,,$$ where the eigenvectors are extended to the whole dataset via padding $B^{\top}U_{L}\Lambda_L^{-1}$, we have $$\begin{aligned}
    \widetilde{M} \coloneqq&\, \widetilde{U}\widetilde\Lambda\widetilde{U}^{\top}
     =
    \begin{bmatrix} 
        U_{L} \\
        B^{\top}U_{L}\Lambda_{L}^{-1}
    \end{bmatrix}
    \Lambda_{L}
    \begin{bmatrix} 
        U^{\top}_{L} & \Lambda_{L}^{-1}U^{\top}_{L}B
    \end{bmatrix}
    \\\nonumber =&\,
    \begin{bmatrix} 
        U_{L}\Lambda_{L} U^{\top}_{L} & B\\
        B^{\top} & B^{\top}U_{L}\Lambda_{L}^{-1}U^{\top}_{L}B
    \end{bmatrix}
 =
    \begin{bmatrix} 
        M_{L} & B\\
        B^{\top} & B^{\top}M_{L}^{-1}B
    \end{bmatrix}\,.\end{aligned}$$ The quality of the approximation is usually analyzed by measuring how well $B^{\top}M_L^{-1}B$ approximates $C$ using Frobenius Norm or 2-norm.  Note that this approach is also called the Nyström low-rank approximation since we have $$\begin{aligned}
    \widetilde{M} =
    \begin{bmatrix} 
        M_{L}\\
        B^{\top}
    \end{bmatrix} 
    \begin{bmatrix} 
        M_{L}
    \end{bmatrix}^{-1}
    \begin{bmatrix} 
        M_{L} & B
    \end{bmatrix}.\end{aligned}$$ The asymptotic spectral behavior and error bounds between the ground truth spectrum of the kernel function, full kernel matrix and the Nyström low-rank approximation is analyzed in [@chang2013asymptotic], and we refer the reader with interest there. 

For the spectral embedding purpose, like the Roseland or DM, we are interested in the eigen-decomposition of the transition matrix defined in or .  Since the transition matrix in general is not symmetric, if we want to apply the Nyström extension, we need to consider its symmetric counterpart defined in .  To simplify notations, without loss of generality we consider the symmetric kernel matrix $M=D^{-1/2}WD^{-1/2}$. Note that if we want to directly apply the Nyström extension, we need to compute the affinity matrix $W$ and the degree matrix $D$, which is expensive and the kNN scheme is needed.  Here, to have a fair comparison, we follow the existing literature [@lafon2006data; @singer2012vector] and apply the following slightly modified Nyström extension. 

Suppose we have $n$ data points.  First, run the DM on a subset of the database, also called the landmark set, which is of size $L=n^\beta$, where $\beta\in (0, 1)$. Denote the affinity matrix associated with this landmark set as $W_{L}$, and run the eigen-decomposition of the matrix $D^{-1/2}_{L}W_{L}D^{-1/2} = U_{L}\Lambda_{L}U_{L}^{\top}$, where $D_{L}$ is the degree matrix associated with $W_{L}$. Let $\widetilde{U}_{L}=D_{L}^{-1/2}U_{L}$ be the eigenvectors of $D^{-1}_{L}W_{L}$, we then extend it to the rest $n-L$ points by: $$\widetilde{U}_{\textup{ext}} = D^{-1}_{n-L}E^{\top}\widetilde{U}_{L}\Lambda_{L}^{-1}\in \mathbb{R}^{(n-L)\times L}\,,$$ where $E\in \mathbb{R}^{(n-L)\times L}$ contains the affinities between the landmark set and the remaining points; that is, $E_{i,j}$ is the affinity between $x_j$ in the landmark set and $x_i$ in the remaining dataset, $D_{n-L}$ is a $(n-L)\times (n-L)$ diagonal matrix such that $D_{n-L}(i,i) = \sum_{j=1}^{L}E^{\top}(i,j)$. Put everything together, the eigenvectors to be used in DM on whole dataset is: $$\widetilde{U} = 
\begin{bmatrix} 
\widetilde{U}_{L}\\
\widetilde{U}_{\textup{ext}}
\end{bmatrix} = 
\begin{bmatrix} 
\widetilde{U}_{L}\\
D^{-1}_{n-L}E^{\top}\widetilde{U}_{L}\Lambda_{L}^{-1}
\end{bmatrix} = 
\begin{bmatrix} 
D^{-1}_{L} &  \\
 & D^{-1}_{n-L}
\end{bmatrix}
\begin{bmatrix} 
W_{L}  \\
E^{\top}
\end{bmatrix}
\widetilde{U}_{L}\Lambda^{-1}_{L}\,.$$ In other words, we run the DM on a small subset of the whole database, and extend it to the whole database.  While it is slightly different from the original Nyström extension, we still call it the Nyström extension.  Note that in practice, we only need to calculate $W_{L}$ and $E^{\top}$ instead of $W$ and $D$, which is more efficient in the sense of both time and spatial complexities.  With the estimated eigenvectors on the whole dataset, we can define the associated embedding and hence the distance just like the above.  Specifically, suppose we have $\tilde{U}=\begin{bmatrix} \tilde{u}_1 & \ldots& \tilde{u}_L \end{bmatrix}\in \mathbb{R}^{N\times L}$ and $\tilde{\Lambda}=\texttt{diag}\begin{bmatrix}\lambda_1 &\ldots &\lambda_L &0&\ldots &0\end{bmatrix}$. Then we can define the associated DM and DD by $$\Phi^{(\textup{Nystr\"om})}_t:x_{i}\mapsto e_i^\top \tilde{U}_{m}\tilde{\Lambda}_{m}^t\,,$$ where $t>0$ is the chosen diffusion time, $\tilde{U}_m\in \mathbb{R}^{n\times m}$ to be a matrix consisting of $\tilde{u}_2,\ldots,\tilde{u}_{m+1}$ and $\tilde{\Lambda}_m:=\texttt{diag}(\lambda_2,\ldots,\lambda_{m+1})$, and the associated DD $$D^{(\textup{Nystr\"om})}_t(x_i,x_j):=\|\Phi^{(\textup{Nystr\"om})}_t(x_{i})-\Phi^{(\textup{Nystr\"om})}_t(x_j)\|_{\mathbb{R}^{m}}\,.$$

Complexity analysis
-------------------

Suppose we have $n$ data points and the landmark set is of size $n^\beta$, where $\beta\leq 1$, both for the Nyström extension and the Roseland.  For the kNN scheme, we set $k=n^\beta$ for a fair comparison.  If we do not carry out the kNN scheme, or if the kernel is not compactly supported, the spatial complexity of the DM is $O(n^2)$. If we construct the affinity matrix for the DM by the kNN scheme, the spatial complexity of the DM becomes $O(n^{1+\beta})$. On the other hand, no matter what kernel we choose, compactly supported or not, the spatial complexity of the Roseland and the Nyström extension is $O(n^{1+\beta})$.

For the computational complexity, it is generally divided into step parts.  The first part is calculating the associated affinity matrix and degree matrix, and the second part is calculating the eigendecomposition or SVD.  In the ordinary DM, the construction of the affinity matrix and the degree matrix is $O(n^2)$. If the kNN construction is considered and the k-d tree based algorithm is applied, the averaged time complexity of constructing the affinity matrix and the degree matrix is $O(n\log(n)+n^{1+\beta})=O(n^{1+\beta})$. In the Nyström extension, the construction of the $W_L$ and hence its degree matrix is $O(n^{2\beta})$ when $L=n^\beta$ for $\beta\leq 1$, while the construction of $E$ and $D_{n-L}$ is $O(n^{1+\beta})$. Thus, the first part complexity for the Nyström extension is $O(n^{2\beta}+n^{1+\beta})$. In the Roseland, the construction of the landmark-set affinity matrix and its associated degree matrix is $O(n^{1+\beta})$ when the size of the landmark set is $n^\beta$. For the second part, it falls in the discussion of the complexity of the general eigendecomposition and SVD.  Given a symmetric kernel matrix $M\in\mathbb{R}^{N\times N}$, the eigen-decomposition complexity is usually $O(N^3)$,[^1] and when $M$ is $k$ sparse, where $k\leq n$, the complexity can be improved to $O(N^{2+\eta'})$ for an arbitrary $\eta'>0$ when $k\leq N^{0.14}$ [@Yuster2004].  In our application, even if we make $k\leq N^{0.14}$, the eigen-decomposition of the $M$ is roughly ${O}(N^{2+\eta'})$. On the other hand, for a matrix of size $N\times N'$, where $N\geq N'$, then the complexity of the SVD for is ${O}(NN'^2)$. By combining these two parts, the overall computational complexity is $O(n^3)$ for the ordinary DM, $O(n^\omega)$ for the DM with the kNN scheme, where $\omega>2$ depends on the chosen $\beta$, $O(n^{1+\beta}+n^{3\beta})$ for the Nyström extension, and $O(n^{1+2\beta})$ for the Roseland.  The complexity of the HKC is the same as that of the Roseland.  In sum, while the Roseland is not faster than the traditional Nyström extension approach, it is comparable, particularly when $\beta$ is small.  On the other hand, both the Nyström extension and the Roseland are more efficient than the ordinary DM and the DM with the kNN scheme. 

Denoising by the Roseland
-------------------------

As has been summarized in Section \[section: related work: robustness\], spectral embedding algorithms are usually challenged by inevitable noise.  If we do not handle it properly, a spectral embedding might result in a bad, or even misleading result [@el2016graph].  While there have been quite a few efforts toward that end; however, there is definitely a space for an improvement.  The proposed Roseland automatically enjoys this robustness property.  Specifically, since each data point first diffuses to [*all*]{} landmarks and then diffuses back, the danger of finding [*wrong*]{} neighbors inherited in the kNN scheme does not exist.  In other words, this step can be viewed as a surrogate of finding true neighbors.  Second, since by construction the graph is complete within 2 steps, the argument in [@el2016graph] can be directly generalized to explain the robustness of the Roseland.  Moreover, under the condition that we are able to collect a clean but small dataset in addition to the large but noisy dataset, the Roseland can help us handle noise by taking the clean but small dataset as the landmark set.  We comment that this is not an artificial condition; instead, it is the situation we commonly encounter in real life.  For example, in the medical field, collecting a clean dataset of high quality is usually labor-intense and expensive.  However, it is relatively easy to collect a large dataset from a rather cheap equipment, in exchange of the data quality. 

Reference set as subset of the data
-----------------------------------

In general, we may not be able to acquire the landmark set in addition to the dataset, and have to select the landmark set from the available dataset; that is the landmark set $\mathcal{Y}\subseteq\mathcal{X}$ and is not independent of the dataset $\mathcal{X}$. In this case, we propose to first sample $m$ landmarks, denoted by $\mathcal{Y}$, from $\mathcal{X}$ so that $\mathcal{Y}$ is independent of $\mathcal{X} \setminus \mathcal{Y}$. Then, we apply the Roseland on $\mathcal{X} \setminus \mathcal{Y}$ using $\mathcal{Y}$, and extend the embedding to $\mathcal{Y}$ by the Nyström extension.  When $|\mathcal{Y}| \ll  |\mathcal{X}|$, the discrepancy of this approach and the independence setup is negligible, and will asymptotically vanish. 

Mathematical setups and notations {#setup}
=================================

Manifold model {#dm_setup}
--------------

Denote our observed data set the point could $\mathcal{X} = \{x_{i}\}_{i=1}^{n}\subseteq\mathbb{R}^{q}$, which are independent and identically distributed (i.i.d.)  sampled from a random vector $X:(\Omega,\mathcal{F},\mathbb{P})\rightarrow\mathbb{R}^{q}$. We assume that the range $X$ is supported on a $d$-dimensional compact smooth Riemannian manifold $(M^{d}, g)$ without boundary that is isometrically embedded in $\mathbb{R}^p$ via $\iota:M^{d}\hookrightarrow \mathbb{R}^p$. Hence, $X$ induces a probability measure on $\iota(M^{d})$, denoted by $\widetilde{\mathbb{P}}_{X}$. Further assume $\widetilde{\mathbb{P}}_{X}$ is absolutely continuous with respect to the Riemannian measure on $\iota(M)$, denoted by $\iota_{\ast}dV(x)$. Then, by the Radon Nikodym theorem, we have $d\widetilde{\mathbb{P}}_{X}(x)=p_{X}(\iota^{-1}(x))\iota_{\ast}dV(x)$. Clearly, $p_X$ is a function defined on $M^{d}$.

We call $p_{X}$ defined above the probability density function (p.d.f.)  associated with $X$. When $p_{X}$ is constant, $X$ is called uniform; otherwise non-uniform. 

We assume $p_{X}$ satisfies $p_{X} \in \mathcal{C}^4(M^{d})$ and $0 < \inf_{x\in M^{d}}p_{X}(x)\leq \sup_{x\in M^{d}}p_{X}(x)$.

\[Definition:kernel\] A kernel function is any non-negative function $K:[0,\infty)\rightarrow\mathbb{R}^{+}$ so that it is $C^{3}$, $K(0)>0$ and decays exponentially fast.  Denote $\mu_{r,l}^{(k)}\coloneqq \int_{\mathbb{R}^{d}}\|x\|^{l}\partial_{k}K^{r}(\|x\|)\,dx$, for $r,l,k=0,1,2,\ldots$, and assume $K$ is normalized so that $\mu_{1,0}^{0}=1.$

Roseland in the manifold setup
------------------------------

For the landmark set $\mathcal{Y} = \{y_{j}\}_{j=1}^{m}\subseteq\mathbb{R}^{q}$, we assume that $y_{j}$’s are i.i.d.  samples from a random vector $Y:(\Omega,\mathcal{F},\mathbb{P})\rightarrow\mathbb{R}^{q}$, whose range is supported on the same manifold $M^{d}$, and has p.d.f.  $p_{Y}$ on $M^{d}$. Moreover, we assume that $Y$ is independent of $X$. We also assume $p_{Y}$ satisfies $p_{Y} \in \mathcal{C}^4(M^{d})$. Clearly, we have $0 < \inf_{x\in M^{d}}p_{Y}(x)\leq \sup_{x\in M^{d}}p_{Y}(x)$.

In the Roseland, the “affinity” between two data points $x_{i}$ and $x_{j}$ is measured via the landmark set $\mathcal{Y} = \{y_{j}\}_{j=1}^{m}$. Specifically, note that the affinity matrix $W^{(\textup{R})}=W^{(\textup{r})}[W^{(\textup{r})}]^{\top}\in\mathbb{R}^{n\times n}$, where $W^{(\textup{r})}\in\mathbb{R}^{n\times m}$ such that $W^{(\textup{r})}(i,k)= K_{\epsilon}(x_{i},y_{k})$ and $[W^{(\textup{r})}]^{\top}$ is the transpose of $W^{(\textup{r})}$.

\[ref\_affinity\] Take the kernel function $K$. The affinity between any two points $x_{i}$ and $x_{j}$ via a landmark set $\mathcal{Y} = \{y_{j}\}_{j=1}^{m}$ is defined by $$K_{\textup{ref},\epsilon,n}(x_{i},x_{j})\coloneqq \frac{1}{m}\sum_{k=1}^{m}K_{\epsilon}(x_{i},y_{k})K_{\epsilon}(y_{k},x_{j}).$$

See Figure \[Figure:refDMillustration\] for an illustration of how this affinity is determined.  To study the asymptotical behavior of the Roseland, we take the following expansion into account.  For $f\in C(M)$, denote its discretization on $\{x_i\}_{i=1}^n$ as $\textbf{\textit{f}}\in\mathbb{R}^{n}$ such that $\textbf{\textit{f}}_{i}=f(x_{i})$. By a direct expansion, we have $$\begin{aligned}
\label{discrete_op2}
    \big[(D^{(\textup{R})})^{-1}W^{(\textup{R})}\textbf{\textit{f}}\;\big](i)=\frac{\sum_{j=1}^{n}W^{(\textup{R})}_{ij}\textbf{\textit{f}}_{j}}{\sum_{j=1}^{n}W^{(\textup{R})}_{ij}}\,.$$ The denominator and numerator ring the bell of the law of large numbers.  We thus expect that (\[discrete\_op2\]) converges to the following integral operators when $n\to\infty$:

\[GL\_cts\] Take $f\in C(M), \epsilon>0$, define $$T_{\textup{ref},\epsilon}f(x)\coloneqq\int_{M}\frac{K_{\textup{ref},\epsilon}(x,y)}{d_{\textup{ref},\epsilon}(x)}f(y)p_{X}(y)\,dV(y),$$ where $K_{\textup{ref},\epsilon}:M\times M\to \mathbb{R}$ defined as $$\begin{aligned}
K_{\textup{ref},\epsilon}(x,y)&\coloneqq\int_{M}K_{\epsilon}(x,z)K_{\epsilon}(z,y)p_{Y}(z)\,dV(z)\end{aligned}$$ is called the [*landmark-kernel*]{} induced by the landmark set, and $$\begin{aligned}
d_{\textup{ref},\epsilon}(x)&\coloneqq\int_{M}K_{\textup{ref},\epsilon}(x,y)p_{X}(y)\,dV(y)\,.\end{aligned}$$

Kernel behavior with the landmark set
-------------------------------------

It is worth a bit more discussion of the landmark-kernel induced by the landmark set.  Recall Definition \[ref\_affinity\].  The affinity between two points is now determined by passing through the landmark set.  A direct consequence is that the kernel function may vary from point to point, depending on how two points are geologically related to the landmark set.  The affinity of a point to itself might by smaller than the the affinity between a point and its close neighbor.  To further illustrate this finding, take the Gaussian as the kernel function; that is, $K_{\epsilon}(x,y)=e^{-\|x-y\|_{\mathbb{R}^{q}}^{2}/\epsilon}$. Then, the landmark-kernel function at each point is a mixture of Gaussian.  Next, take a set of equally spaced samples from $S^{1}$, denoted as $\mathcal{X}$, and order them by their angles.  Let the landmark set $\mathcal{Y}\subseteq\mathcal{X}$ contains $5\%$ equally spaced points of $\mathcal{X}$, also ordered by angle to its center.  See Figure \[S1ker\] $(a)$. By the symmetry of $\mathcal{X}$ and $\mathcal{Y}$, $|\mathcal{Y}|=5\%\times|\mathcal{X}|$ means we have $20$ difference kernel functions.  Indeed, we have 20 data points between 2 consecutive landmarks, so there are totally 20 different geological relationships between the dataset and the landmark set.  See Figure \[S1ker\] $(b)$ for plots of the $20$ landmark-kernels at the data points (in order) between two consecutive landmarks.  Similarly, we plot the kernel functions when $|\mathcal{Y}|=a\times|\mathcal{X}|$, for $a=10\%,20\%$, see Figure \[S1ker\] $(c)$ and $(d)$. Note that when $N$ is fixed and $M$ increases, it is not surprising that the kernel looks more like a Gaussian. 

![ Top left: $S^{1}$ with 5% equally spaced landmarks.  Top right: $20$ different landmark-kernel functions from 5% equal spaced landmarks as in top left.  Bottom left: $10$ different landmark-kernels if choose 10% equally spaced landmarks.  Bottom right: $5$ different landmark-kernels if choose 20% equal spaced landmarks.[]{data-label="S1ker"}](data_ref.pdf "fig:"){width="4.5cm"} ![ Top left: $S^{1}$ with 5% equally spaced landmarks.  Top right: $20$ different landmark-kernel functions from 5% equal spaced landmarks as in top left.  Bottom left: $10$ different landmark-kernels if choose 10% equally spaced landmarks.  Bottom right: $5$ different landmark-kernels if choose 20% equal spaced landmarks.[]{data-label="S1ker"}](ker1.pdf "fig:"){width="5.5cm"}\
![ Top left: $S^{1}$ with 5% equally spaced landmarks.  Top right: $20$ different landmark-kernel functions from 5% equal spaced landmarks as in top left.  Bottom left: $10$ different landmark-kernels if choose 10% equally spaced landmarks.  Bottom right: $5$ different landmark-kernels if choose 20% equal spaced landmarks.[]{data-label="S1ker"}](ker2.pdf "fig:"){width="5cm"} ![ Top left: $S^{1}$ with 5% equally spaced landmarks.  Top right: $20$ different landmark-kernel functions from 5% equal spaced landmarks as in top left.  Bottom left: $10$ different landmark-kernels if choose 10% equally spaced landmarks.  Bottom right: $5$ different landmark-kernels if choose 20% equal spaced landmarks.[]{data-label="S1ker"}](ker3.pdf "fig:"){width="5cm"}

Asymptotical behavior of the Roseland {#Section:PointConvergence}
=====================================

Before stating our main asymptotical results of the Roseland, we briefly summarize existing literature about the DM.  A celebrated spectral embedding [@berard1994embedding] gives a solid foundation of various spectral based unsupervised learning algorithms, particular the DM.  It says that one can embed any given smooth closed $n$-dimensional Riemannian manifold by the eigenfunctions of its Laplace–Beltrami operator, and the embedding can be “tuned” to be as isometric as possible.  However, this spectral embedding needs [*all*]{} eigenfunctions, which is not numerically affordable.  To resolve this issue, it is proved in [@portegies2016embeddings] that for a given tolerable metric recovery error, we can achieve an almost isometric embedding with that tolerable error with a finite number of eigenfunctions of the Laplace–Beltrami operator, where the number only depends on the geometric bounds and the dimension.  The above two results are on the continuous setup.  To utilize these results, we need to link numerical finite sampling dataset to the continuous manifold setup.  Specifically, if we are able to prove how the eigenvectors of the GL converges in the spectral sense to the eigenfunctions of the Laplace-Beltrami operator, we can apply the above-mentioned results in the continuous setup.  In [@singer2016spectral; @trillos2018error], the authors provide the spectral convergence of the GL constructed from random samples to the Laplace–Beltrami operator, and a convergence rate is also discussed in [@trillos2018error].  To handle the inevitable noise in practice, the robustness of the DM, and more generally the GL, is shown in [@el2016graph].  In this section we show that the eigenvalues and eigenvectors of the Roseland converge to the eigenvalues and eigenfunctions of the Laplace-Beltrami operator, and quantify the convergence rate, both pointwisely and spectrally.  At the first glance, it might be expected that the proofs are similar to those shown in [@singer2016spectral; @trillos2018error].  However, as we will see below, we run into the dependence issue when the landmark set is considered, and extra efforts and new technical tools are needed. 

Pointwise convergence
---------------------

We first state the bias analysis of the Roseland. 

\[biasthm\](Bias analysis) Take $f\in C^{4}(M^{d})$. Then, for all $x\in M^{d}$ we have $$\begin{aligned}
    T_{\textup{ref},\epsilon}f(x) - f(x) =&\, \frac{\epsilon\mu_{1,2}^{(0)}}{d}\left(\frac{2\nabla p_{X}(x)}{p_{X}(x)}+\frac{\nabla p_{Y}(x)}{p_{Y}(x)}\right)\cdot\nabla f(x) \\
    &+ \frac{\epsilon\mu_{1,2}^{(0)}}{d}\Delta f(x) + \mathcal{O}(\epsilon^{2})\,.\nonumber\end{aligned}$$

The proof is postponed to Section \[biasproof\]. 

\[remark: dm\_thm\] We compare the obtained result with the existing theorems for the DM shown in [@Coifman2006].  Take $f\in C^{4}(M^{d})$. Recall the definition: $$T_{\epsilon,\alpha}f(x)\coloneqq\int_{M}\frac{K_{\epsilon,\alpha}(x,y)}{d_{\epsilon,\alpha}(x)}f(y)p_{X}(y)\,dV(y)\,,$$ where $K_{\epsilon,\alpha}(x,y)\coloneqq\frac{K_{\epsilon}(x,y)}{p_{X,\epsilon}^{\alpha}(x)p_{X,\epsilon}^{\alpha}(y)}$, $p_{X,\epsilon}(x)\coloneqq\int_{M}K_{\epsilon}(x,y)p_{X}(y)\,dV(y)$, and $d_{\epsilon,\alpha}(x)\coloneqq\int_{M}K_{\epsilon,\alpha}(x,y)p_{X}(y)\,dV(y)$. Then for all $x\in M^{d}$ we have: $$T_{\epsilon,\alpha}f(x) - f(x) = \frac{\epsilon\mu_{1,2}^{(0)}}{2d}\left(\Delta f(x) + \frac{2\nabla f(x)\cdot\nabla p_{X}^{1-\alpha}(x)}{p_{X}^{1-\alpha}(x)}\right) + \mathcal{O}(\epsilon^{2})\,.$$ Based on this result, we see that in the DM, we can remove the impact of the non-uniformly sampling of the data set by letting $\alpha=1$ in the $\alpha$-normalization step.  In the Roseland, if $\frac{2\nabla p_{X}(x)}{p_{X}(x)}+\frac{\nabla p_{Y}(x)}{p_{Y}(x)}=0$, then we remove the impact of the non-uniformly sampling and recover the Laplace-Beltrami operator.  Note that $\frac{2\nabla p_{X}(x)}{p_{X}(x)}+\frac{\nabla p_{Y}(x)}{p_{Y}(x)}=0$ suggests that we may want to have the landmark set following $p_{Y}(x)\propto\frac{1}{p^{2}_{X}(x)}$. This serves us as the guidance of how to design the landmark set. 

\[variancethm\](Variance analysis) Take $\mathcal{X}=\{x_i\}_{i=1}^n$ and $\mathcal{Y}=\{y_j\}_{i=1}^m$, where $m=[n^{\beta}]$ for some $0<\beta\leq 1$ and $[x]$ is the nearest integer of $x\in\mathbb{R}$. Take $f\in C(M^{d})$ and denote $\textbf{\textit{f}}\in\mathbb{R}^{n}$ such that $\textbf{\textit{f}}_{i}=f(x_{i})$. Let $\epsilon = \epsilon(n)$ so that $\frac{\sqrt{\log n}}{n^{\beta/2}\epsilon^{d/2+1/2}}\rightarrow 0$ and $\epsilon\rightarrow 0$ when $n\rightarrow\infty$. Then with probability higher than $1 - \mathcal{O}(1/n^{2})$, we have $$\big[(I-(D^{(\textup{R})})^{-1}W^{(\textup{R})})\textbf{\textit{f}}\;\big](i) = f(x_{i}) - T_{\textup{ref},\epsilon}f(x_{i}) + \mathcal{O}\Big(\frac{\sqrt{\log n}}{n^{\beta/2}\epsilon^{d/2-1/2}}\Big)$$ for all $i=1,2,\ldots,n$.

The proof is postponed to Section \[biasproof\]. 

\[Remark about the variance convergence rate\] We compare the obtained result with the existing theorems shown in [@Coifman2006; @singer2016spectral].  For the variance analysis for the DM, we have the following result when there are $n$ data points $\mathcal{X}=\{x_i\}_{i=1}^n$ from the manifold.  Take $f\in C(M^{d})$.

-   For $0<\alpha\leq 1$, let $\epsilon = \epsilon(n)$ so that $\frac{\sqrt{\log n}}{n^{1/2}\epsilon^{d/4+1/2}}\rightarrow 0$ and $\epsilon\rightarrow 0$ when $n\rightarrow\infty$. Then with probability higher than $1 - \mathcal{O}(1/n^{2})$, for all $i=1,2,\ldots,n$, we have $$\big[(I-(D^{(\alpha)})^{-1}W^{(\alpha)})\textbf{\textit{f}}\;\big](i) = T_{\epsilon,\alpha}f(x_{i}) - f(x_{i}) + \mathcal{O}\Big(\frac{\sqrt{\log n}}{n^{1/2}\epsilon^{d/4}}\Big)\,.$$

-   For $\alpha=0$, let $\epsilon = \epsilon(n)$ so that $\frac{\sqrt{\log n}}{n^{1/2}\epsilon^{d/4+1/2}}\rightarrow 0$ and $\epsilon\rightarrow 0$, when $n\rightarrow\infty$. Then with probability higher than $1 - \mathcal{O}(1/n^{2})$, for all $i=1,2,\ldots,n$, we have $$\big[(I-(D^{(0)})^{-1}W^{(0)})\textbf{\textit{f}}\;\big](i) = T_{\epsilon,\alpha}f(x_{i}) - f(x_{i}) + \mathcal{O}\Big(\frac{\sqrt{\log n}}{n^{1/2}\epsilon^{d/4-1/2}}\Big)\,.$$

Clearly, unlike the DM, in the Roseland, its convergence rate depends on $n^{\beta}$, which is the size of the landmark set.  This pointwise convergence result tells us that the smaller the landmark set is, the faster the algorithm, but the slower the convergence rate to the Laplace-Beltrami operator.  We should compare the rate of Roseland with the rate of DM when the alpha normalization is $0$. The error term in the DM is of order $ \mathcal{O}\Big(\frac{\sqrt{\log n}}{n^{1/2}\epsilon^{d/4-1/2}}\Big)$ while the error term in the Roseland is of order $\mathcal{O}\Big(\frac{\sqrt{\log n}}{n^{\beta/2}\epsilon^{d/2-1/2}}\Big)$, where $n$ is the size of data, $n^{\beta}$ is the size of landmark set.  Note that even when we let $\beta=1$, the convergence rate of Roseland still does not recover the convergence rate of DM, where they differ by a factor of $\epsilon^{\frac{d}{4}}$. This is because the Roseland introduces dependence relation among data points by diffusing through the landmark set, and this dependence relation results in a larger variance of the random variable to be analyzed.  This fact slows down the convergence rate when we apply the large deviation bound.  See section \[Discussion optimal variance convergence discrepancy\] for more details. 

Idea of analyzing the variance {#grid_samp_def}
------------------------------

Let $X$ and $Y$ be two independent random variables and $f:(X,Y)\rightarrow\mathbb{R}$. One way to compute $\mathbb{E}(f(X,Y))$ numerically is by i.i.d.  sampling $n$ pairs of points $\{(x_{i},y_{i})\}_{i=1}^{n}$ from the joint distribution of $(X,Y)$. Then we have $\frac{1}{n}\sum_{i=1}^{n}f(x_{i},y_{i})\rightarrow\mathbb{E}(f(X,Y))$ almost surely by the law of large numbers.  And there are standard techniques available to compute its convergence rate.  Due to the nature of landmark set, this approach does not hold.  Indeed, note that if we expand , we have $$\begin{aligned}
\label{Expansion:refDM in the summation formate}
    \big[D^{(\textup{R})})^{-1}W^{(\textup{R})}\textbf{\textit{f}}\big](i)=\frac{\frac{1}{n}\sum_{j=1}^{n}\big[\sum_{k=1}^{m}K_{\epsilon}(x_{i},y_{k})K_{\epsilon}(y_{k},x_{j})\big]\textbf{\textit{f}}_{j}}{\frac{1}{n}\sum_{j=1}^{n}\big[\sum_{k=1}^{m}K_{\epsilon}(x_{i},y_{k})K_{\epsilon}(y_{k},x_{j})\big]}\,,\end{aligned}$$ which generates dependence among the summands.  Specifically, we have the following definition:

\[gridsamp\] Let $X$ and $Y$ be two independent random variables.  We call $\{(x_i,y_j)|\;i=1,\ldots,n,\, j=1,\ldots,m\}$ a grid sampling if $\{x_{i}\}_{i=1}^{n}$ is i.i.d.  sampled from $X$ and $\{y_{j}\}_{j=1}^{m}$ is i.i.d.  sampled from $Y$.

Clearly, the grid samples are not independent, and we know that $(x_{i_{1}},y_{j_{1}})$ is independent of $(x_{i_{2}},y_{j_{2}})$ if and only if $i_{1}\neq i_{2}$ and $j_{1}\neq j_{2}$. In general, can be formulated in the following way.  Given $f:(X,Y)\rightarrow\mathbb{R}$ and a grid sampling $\{(x_i,y_j)|\;i=1,\ldots,n,\, j=1,\ldots,m\}$, we ask how well we can approximate $\mathbb{E}(f(X,Y))$ from the sampling grid; that is, what is the convergence rate of $$\begin{aligned}
\label{gridconvg}
    \frac{1}{mn}\sum_{i=1}^{n}\sum_{j=1}^{m}f(x_{i},y_{j})\longrightarrow\mathbb{E}(f(X,Y))\,.\end{aligned}$$ Clearly, we need to handle the dependence on the grid sampling. 

To answer this question, we consider the work in [@janson2004large], which provides a method of computing the convergence rate of this kind of sampling.  In general, we consider the random variable of the form $$\begin{aligned}
\label{dep_sum}
    X=\sum_{\alpha\in\mathcal{A}}Y_{\alpha}\end{aligned}$$ where $Y_{\alpha}$ are random variables, independent of not, and $\alpha$ ranging over some index set $\mathcal{A}$. We have the following definition. 

Given an index set $\mathcal{A}$ and $\{Y_{\alpha}\}_{\alpha\in\mathcal{A}}$, we make the following definitions. 

-   A subset $\mathcal{A}'$ of $\mathcal{A}$ is independent if the corresponding random variables $\{Y_{\alpha}\}_{\alpha\in\mathcal{A}'}$ are independent. 

-   A family $\{\mathcal{A}_{j}\}_{j}$ of subsets of $\mathcal{A}$ is a cover of $\mathcal{A}$ if $\bigcup_{j}\mathcal{A}_{j}=\mathcal{A}$.

-   A cover is proper if each set $\mathcal{A}_{j}$ in it is independent. 

-   $\chi(\mathcal{A})$ is the size of the smallest proper cover of $\mathcal{A}$.

Then we have the first Hoeffding-like concentration inequality. 

\[hoeffding\] Suppose X is defined in with $a_{\alpha}\leq Y_{\alpha}\leq b_{\alpha}$ for every $\alpha\in\mathcal{A}$, where $a_{\alpha},b_{\alpha}\in \mathbb{R}$. Then for all $t>0$, $$\begin{aligned}
    \nonumber\mathbb{P}(X-\mathbb{E}(X)\geq t)\leq\textup{exp}\left(\frac{-2t^{2}}{\chi(\mathcal{A})\sum_{\alpha\in\mathcal{A}}(b_{\alpha}-a_{\alpha})^{2}}\right).\end{aligned}$$ The same estimate holds for $\mathbb{P}(X-\mathbb{E}(X)\leq -t)$.

When $Y_{\alpha}$’s have variances that are substantially smaller than $(b_{\alpha}-a_{\alpha})^{2}/4$, we can improve theorem \[hoeffding\] to the Bernstein-like concentration inequality. 

\[bernstein\] Suppose X is defined in with $Y_{\alpha}-\mathbb{E}(Y_{\alpha})\leq b$ for some $b>0$ for all $\alpha\in\mathcal{A}$. Suppose $S\coloneqq\sum_{\alpha\in\mathcal{A}}\textup{Var}\,Y_{\alpha}<\infty$. Then, for all $t>0$, $$\begin{aligned}
    \nonumber\mathbb{P}(X-\mathbb{E}(X)\geq t)\leq\textup{exp}\left(\frac{-8t^{2}}{25\chi(\mathcal{A})(S+bt/3)}\right).\end{aligned}$$ The same estimates holds for $\mathbb{P}(X-\mathbb{E}(X)\leq -t)$ if also $Y_{\alpha}-\mathbb{E}(Y_{\alpha})\geq -b$ for $b>0$.

With this general theory, we now come back to our setup.  In our grid sampling scheme in the Roseland, $\mathcal{A} = \{(j,k)\}_{j=1,k=1}^{n,m}$. We now claim that $$\label{Optimal bound of the independent set A}
\chi(\mathcal{A})\leq m+n-1=\mathcal{O}(\max(m,\,n))\,.$$ The easiest way of seeing it is by the following grid, where the coordinate $(j_{a}, k_{b})$ corresponds to the random variable $(x_{j_{a}}, y_{k_{b}})$. Clearly, we know that $\{\mathcal{A}_{x}\}_{x=1}^{n+m-1}$ is a proper cover of $\mathcal{A}$. That means the convergence rate in (\[gridconvg\]) should be the same as that of $\frac{1}{\min(n,m)}\sum_{i=1}^{\min(n,m)}f(X_{i},Y_{i})\rightarrow\mathbb{E}(f(X,Y))$, and hence the rate is dominated by $\min(m,n)$. Note that the bound $m+n$ is not optimal, but $\chi(\mathcal{A})=\mathcal{O}(\max(m,n))$ is, which can be seen by the special case where $m=1$.

![ Illustration of the grid argument, note the grid samples within each circle are independent. ]( cover.jpg){width="11cm"}

The challenge we encounter with the grid sampling is directly related to the U-statistic or V-statistic.  Take the “kernel” $h$ of $r$ variables, where $r\in \mathbb{N}$. For the dataset $x_{1},\ldots,x_{n}$, where $n\geq r$, an U-statistic of order $r\in \mathbb{N}$ is defined as $$\begin{aligned}
    U_{r}: = \frac{1}{{n\choose r}}\sum_{(i_1, i_2,\ldots, i_r)\in \langle n\rangle}h(x_{i_{1}},\ldots,x_{i_{r}})\,,\label{definition U statistics}\end{aligned}$$ where $\langle n\rangle$ is the set of all permutations of $\{1,\ldots,n\}$ and $h$ is symmetric in its arguments.  For example, when the kernel $h$ is of $1$ variable and $h(x)=x$, then the U-statistic $U_{1}(x)=(x_{1}+\ldots+x_{n})/n$ is the sample mean $\overline{x}$. In practice, this statistic has a representation as the V-statistic: $$\begin{aligned}
    V_{m,n} := \frac{1}{n^m}\sum_{i_{1}=1}^{n}\cdot\cdot\cdot\sum_{i_{m}=1}^{n}h(x_{i_{1}},\ldots,x_{i_{m}})\,,\end{aligned}$$ where $h$ is a symmetric kernel function.  We call $V_{mn}$ a V-statistic of degree $m$. A typical example of a degree-2 V-statistic is the second central moment; that is, take $h(x,y):=(x-y)^{2}/2$, then $V_{2,n} = \frac{1}{n^{2}}\sum_{i=1}^{n}\sum_{j=1}^{n}\frac{1}{2}(x_{i}-x_{j})^{2} = \frac{1}{n}\sum_{i=1}^{n}(x_{i}-\overline{x})^{2}$. We should notice the difference between U-statistic, V-statistic and grid sampling.  In the U-statistic, we need to take average over distinct ordered samples of size $r$ taken from $\{1,\ldots,n\}$; in the V-statistic, each argument of the function $h$ will run over all sample points.  Clearly, both are not the case in the grid sampling we run into in the Roseland.  Moreover, in the grid sampling, $h$ need not to be symmetric.  However, the U-statistic and the V-statistic are both special cases of the form in equation (\[dep\_sum\]). 

Spectral Convergence {#Section:spectral convergence}
--------------------

The point-wise convergence of the Roseland to the Laplace-Beltrami operator in Section \[Section:PointConvergence\] does not guarantee the spectral convergence.  To fully understand the spectral based methods, we need to establish the spectral convergence, and this is the focus of this section.  Let $\{v_{n}\}_{n\in\mathbb{N}}$ be a set of eigenvectors of the transition matrix $(D^{(\textup{R})})^{-1}W^{(\textup{R})}$ associated with the point cloud $\mathcal{X}=\{x_{i}\}_{i=1}^{n}\subseteq\mathbb{R}^{q}$. We would like to study when $n\to \infty$, how will the eigenvectors $\{v_{n}\}_{n\in\mathbb{N}}$ converge to the eigenfunctions of the Laplace-Beltrami operator.  Note that the vectors $v_{n}$ are in difference Euclidean spaces for different $n$, and the eigenfunctions of the Laplace-Beltrami operator are smooth functions on $M$. Clearly, they cannot be compared directly and we need to manipulate those quantities a bit so that we can compare them.  In brief, we will find a sequence of functions $f_{n}\in C(M)$, such that the restriction of $f_{n}$ on the point cloud $\mathcal{X}$ equals to entries of $v_{n}$; that is, $f_{n}(x_{i}) = v_{n}(i)$, for $i=1,\ldots,n$. Then, we study the convergence of $\{f_{n}\}$ as $n\rightarrow\infty$. To state our spectral convergence theorem, we need the following definitions and results. 

\[Definition various operators discrete\] Take $\mathcal{X}=\{x_{i}\}_{i=1}^{n}\subset M$. Define the following functions $$\begin{aligned}
\widehat{K}_{\textup{ref},\epsilon,n}(x,y)&\,:=\frac{1}{m}\sum_{j}^{m}K_{\epsilon}(x,z_{j})K_{\epsilon}(z_{j},y)\in C(M\times M),\nonumber\\
\widehat{d}_{\textup{ref},\epsilon,n}(x)&\,:=\frac{1}{n}\sum_{i=1}^{n}\widehat{K}_{\textup{ref},\epsilon,n}(x,x_{i})\in C(M),\\
\widehat{M}_{\textup{ref},\epsilon,n}(x,y)&\,:=\frac{\widehat{K}_{\textup{ref},\epsilon,n}(x,y)}{\widehat{d}_{\textup{ref},\epsilon,n}(x)}\in C(M\times M)\,.\nonumber\end{aligned}$$ Also define the following operator: $$\begin{aligned}
\widehat{T}_{\textup{ref},\epsilon,n}f(x)&\,:=\frac{1}{n}\sum_{i=1}^{n}\widehat{M}_{\textup{ref},\epsilon,n}(x,x_{i})f(x_{i})\,.\end{aligned}$$ Moreover, define the restriction operator $\rho_{n}:C(M)\rightarrow\mathbb{R}^{n}$ $$\rho_{n}:f\mapsto [f(x_{1}), f(x_{2}),\ldots, f(x_{n})]^{\top}\,.$$

The following Lemma describes the relationship between $(D^{(\textup{R})})^{-1}W^{(\textup{R})}$ and the integral operator $\widehat{T}_{\textup{ref},\epsilon,n}$.

\[one to one corresp\] Let $U_{n}:=(D^{(\textup{R})})^{-1}W^{(\textup{R})}$, then $U_{n}\circ\rho_{n}=\rho_{n}\circ\widehat{T}_{\textup{ref},\epsilon,n}$. Moreover, we have the following one to one correspondence. 

1.   If $f\in C(M)$ is an eigenfunction of $\widehat{T}_{\textup{ref},\epsilon,n}$ with the eigenvalue $\lambda$, then the vector $v\coloneqq\rho_{n}f$ is an eigenvector of $U_{n}$ with the eigenvalue $\lambda$. Moreover, suppose $\lambda\neq 0$ is an eigenvalue of $\widehat{T}_{\textup{ref},\epsilon,n}$ with the eigenfunction $f$. If we let $v\coloneqq\rho_{n}f$, then $f$ satisfies $$\label{extrapolation of eigenvector to the whole manifold}
        f(x)=\frac{\sum_{j=1}^{n}\widehat{K}_{\textup{ref},\epsilon,n}(x,x_{j})v_{j}}{\lambda\sum_{j=1}^{n}\widehat{K}_{\textup{ref},\epsilon,n}(x,x_{j})}\,.$$

2.   If $v$ is an eigenvector of $U_{n}$ with the eigenvalue $\lambda\neq 0$, then $f$ defined in is an eigenfunction of $\widehat{T}_{\textup{ref},\epsilon,n}$ with the eigenvalue $\lambda$.

The proof of this proposition is standard, and can be found in, for example [@von2008consistency; @singer2016spectral], so we omit it.  With this Lemma, we now can simply study how the eigenstructure of $T_{\textup{ref},\epsilon,n}$ converges to that of $-\Delta$.

Denote $(\lambda_{i}, u_{i})$ to be the $i$-th eigenpair of $-\Delta$, where $\lambda_i$ is the $i$-th smallest eigenvalue.  Note that under our manifold setup, by the well known elliptic theory, the spectrum of $-\Delta$ is discrete with $\infty$ as the only accumulation point, and each eigenspace is of finite dimension.  Also, denote $(\lambda_{\epsilon,n,i}, u_{\epsilon, n,i})$ to be the $i$-th eigenpair of $\frac{I-T_{\textup{ref},\epsilon,n}}{\epsilon}$, where $\lambda_{\epsilon,n,i}$ is the $i$-th smallest eigenvalue.  We assume that both $u_{i}$ and $u_{\epsilon, n,i}$ are normalized in the $L^2$ norm.  With the above preparation, we are now ready to state the main theorem. 

\[spec\_cong\]\[Spectral convergence\] Fix $K\in \mathbb{N}$. Suppose the kernel is Gaussian; that is, $K_{\epsilon}(x,y)=e^{-{\left\lVertx-y\right\rVert}^2/\epsilon}$. Suppose $\lambda_{i}$ is simple.  Suppose $m=n^\beta$, where $\beta\in (0,1)$, $\epsilon=\epsilon(n)$ so that $\epsilon \to 0$ and $\frac{\sqrt{-\log\epsilon}+\sqrt{\log m}}{\sqrt{m}\epsilon^{d}} \rightarrow 0$, as $n \rightarrow \infty$, and $\sqrt{\epsilon} \leq \mathcal{K}_1 \min \Bigg(\bigg(\frac{\min(\mathsf\Gamma_K,1)}{\mathcal{K}_2+\lambda_K^{d/2+5}}\bigg)^2, \frac{1}{(2+\lambda_K^{d+1})^2}\Bigg)$, where $\mathsf\Gamma_K$, $\mathcal{K}_1$ and $\mathcal{K}_2>1$ are introduced in Proposition \[T epsilon and Delta\].  Then, there exists a sequence of signs $\{a_{n}\}$ such that with probability $1-\mathcal{O}(m^{-2})$, for all $i<K$, we have $$\begin{aligned}
{\left\lVerta_{n}u_{\epsilon,n,i}-u_{i}\right\rVert}_{L^\infty}&=\mathcal{O}(\epsilon^{1/2})+\mathcal{O}\left(\frac{\sqrt{-\log \epsilon}+\sqrt{\log m}}{\sqrt{m}\epsilon^{2d+2}}\right)\,,\\
{\left |\lambda_{\epsilon,n,i}-\lambda_{i}\right |}&=\mathcal{O}(\epsilon^{3/4})+\mathcal{O}\left(\frac{\sqrt{-\log \epsilon}+\sqrt{\log m}}{\sqrt{m}\epsilon^{2d+2}}\right)\nonumber\,,\end{aligned}$$ where the implied constants depend on the kernel, the curvature of $M$, $p_X$ and $p_Y$

We mention that the obtained convergence rate is just the worst case bound, and it should not be the optimal bound.  In fact, it is much more slowly than what we observed numerically.  Also, as we will show below, it seems that the eigenvector convergence should be faster than the eigenvalue convergence, but this is not reflected by the above spectral convergence rate.  How to obtain the “correct” convergence rate is however out of the scope of this paper, and will be explored in our future work. 

Numerical Results {#numerical}
=================

To illustrate how the Roseland performs, in addition to showing the dimension reduction and geometric recovery results, we also compare the results with the Nyström extension and the HKC.  For a fair comparison, in all the following simulations, the subset used in the Nyström extension and the reference set used in the HKC to embed the dataset are the same as the landmark set used in the Roseland.  As a result, the ranks of the matrices associated with the Roseland, the HKC and the Nyström extension are the same. 

Scalability of the Roseland – dimensional reduction and visualization
---------------------------------------------------------------------

In this subsection, we compare performance of the Nyström method and the Roseland from the aspect of dimension reduction and visualization. 

We take the dataset consists of random projections of the two-dimensional Shepp-Logan phantom.  A phantom is a $2$-dim image function $\psi$ compactly supported on $\mathbb{R}^{2}$. It is commonly applied in medical imaging society as a benchmark.  See Figure \[Figure:phantom\] for an example.  We generate a high dimensional data set by taking the Radon transform of $\psi$ followed by discretizing the projection image into $p\in \mathbb{N}$ points, and the projection angle is sampled uniformly from $S^{1}$. Without any symmetry assumption, the generated data set $\mathcal{X}\subset \mathbb{R}^p$ is a nonuniform sampling set on $R_\phi(S^1)\subset \mathbb{R}^p$, where $R_\phi$ is the Radon transform of $\phi$. In this simulation, we fix the number of discretization points $p=128$ and let the number of projections $n$ vary.  We run the DM, the Roseland, the HKC and the Nyström extension with $n=10,000$ and $m=n^\beta$, where $\beta=0.5$, and show the 3-dim embedding of $\mathcal{X}$ in Figure \[Figure:phantom\].  Clearly, both the DM and the Roseland recover the $S^1$ structure, while the Roseland is distorted.  On the other hand, the HKC and the Nyström extension are confused and lead to erraneous embeddings.  The computational times of different algorithms with $\beta=0.3$ and $\beta=0.5$ are also shown for a comparison.  Note that when $n=1,280,000$, the Roseland can finish in about 2.5 minutes in an ordinary laptop.  More visualization results can be found in Section \[More numerical results torus\]. 

![ The random projection data from the phantom image, where the data size is $n=10,000$ and the dimension is $p=128$. In the top row, we take $\beta =0.5$ for the Roseland, the Nyström, and the HKC.  In the bottom row, we show the relationship between the computational time and the data size.  Top row, from left to right: the DM embedding, the Roseland, the HKC, and the Nyström extension.  All embeddings are 3-dim, and have been rotated to optimize the visualization.  Bottom row: the runtime comparison of various algorithms when the data size varies, where $\beta=0.3$ on the left and $\beta=0.5$ on the right.  The x-axis is by the natural log, and the largest database size is $1,280,000$.\[Figure:phantom\]](phan10000DM_3d.jpg "fig:"){width="3cm"} ![ The random projection data from the phantom image, where the data size is $n=10,000$ and the dimension is $p=128$. In the top row, we take $\beta =0.5$ for the Roseland, the Nyström, and the HKC.  In the bottom row, we show the relationship between the computational time and the data size.  Top row, from left to right: the DM embedding, the Roseland, the HKC, and the Nyström extension.  All embeddings are 3-dim, and have been rotated to optimize the visualization.  Bottom row: the runtime comparison of various algorithms when the data size varies, where $\beta=0.3$ on the left and $\beta=0.5$ on the right.  The x-axis is by the natural log, and the largest database size is $1,280,000$.\[Figure:phantom\]](phan10000ref_3d.jpg "fig:"){width="3cm"} ![ The random projection data from the phantom image, where the data size is $n=10,000$ and the dimension is $p=128$. In the top row, we take $\beta =0.5$ for the Roseland, the Nyström, and the HKC.  In the bottom row, we show the relationship between the computational time and the data size.  Top row, from left to right: the DM embedding, the Roseland, the HKC, and the Nyström extension.  All embeddings are 3-dim, and have been rotated to optimize the visualization.  Bottom row: the runtime comparison of various algorithms when the data size varies, where $\beta=0.3$ on the left and $\beta=0.5$ on the right.  The x-axis is by the natural log, and the largest database size is $1,280,000$.\[Figure:phantom\]](phan10000haddad_3d.jpg "fig:"){width="3cm"} ![ The random projection data from the phantom image, where the data size is $n=10,000$ and the dimension is $p=128$. In the top row, we take $\beta =0.5$ for the Roseland, the Nyström, and the HKC.  In the bottom row, we show the relationship between the computational time and the data size.  Top row, from left to right: the DM embedding, the Roseland, the HKC, and the Nyström extension.  All embeddings are 3-dim, and have been rotated to optimize the visualization.  Bottom row: the runtime comparison of various algorithms when the data size varies, where $\beta=0.3$ on the left and $\beta=0.5$ on the right.  The x-axis is by the natural log, and the largest database size is $1,280,000$.\[Figure:phantom\]](phan10000nys_3d.jpg "fig:"){width="3cm"}\
![ The random projection data from the phantom image, where the data size is $n=10,000$ and the dimension is $p=128$. In the top row, we take $\beta =0.5$ for the Roseland, the Nyström, and the HKC.  In the bottom row, we show the relationship between the computational time and the data size.  Top row, from left to right: the DM embedding, the Roseland, the HKC, and the Nyström extension.  All embeddings are 3-dim, and have been rotated to optimize the visualization.  Bottom row: the runtime comparison of various algorithms when the data size varies, where $\beta=0.3$ on the left and $\beta=0.5$ on the right.  The x-axis is by the natural log, and the largest database size is $1,280,000$.\[Figure:phantom\]](runtime_log_beta03.pdf "fig:"){width="6cm"} ![ The random projection data from the phantom image, where the data size is $n=10,000$ and the dimension is $p=128$. In the top row, we take $\beta =0.5$ for the Roseland, the Nyström, and the HKC.  In the bottom row, we show the relationship between the computational time and the data size.  Top row, from left to right: the DM embedding, the Roseland, the HKC, and the Nyström extension.  All embeddings are 3-dim, and have been rotated to optimize the visualization.  Bottom row: the runtime comparison of various algorithms when the data size varies, where $\beta=0.3$ on the left and $\beta=0.5$ on the right.  The x-axis is by the natural log, and the largest database size is $1,280,000$.\[Figure:phantom\]](runtime_log_beta05.pdf "fig:"){width="6cm"}

Robustness of the Roseland – Geometric structure recovery {#numerical section noisy data and subeset}
---------------------------------------------------------

In this subsection, we show the robustness of the Roseland.  We compare performance of the Nyström extension, the HKC, and the Roseland from the aspect of recovering the geometric structure.  We consider the standard $S^1$ model, which is the one-dimensional canonical $S^1$ embedded in the first two coordinates of $\mathbb{R}^{100}$, since all ground truths can be analytically calculated.  Specifically, we uniformly sample $n=90,000$ points from the $S^1$ to be the dataset and independently sample another $m=3 00$ points uniformly to be the landmark set; that is, $\beta=0.5$. Then, embed all points to $\mathbb{R}^{p}$, where $p=100$, and add independent Gaussian noise $\epsilon_{i}$ to both the dataset and the landmark set, where $\epsilon_{i}$ are i.i.d.  sampled from $\mathcal{N}(0, \frac{1}{\sqrt{p}} I_{p\times p})$. For a fair comparison, we follow the following criterion used throughout this subsection. 

-   For the eigenvalues, we compare them with the ground truth and calculate the relative error.  Note that the non-trivial eigenvalues of the Laplace Beltromi operator of $S^{1}$ are $-k^{2}$, where $k=1,2,\ldots$, each with multiplicity $2$.

-   For the eigenvectors, note that the eigenfunctions related to the eigenvalue $-k^2$ are $\sin{k\pi\theta}$ and $\cos{k\pi\theta}$. Due to the non-trivial multiplicity, the eigen-decompsition is free up to rotations and reflections.  We carry out the following procedure:

    -   normalize all eigenvectors to have the unit $L^{2}$ norm. 

    -   For each positive odd $i$, the $i^{\textup{th}}$ and $(i+1)^{\textup{th}}$ non-trivial eigenvectors, denoted as $\textbf{u}_{i}$ and $\textbf{u}_{i+1}$ respectively, form a pair, because they correspond to the same true eigenvalue with multiplicity $2$. We will align them with $\sin{(i+1)\theta/2}$ and $\cos{(i+1)\theta/2}$ by searching for the phase shift $\phi$ and $a,b\in\{\pm 1\}$ that minimizes the following objective function: $${\left\lVerta\textbf{u}_{i} - \sin{(i+1)(\theta+\phi)/2}\right\rVert}_{\infty} + {\left\lVertb\textbf{u}_{i+1} - \cos{(i+1)(\theta+\phi)/2}\right\rVert}_{\infty}.$$ Denote $a^{\ast}$, $b^{\ast}$, and $\phi^{\ast}$ be the solutions. 

    -   Use $a^{\ast}, b^{\ast}, \phi^{\ast}$ to compute relative $L^{\infty}$ and $L^{2}$ errors between each eigenvector and the aligned ground truth eigenfunctions. 

-   For visualization purpose, we also plot the top 12 non-trivial eigenvectors of both methods and superimposed with the ground truth. 

Since the data sampling is uniform, for the HKC we apply the algorithm in the form shown in [@haddad2014texture Proposition 3.1]. 

The visualization results are shown in Figures \[S1noisy-noisy\_embedding\].  Clearly, while the Nyström extension and the HKC embed $S^1$ successfully, the embedding by the Roseland is cleaner.  Next, the recovered eigenvectors are shown in Figures \[S1noisy-noisy\_eigenvector\] and \[S1noisy-noisy\_eigenvalue\].  Clearly, the Nyström method can only successfully recover first few eigenfunctions of the Laplace-Beltrami operator (visually, only the first 8 look reasonably), while the HKC and the Roseland can recover more eigenfunctions (visually, the first 10 are reasonably well).  Also, note that there seems to be a sudden change, or “transition phenomenon”, from the 10-th to the 11-th eigenfunctions of the HKC and the Roseland, which cannot be explained by our current theory.  We will explore it in our future work.  Also note that “visually” the qualities of the first two non-trivial eigenfunctions of the Nyström extension and Roseland are similar, and they are both better than those of the HKC, but the qualities of embeddings are different.  To understand this discrepancy, we consider the following quantities.  Note that the first two non-trivial eigenvectors, $v_1,v_2\in \mathbb{R}^n$ from either the Nyström extension, the HKC or Roseland, if successfully recovered the eigenfunctions of the Laplace-Beltrami operator, should be $\sin(\theta+\phi)$ and $\cos(\theta+\phi)$ for some $\phi\in(0,2\pi]$ respectively.  Here, the phase $\phi$ comes from the uncertainty nature of the spectral embedding methods.  We then plot $\arctan(v_1(i)/v_2(i))$ and $\sqrt{v_1(i)^2+v_2(i)^2}$ against $\theta_i$, where $\theta_i$ is the angle of the $i$-th sampled point.  The results are shown in Figure \[S1noisy-noisy\_eigenvector\_phase\].  Clearly, the amplitude eigenvectors of the Nyström extension and the HKC fluctuates more than those of the Roseland, while the phase recovery qualities are similar.  Finally, we further quantify the difference by evaluating the relative error of eigenvalue and eigenfunction.  The results are shown in Figure \[S1noisy-noisy\_eigenvalue\].  Since the HKC is not designed for this purpose, it is not clear how to recover the eigenvalues.  So, we do not take it into consideration.  It is as expected that recovering the eigenvalues is challenging.  It is limited by the bandwidth we choose, since the largest eigenvalues we can recover is bounded by the inverse of the bandwidth.  However, the Roseland leads to a smaller relative error in eigenfunctions recovery.  Overall, this result indicates the potential of the proposed Roseland algorithm. 

More numerical results, including the clean dataset and landmark sets, and the noisy dataset and clean landmark set, are shown in Section \[More numerical results\]. 

![ Left: noisy data and noisy subset (only the first two coordinately are shown).  Middle left: then Nyström embedding.  Middle right: the HKC embedding.  Right: the Roseland embedding.  \[S1noisy-noisy\_embedding\]](noisy_S1_noisy_subset.jpg "fig:"){width="3cm"} ![ Left: noisy data and noisy subset (only the first two coordinately are shown).  Middle left: then Nyström embedding.  Middle right: the HKC embedding.  Right: the Roseland embedding.  \[S1noisy-noisy\_embedding\]](noisy_S1_noisy_subset_nys.jpg "fig:"){width="3cm"} ![ Left: noisy data and noisy subset (only the first two coordinately are shown).  Middle left: then Nyström embedding.  Middle right: the HKC embedding.  Right: the Roseland embedding.  \[S1noisy-noisy\_embedding\]](noisy_S1_noisy_subset_haddad.jpg "fig:"){width="3cm"}![Left: noisy data and noisy subset (only the first two coordinately are shown).  Middle left: then Nyström embedding.  Middle right: the HKC embedding.  Right: the Roseland embedding.  \[S1noisy-noisy\_embedding\]](noisy_S1_noisy_subset_ref.jpg "fig:"){width="3cm"}

![ Noisy data set and noisy subset.  Superimpose the top 12 non-trivial eigenvectors by the Nyström, the HKC and the Roseland with the ground truth (superimposed in red).  Top three rows: the top 6 eigenvectors; bottom three rows: the $7^{\textup{th}}$ to the $12^{\textup{th}}$ eigenvectors.  \[S1noisy-noisy\_eigenvector\]](noisyS1_noisy_sub_eigvec1-6.jpg "fig:"){width="14cm"} ![ Noisy data set and noisy subset.  Superimpose the top 12 non-trivial eigenvectors by the Nyström, the HKC and the Roseland with the ground truth (superimposed in red).  Top three rows: the top 6 eigenvectors; bottom three rows: the $7^{\textup{th}}$ to the $12^{\textup{th}}$ eigenvectors.  \[S1noisy-noisy\_eigenvector\]](noisyS1_noisy_sub_eigvec7-12.jpg "fig:"){width="14cm"}

![ Noisy data set and noisy subset.  Top left: the phase of the embedding by the Roseland.  Top middle: the phase of the embedding by the HKC.  Top right: the phase of the embedding by the Nyström extension.  Bottom left: the amplitude of the embedding by the Roseland and the Nyström.  Bottom right: the amplitude of the embedding by the Roseland and the HKC.  The phase of the embedding is determined by $\arctan(v_1(i)/v_2(i))$ and the amplitude is determined by $\sqrt{v_1(i)^2+v_2(i)^2}$, where $v_1$ and $v_2$ are the first non-trivial eigenvectors determined by the Nyström extension or Roseland.  \[S1noisy-noisy\_eigenvector\_phase\]](noisyS1_noisy_sub_nys_phase.jpg "fig:"){width="4cm"}![Noisy data set and noisy subset.  Top left: the phase of the embedding by the Roseland.  Top middle: the phase of the embedding by the HKC.  Top right: the phase of the embedding by the Nyström extension.  Bottom left: the amplitude of the embedding by the Roseland and the Nyström.  Bottom right: the amplitude of the embedding by the Roseland and the HKC.  The phase of the embedding is determined by $\arctan(v_1(i)/v_2(i))$ and the amplitude is determined by $\sqrt{v_1(i)^2+v_2(i)^2}$, where $v_1$ and $v_2$ are the first non-trivial eigenvectors determined by the Nyström extension or Roseland.  \[S1noisy-noisy\_eigenvector\_phase\]](noisyS1_noisy_sub_haddad_phase.jpg "fig:"){width="4cm"} ![ Noisy data set and noisy subset.  Top left: the phase of the embedding by the Roseland.  Top middle: the phase of the embedding by the HKC.  Top right: the phase of the embedding by the Nyström extension.  Bottom left: the amplitude of the embedding by the Roseland and the Nyström.  Bottom right: the amplitude of the embedding by the Roseland and the HKC.  The phase of the embedding is determined by $\arctan(v_1(i)/v_2(i))$ and the amplitude is determined by $\sqrt{v_1(i)^2+v_2(i)^2}$, where $v_1$ and $v_2$ are the first non-trivial eigenvectors determined by the Nyström extension or Roseland.  \[S1noisy-noisy\_eigenvector\_phase\]](noisyS1_noisy_sub_ref_phase.jpg "fig:"){width="4cm"}\
![ Noisy data set and noisy subset.  Top left: the phase of the embedding by the Roseland.  Top middle: the phase of the embedding by the HKC.  Top right: the phase of the embedding by the Nyström extension.  Bottom left: the amplitude of the embedding by the Roseland and the Nyström.  Bottom right: the amplitude of the embedding by the Roseland and the HKC.  The phase of the embedding is determined by $\arctan(v_1(i)/v_2(i))$ and the amplitude is determined by $\sqrt{v_1(i)^2+v_2(i)^2}$, where $v_1$ and $v_2$ are the first non-trivial eigenvectors determined by the Nyström extension or Roseland.  \[S1noisy-noisy\_eigenvector\_phase\]](noisyS1_noisy_sub_nys_ref_amplitude.pdf "fig:"){width="4cm"} ![ Noisy data set and noisy subset.  Top left: the phase of the embedding by the Roseland.  Top middle: the phase of the embedding by the HKC.  Top right: the phase of the embedding by the Nyström extension.  Bottom left: the amplitude of the embedding by the Roseland and the Nyström.  Bottom right: the amplitude of the embedding by the Roseland and the HKC.  The phase of the embedding is determined by $\arctan(v_1(i)/v_2(i))$ and the amplitude is determined by $\sqrt{v_1(i)^2+v_2(i)^2}$, where $v_1$ and $v_2$ are the first non-trivial eigenvectors determined by the Nyström extension or Roseland.  \[S1noisy-noisy\_eigenvector\_phase\]](noisyS1_noisy_sub_Haddad_ref_amplitude.pdf "fig:"){width="4cm"}

![ Illustration of the Nyström method and the Roseland on the noisy dataset and noisy landmark set.  Left: the top 18 non-trivial eigenvalues by the Nyström and Roseland with the ground truth.  Middle: relative error of eigenvalues.  Right: relative $L^{\infty}$ and $L^{2}$ error of the top 12 non-trivial eigenvectors by the Nyström and Roseland with the ground truth.  \[S1noisy-noisy\_eigenvalue\]](noisyS1_noisy_sub_eigval.pdf "fig:"){width="3.9cm"} ![ Illustration of the Nyström method and the Roseland on the noisy dataset and noisy landmark set.  Left: the top 18 non-trivial eigenvalues by the Nyström and Roseland with the ground truth.  Middle: relative error of eigenvalues.  Right: relative $L^{\infty}$ and $L^{2}$ error of the top 12 non-trivial eigenvectors by the Nyström and Roseland with the ground truth.  \[S1noisy-noisy\_eigenvalue\]](noisyS1_noisy_sub_eigval_err.pdf "fig:"){width="4cm"} ![ Illustration of the Nyström method and the Roseland on the noisy dataset and noisy landmark set.  Left: the top 18 non-trivial eigenvalues by the Nyström and Roseland with the ground truth.  Middle: relative error of eigenvalues.  Right: relative $L^{\infty}$ and $L^{2}$ error of the top 12 non-trivial eigenvectors by the Nyström and Roseland with the ground truth.  \[S1noisy-noisy\_eigenvalue\]](noisyS1_noisy_sub_err.pdf "fig:"){width="4.3cm"}

Geometric structure recovery – Geodesic distance estimation
-----------------------------------------------------------

We now show that if we want to recover the geodesic distance by the DM, then the Roseland outperforms the Nyström method.  Note that since the HKC is not designed for this purpose, we do not compare it here.  First, we describe the comparison methodology. 

-   Uniformly and independently sample $2,500$ points from $S^{1}$ as the dataset.  Uniformly and independently sample another $m=n^{\beta}$ points from $S^{1}$ as the landmark set for the Roseland.  In this experiment, we chose $m=50$ and $m=300$, where $\beta=0.5$ and $\beta\approx 0.72$ respectively. 

-   Fix $K\in \mathbb{N}$. Denote $\{x_{i}\}_{i=1}^{2,500}$ to be the dataset, and denote $\{y_{i}\}_{i=1}^{2,500}\subset \mathbb{R}^l$ to be the embedded dataset, where $l\in \mathbb{N}$ is the dimension of the Roseland.  Denote $x_i^{(K)}$ (resp.  $y_i^{(K)}$) to be the $K$-th nearest neighbor of $x_i$ (resp.  $y_i$).  The relative errors of the geodesic distance between $x_i$ and its $K$-th nearest neighbor is calculated by $$\frac{\big|D_t(x_{i},x_{i}^{(K)}) - d(x_{i},x_{i}^{(K)})\big|}{d(x_{i},x_{i}^{(K)})}\,,$$ where $d(x_{i},x_{i}^{(K)})$ is the ground truth geodesic distance between $x_{i}$ and $x_{i}^{(K)}$, and $D_t(x_{i},x_{i}^{(K)})$ can be the DD determined by the DM, the Roseland, or the Nyström extension. 

-   We compare the relative errors of the geodesic distance between $x_i$ and its $K$-th nearest neighbor by embeddings from the DM, the Nyström method and the Roseland. 

The results are shown in Figures \[geo\_recov\_300\] and \[geo\_recov\_50\].  Clearly, the geodesic distance can be well recovered by the Roseland, but the Nyström method is limited in doing so, in particularly when the landmark set (or $\beta$) is small.  This result is more dramatic if we use the ground truth eigenvalues to estimate the DD.  This result indicates that it is critical to recover the eigenvalues. 

The result is not surprising from the theoretical standpoint.  Note that it is shown in [@portegies2016embeddings] that one can obtain an almost isometric embedding of the manifold as long as one has enough number of eigenvalues and eigenfunctions of the Laplace-Beltromi operator.  Specifically, let $\varepsilon>0$ be any tolerable error given, then there exists some $t_{0}$, which depends on the manifold’s intrinsic dimension, Ricci curvature, injectivity radius and $\varepsilon$ such that for all $0<t<t_{0}$, there is a $N_{E}$ depends on the manifold’s intrinsic dimension, Ricci curvature, injectivity radius, volume, $\varepsilon$ and $t$ such that if $N>N_{E}$, the spectral embedding $$\begin{aligned}
x\mapsto 2t^{(d+2)/4}\sqrt{2}(4\pi)^{d/4} \begin{bmatrix}e^{-\lambda_{1}t}\phi_{1}(x)&\ldots& e^{-\lambda_{N}t}\phi_{N}(x)\end{bmatrix}^\top\end{aligned}$$ is almost isometric with the error controlled by $\varepsilon$, where $\lambda_{i}$ and $\phi_{i}$ is the $i^{\textup{th}}$ eigenpair of the Laplace Beltromi operator of the manifold.  The above theorem essentially says that we need to recover enough eigenfunctions if we want to have an accurate geodesic distance estimate by the DD.  Since the diffusion structure is not taken into account in the Nyström extension, it is limited in recovering higher order eigenfunctions.  On the other hand, the Roseland preserves the diffusion property, and hence the eigenfunctions of the Laplace-Beltrami operator. 

![ Subset size = 300, so $\beta= 0.5$. Left: relative errors of geodesic recovery by the DM, the Nyström extension and the Roseland using their own eigenpairs.  Right: relative errors of geodesic recovery by the DM, the Nyström extension and the Roseland using their own eigenvectors and the ground truth eigenvalues.\[geo\_recov\_300\]](geo_recover_sub300_err.pdf "fig:"){width="6cm"} ![ Subset size = 300, so $\beta= 0.5$. Left: relative errors of geodesic recovery by the DM, the Nyström extension and the Roseland using their own eigenpairs.  Right: relative errors of geodesic recovery by the DM, the Nyström extension and the Roseland using their own eigenvectors and the ground truth eigenvalues.\[geo\_recov\_300\]](geo_recover_sub300_err_true_eigenval.pdf "fig:"){width="6cm"}

![ Subset size = 50, so $\beta\approx 0.435$. Left: relative errors of geodesic recovery by the DM, the Nyström extension and the Roseland using their own eigenpairs.  Right: relative errors of geodesic recovery by the DM, the Nyström extension and the Roseland using their own eigenvectors and the ground truth eigenvalues.  \[geo\_recov\_50\]](geo_recover_sub50_err.pdf "fig:"){width="6cm"} ![ Subset size = 50, so $\beta\approx 0.435$. Left: relative errors of geodesic recovery by the DM, the Nyström extension and the Roseland using their own eigenpairs.  Right: relative errors of geodesic recovery by the DM, the Nyström extension and the Roseland using their own eigenvectors and the ground truth eigenvalues.  \[geo\_recov\_50\]](geo_recover_sub50_err_true_eigenval.pdf "fig:"){width="6cm"}

Control non-uniform sampling by designing the landmark set
----------------------------------------------------------

Recall Remark \[remark: dm\_thm\] – in the Roseland, if we can design the sampling scheme for the landmark set so that $\frac{2\nabla p_{X}(x)}{p_{X}(x)}+\frac{\nabla p_{Y}(x)}{p_{Y}(x)}=0$, then we remove the impact of the non-uniformly sampling and recover the Laplace-Beltrami operator.  The condition $\frac{2\nabla p_{X}(x)}{p_{X}(x)}+\frac{\nabla p_{Y}(x)}{p_{Y}(x)}=0$ suggests that we may want to sample the landmark set following the density function $p_{Y}(x)\propto\frac{1}{p^{2}_{X}(x)}$.

To illustrate this fact, we use a dataset non-uniformly sampled from the canonical $S^{1}$ as an example.  See Figure \[Figure:nonuniform-design\] for the result.  It is clear that if we do not design the landmark set, the first two nontrivial eigenfunctions estimated from the Roseland are deviated from the ground truth.  However, if the landmark set is well designed according to the developed theory, or could incorporate the background knowledge (like the HKC algorithm designed for the texture separation problem), we may better recover the desired ground truth.  This preliminary result warrants a further consideration of this “design problem” in our future work. 

![ Top row: the left subplot shows the true p.d.f.  $p(\theta)$ we use to sample points from $S^{1}$, which is parametrized by $\theta\in [0,2\pi)$, the left middle subplot shows the histogram of the sampled data, the right middle subplot shows the designed p.d.f.  $q(\theta)$ that is proportional to $1/p^2(\theta)$, and the right subplot shows the histogram of the landmark set sampled $\theta$ according to $q(\theta)$. Bottom row: the left subplot is the embedding determined by the Roseland, where the landmark set is chosen randomly, the left middle subplot shows the scattering plot of the recovered angles from the Roseland with the randomly chosen landmarks, where the x-axis is the estimated angles of $S^1$, and the y-axis is the associated true angles, the right middle subplot is the embedding determined by the Roseland with the designed landmark set following the law of $q(\theta)$, and the right subplot is the scattering plot of the recovered angles from the Roseland with the designed landmark set, where the x-axis is the estimated angles of $S^1$, and the y-axis is the associated true angles.  \[Figure:nonuniform-design\]](pdf_data.pdf "fig:"){width="3cm"} ![ Top row: the left subplot shows the true p.d.f.  $p(\theta)$ we use to sample points from $S^{1}$, which is parametrized by $\theta\in [0,2\pi)$, the left middle subplot shows the histogram of the sampled data, the right middle subplot shows the designed p.d.f.  $q(\theta)$ that is proportional to $1/p^2(\theta)$, and the right subplot shows the histogram of the landmark set sampled $\theta$ according to $q(\theta)$. Bottom row: the left subplot is the embedding determined by the Roseland, where the landmark set is chosen randomly, the left middle subplot shows the scattering plot of the recovered angles from the Roseland with the randomly chosen landmarks, where the x-axis is the estimated angles of $S^1$, and the y-axis is the associated true angles, the right middle subplot is the embedding determined by the Roseland with the designed landmark set following the law of $q(\theta)$, and the right subplot is the scattering plot of the recovered angles from the Roseland with the designed landmark set, where the x-axis is the estimated angles of $S^1$, and the y-axis is the associated true angles.  \[Figure:nonuniform-design\]](sample_pdf_data.pdf "fig:"){width="3cm"} ![ Top row: the left subplot shows the true p.d.f.  $p(\theta)$ we use to sample points from $S^{1}$, which is parametrized by $\theta\in [0,2\pi)$, the left middle subplot shows the histogram of the sampled data, the right middle subplot shows the designed p.d.f.  $q(\theta)$ that is proportional to $1/p^2(\theta)$, and the right subplot shows the histogram of the landmark set sampled $\theta$ according to $q(\theta)$. Bottom row: the left subplot is the embedding determined by the Roseland, where the landmark set is chosen randomly, the left middle subplot shows the scattering plot of the recovered angles from the Roseland with the randomly chosen landmarks, where the x-axis is the estimated angles of $S^1$, and the y-axis is the associated true angles, the right middle subplot is the embedding determined by the Roseland with the designed landmark set following the law of $q(\theta)$, and the right subplot is the scattering plot of the recovered angles from the Roseland with the designed landmark set, where the x-axis is the estimated angles of $S^1$, and the y-axis is the associated true angles.  \[Figure:nonuniform-design\]](pdf_ref.pdf "fig:"){width="3cm"} ![ Top row: the left subplot shows the true p.d.f.  $p(\theta)$ we use to sample points from $S^{1}$, which is parametrized by $\theta\in [0,2\pi)$, the left middle subplot shows the histogram of the sampled data, the right middle subplot shows the designed p.d.f.  $q(\theta)$ that is proportional to $1/p^2(\theta)$, and the right subplot shows the histogram of the landmark set sampled $\theta$ according to $q(\theta)$. Bottom row: the left subplot is the embedding determined by the Roseland, where the landmark set is chosen randomly, the left middle subplot shows the scattering plot of the recovered angles from the Roseland with the randomly chosen landmarks, where the x-axis is the estimated angles of $S^1$, and the y-axis is the associated true angles, the right middle subplot is the embedding determined by the Roseland with the designed landmark set following the law of $q(\theta)$, and the right subplot is the scattering plot of the recovered angles from the Roseland with the designed landmark set, where the x-axis is the estimated angles of $S^1$, and the y-axis is the associated true angles.  \[Figure:nonuniform-design\]](sample_pdf_ref.pdf "fig:"){width="3cm"} ![ Top row: the left subplot shows the true p.d.f.  $p(\theta)$ we use to sample points from $S^{1}$, which is parametrized by $\theta\in [0,2\pi)$, the left middle subplot shows the histogram of the sampled data, the right middle subplot shows the designed p.d.f.  $q(\theta)$ that is proportional to $1/p^2(\theta)$, and the right subplot shows the histogram of the landmark set sampled $\theta$ according to $q(\theta)$. Bottom row: the left subplot is the embedding determined by the Roseland, where the landmark set is chosen randomly, the left middle subplot shows the scattering plot of the recovered angles from the Roseland with the randomly chosen landmarks, where the x-axis is the estimated angles of $S^1$, and the y-axis is the associated true angles, the right middle subplot is the embedding determined by the Roseland with the designed landmark set following the law of $q(\theta)$, and the right subplot is the scattering plot of the recovered angles from the Roseland with the designed landmark set, where the x-axis is the estimated angles of $S^1$, and the y-axis is the associated true angles.  \[Figure:nonuniform-design\]](embed_rand_ref.jpg "fig:"){width="3cm"} ![ Top row: the left subplot shows the true p.d.f.  $p(\theta)$ we use to sample points from $S^{1}$, which is parametrized by $\theta\in [0,2\pi)$, the left middle subplot shows the histogram of the sampled data, the right middle subplot shows the designed p.d.f.  $q(\theta)$ that is proportional to $1/p^2(\theta)$, and the right subplot shows the histogram of the landmark set sampled $\theta$ according to $q(\theta)$. Bottom row: the left subplot is the embedding determined by the Roseland, where the landmark set is chosen randomly, the left middle subplot shows the scattering plot of the recovered angles from the Roseland with the randomly chosen landmarks, where the x-axis is the estimated angles of $S^1$, and the y-axis is the associated true angles, the right middle subplot is the embedding determined by the Roseland with the designed landmark set following the law of $q(\theta)$, and the right subplot is the scattering plot of the recovered angles from the Roseland with the designed landmark set, where the x-axis is the estimated angles of $S^1$, and the y-axis is the associated true angles.  \[Figure:nonuniform-design\]](angle_rand_ref.jpg "fig:"){width="3cm"} ![ Top row: the left subplot shows the true p.d.f.  $p(\theta)$ we use to sample points from $S^{1}$, which is parametrized by $\theta\in [0,2\pi)$, the left middle subplot shows the histogram of the sampled data, the right middle subplot shows the designed p.d.f.  $q(\theta)$ that is proportional to $1/p^2(\theta)$, and the right subplot shows the histogram of the landmark set sampled $\theta$ according to $q(\theta)$. Bottom row: the left subplot is the embedding determined by the Roseland, where the landmark set is chosen randomly, the left middle subplot shows the scattering plot of the recovered angles from the Roseland with the randomly chosen landmarks, where the x-axis is the estimated angles of $S^1$, and the y-axis is the associated true angles, the right middle subplot is the embedding determined by the Roseland with the designed landmark set following the law of $q(\theta)$, and the right subplot is the scattering plot of the recovered angles from the Roseland with the designed landmark set, where the x-axis is the estimated angles of $S^1$, and the y-axis is the associated true angles.  \[Figure:nonuniform-design\]](embed_design_ref.pdf "fig:"){width="3cm"} ![ Top row: the left subplot shows the true p.d.f.  $p(\theta)$ we use to sample points from $S^{1}$, which is parametrized by $\theta\in [0,2\pi)$, the left middle subplot shows the histogram of the sampled data, the right middle subplot shows the designed p.d.f.  $q(\theta)$ that is proportional to $1/p^2(\theta)$, and the right subplot shows the histogram of the landmark set sampled $\theta$ according to $q(\theta)$. Bottom row: the left subplot is the embedding determined by the Roseland, where the landmark set is chosen randomly, the left middle subplot shows the scattering plot of the recovered angles from the Roseland with the randomly chosen landmarks, where the x-axis is the estimated angles of $S^1$, and the y-axis is the associated true angles, the right middle subplot is the embedding determined by the Roseland with the designed landmark set following the law of $q(\theta)$, and the right subplot is the scattering plot of the recovered angles from the Roseland with the designed landmark set, where the x-axis is the estimated angles of $S^1$, and the y-axis is the associated true angles.  \[Figure:nonuniform-design\]](angle_design_ref.jpg "fig:"){width="3cm"}

Discussion and conclusion {#discuss}
=========================

In this paper, we introduce a new algorithm based on the landmark set to accelerate the DM algorithm.  In addition to providing a series of theoretical justification, we also provide a series of numerical examples to support the potential of the algorithm. 

Relationship with the alternating diffusion algorithms
------------------------------------------------------

Note that the Roseland is related to the recently developed alternating DM (ADM) algorithm [@lederman2015learning].  The ADM is developed to deal with the sensor fusion problem, when we have multiple data-sets simultaneously acquired by multimodal sensors.  In short, suppose we have 2 aligned data sets $\mathcal{X}^{(1)}=\{x_{i}^{(1)}\}_{i=1}^{n}$ and $\mathcal{X}^{(2)}=\{x_{i}^{(2)}\}_{i=1}^{n}$, then we build two transition matrices $K^{(j)}=D^{-1}W$ from $\mathcal{X}^{(j)}$, where $j=1,2$. Next the alternating diffusion operator is defined as $K^{(1)}K^{(2)}$. It is easy to check that $K^{(1)}K^{(2)}$ is also row stochastic, so it can be considered as a transition probability matrix of a new Markov chain that alternates between the two data sets.  The idea of the Roseland is closely related to the ADM, in the sense that we can consider the landmark set as a new data set on its own right.  So, the idea of measuring similarity between two data points via the landmark set can be understood as diffusing between the original data set and the landmark set from the ADM point of view. 

Application of the Roseland
---------------------------

The idea of landmark set considered in the Roseland have several applications.  Here we mention two of them.  The VDM [@singer2012vector] is a generalization of DM that aims to encode the group structure when comparing objects.  The VDM suffers from the expensive computational cost more than the DM, since the group structure is usually represented as a matrix, which inflates the matrix size.  Specifically, if the group structure is represented as a $q\times q$ matrix and we have $n$ objects to compare, then we need to eigendecompse a $nq\times nq$ kernel matrix in the VDM.  We remark that the idea of the landmark set can be generalized to the VDM, and hence for the application.  We will report an application of generalizing the landmark set to accelerate VDM in the future paper.  Spectral clustering methods are known to perform well when the classical clustering methods such as $k$-means and linkage fail [@alzate2008multiway].  In real applications, however, we usually run into the multiway spectral clustering problem.  It is well known that the more clusters we need to determine, the more eigenvectors we need [@alzate2008multiway; @lee2014multiway].  As is shown in the numerical section, the Roseland has the ability to recover more and better eigenvectors, at least compared with the Nystöm extension.  This shows the potential of applying the Roseland for the multiway spectral clustering purpose. 

Optimal variance control {#Discussion optimal variance convergence discrepancy}
------------------------

Note that in the pointwise convergence, the variance is bounded by the large deviation, while the bound might not be the optimal one.  We discuss this problem from a theoretical aspect and a numerical aspect.  Recall the definition of the U-statistics in , which is a special case of (\[dep\_sum\]).  Hoeffding [@hoeffding1994probability] proved that for all $t>0$, we have the following bound for the $r$-degree U-statistics: $$\mathbb{P}(|U_{r}-\mathbb{E}(U_{r})|\geq t)\leq 2\textup{exp}\left(\frac{-nt^{2}}{r}\right).$$ Later, when $\sigma=\textup{Var}(h(X_{1},\ldots,X_{r}))$ is finite, Arcones [@arcones1995bernstein] refined it to a Bernstein-like inequalities; that is, for all $t>0$, $$\mathbb{P}(|U_{r}-\mathbb{E}(U_{r})|\geq t)\leq a\textup{exp}\left(\frac{-(n/r)t^{2}}{2\sigma^{2}+bt}\right),$$ where $a,b$ are some constants.  Note that this bound is not better than the one shown in Theorem \[bernstein\] [@janson2004large].  We comment that since the U-statistic is a special form of equation (\[dep\_sum\]), we would expect the order of the convergence rate provided in Theorems \[hoeffding\] and \[bernstein\] be good enough since they are of the same order as those for the U-statistics.  While showing the optimal bound, either order or constant, for the U-statistics or the more general form like (\[dep\_sum\]) is not the focus of this paper, we could numerically evaluate the quality of the bound provided in Theorems \[bernstein\] [@janson2004large].  See Figure \[largedev\] for a simulation, where we observe that the empirical convergence rate of $M\times N$ grid samples is “faster” than $M$ i.i.d.  samples.  This at least suggests that we could expect a better constant in front of the convergence order term.  This theoretical work is critical to further understand the algorithm, and we will explore this large deviation rate in our future work. 

![ We run two simulations to evaluate the optimality of the convergence rate and plot the relative errors.  In both cases, we let $M=N^{1/2}$ and compare the convergence of the empirical mean using i.i.d.  $M$ samples, i.i.d.  $N$ samples, and $N\times M$ grid samples.  Left: $X,Y\sim U[0,6]$ and $f(X,Y)=XY$. Right: $X,Y\sim U[0,1]$ and $f(X,Y)=100e^{-100[(Y-0.5)^{2}+(X-Y)^{2}]}$.[]{data-label="largedev"}](large_dev1.pdf "fig:"){width="6.2cm"} ![ We run two simulations to evaluate the optimality of the convergence rate and plot the relative errors.  In both cases, we let $M=N^{1/2}$ and compare the convergence of the empirical mean using i.i.d.  $M$ samples, i.i.d.  $N$ samples, and $N\times M$ grid samples.  Left: $X,Y\sim U[0,6]$ and $f(X,Y)=XY$. Right: $X,Y\sim U[0,1]$ and $f(X,Y)=100e^{-100[(Y-0.5)^{2}+(X-Y)^{2}]}$.[]{data-label="largedev"}](large_dev2.pdf "fig:"){width="6.2cm"}

Another natural issue we shall further elaborate is the discrepancy raised in Remark \[Remark about the variance convergence rate\].  If we traced the proof of the variance analysis stated in Theorem \[variancethm\], we know that the “slower” convergence rate incurred by the worse variance $\mathcal{O}\big(\frac{\sqrt{\log(n)}}{n^{\beta/2}\epsilon^{d/2-1/2}}\big)$ comes from Step 2 in the proof.  Specifically, the variance shown in \[Proof Variance Theorem key variance quantity\], which is $\epsilon^{-d/2}$ worse than that of the GL.  As discussed in Remark \[Remark about the variance convergence rate\], even if $\beta=1$; that is, we have the same number of independent points in the landmark set, the variance is still worse than that of the ordinary GL according to our theory.  This counterintuitive result comes from the fact that the “kernel” cannot be well approximated.  To be more specific, note that we would expect the quantity $\frac{1}{m}\sum_{k=1}^{m}\epsilon^{-d} K_{\epsilon}(x_{i},y_{k})K_{\epsilon}(y_{k},x_{j})\approx \epsilon^{-d} \int_M K_{\epsilon}(x_{i},y)K_{\epsilon}(y,x_{j})p_Y(y)dV(y)=:\epsilon^{-d/2} \tilde{K}_{\epsilon}(x_i,x_j)$ when $m$ is big, like $m=n$, where $\tilde{K}$ is another kernel.  If this approximation is accurate, then we are reduced to the GL case.  However, by a direct calculation, we see that this approximation is not accurate with a non-negligible variance.  It is this non-negligible error that leads to the worse variance in the Roseland algorithm.  In Figure \[m=n figures\], we show a quantitative result for a comparison.  It is clear that even if $m=n=2,500$, the performance of the Roseland is worse than that of the DM.  This clearly indicates that when $m$ is large, like $m=n$, the Roseland does not gain any benefit. 

![ The comparison of the DM with $n=2,500$ uniform samples from $S^1$ and the Roseland with the same $n$ samples and $m=n$ landmark sets.  Left: the relative errors of eigenvalues; right: the relative errors or eigenfunctions.  Note that even when $m=n$, the performance of the Roseland is worse than the DM.[]{data-label="m=n figures"}](cleanS1_clean_sub_m=n_err.pdf "fig:"){width="6.2cm"} ![ The comparison of the DM with $n=2,500$ uniform samples from $S^1$ and the Roseland with the same $n$ samples and $m=n$ landmark sets.  Left: the relative errors of eigenvalues; right: the relative errors or eigenfunctions.  Note that even when $m=n$, the performance of the Roseland is worse than the DM.[]{data-label="m=n figures"}](cleanS1_clean_sub_m=n_eigval_err.pdf "fig:"){width="6.2cm"}

Optimal spectral convergence rate
---------------------------------

Another relevant topic we need to discuss is the spectral convergence.  To the best of our knowledge, there have been several papers studying the spectral convergence rate, for example [@von2008consistency; @wang2015spectral; @shi2015convergence; @trillos2018error; @calder2019improved], but the “optimal” spectral convergence rate is not yet known.  Among these papers, our approach and setup are closest to [@von2008consistency; @wang2015spectral], where we studied the connection between the graph laplacian and the heat kernel, and utilized the covering number of the kernel function.  In [@shi2015convergence], the authors consider a different setup.  In [@trillos2018error; @calder2019improved], the setup is similar, but the approach is different.  The optimal transport approach considered in [@trillos2018error] and an improvement in [@calder2019improved] both have a potential to be combined with the analysis strategy considered in this paper.  See [@DunsonWuWu2019] for more discussion.  It is interesting to ask if we are able to determine the “optimal” spectral convergence rate under the manifold setup, and design an algorithm to achieve it.  Last but not the least, we shall mention that without the knowledge of “optimal spectral convergence rate”, we are not able to select bandwidth to recover the spectral structure.  In our numerical simulation, we provide a fair comparison based on a naive bandwidth setup.  However, we consistently find that if we select a smaller bandwidth for the Roseland, then the eigenvalue reconstruction is better, while the Nystöm extension fails.  To sum up, this topic is critical to further understand the algorithm and for the purpose of statistical inference, and we will report our exploration in our future work. 

Acknowledgement
===============

Chao Shen thanks Xiucai Ding and Nan Wu for various discussion of the topic. 

Technical background for the proof {#lemma}
==================================

In this section we prepare some known technical lemmas and necessary results for our proof.  To be succinct, detailed proofs will be skipped and we refer readers to the relative papers. 

Some basic differential geometry facts
--------------------------------------

The following Lemma is critical for us to control the size of the eigenfunction. 

\[lemma hormander\] Fix a compact Riemannian manifold $(M,g)$. For the $l$-th pair of eigenvalue $\lambda_l$ and eigenfunction $\phi_l$ of the Laplace-Beltrami operator, where $\|\phi_l\|_2=1$, we have $$\begin{aligned}
\|\phi_l\|_\infty \leq  C_1 \lambda_l^{\frac{d-1}{4}}\,,   \nonumber \end{aligned}$$ where $C_1$ is a constant depending on the injectivity radius and sectional curvature of the manifold $M$.

The following lemma is the well-known Weyl’s law, which controls the eigenvalue growth. 

\[lemma weyls law\] Fix a compact and connected Riemannian manifold $(M,g)$. The eigenvalues of the Laplace-Beltrami operator, denoted as $0=\lambda_1< \lambda_2\leq\ldots$, satisfy $$\begin{aligned}
l=C_2 \lambda^{d/2}_l+O\big(\lambda^{\frac{d-1}{2}}_l \log \lambda_l\big)\,,  \nonumber \end{aligned}$$ where $C_2$ is a constant depending on the volume of the manifold. 

The proof of the following truncation lemma can be found in, for example, [@Coifman2006]. 

\[trunc\] Suppose $f\in L^{\infty}(M^{d})$ and $0<\gamma<1/2$. Then for any $x\in M^{d}$, when $\epsilon$ is sufficiently small, we have $$\left|\int_{M^{d}\backslash\widetilde{\mathcal{B}}_{\epsilon^{\gamma}}(x)}\epsilon^{-d/2}K_{\epsilon}(x,y)f(y)\,dV(y)\right|=\mathcal{O}(\epsilon^{2})\,,$$ where $\mathcal{O}(\epsilon^{2})$ depends on $\|f\|_{\infty}$, and $\widetilde{\mathcal{B}}_{\epsilon^{\gamma}}(x)\coloneqq\iota^{-1}(\mathcal{B}_{\epsilon^{\gamma}}(x)\cap\iota(M^{d}))\subseteq M^{d}$, where $\mathcal{B}_{\epsilon^{\gamma}}(x)$ is the Euclidean ball with radius $\epsilon^{\gamma}$ centered at $x$.

The proof of the following approximation of identity type lemma follows the standard argument, and can be found in [@singer2016spectral Lemma B.3]. 

\[for\_numerator\_est\] If $f\in C^{4}(M^{d})$, then for all $x\in M^{d}$, we have $$\int_{M^{d}}\epsilon^{-d/2}K_{\epsilon}(x,y)f(y)\,dV(y)=f(x)+\frac{\epsilon\mu_{1,2}^{(0)}}{2d}(\Delta f(x)-w(x)f(x))+\mathcal{O}(\epsilon^{2})\,,$$ where $w(x)=\frac{1}{3}s(x)-\frac{d}{12|S^{d-1}|}\frac{\mu_{1,3}^{(1)}}{\mu_{1,2}^{(0)}}\int_{S^{d-1}}{\textup{I}\!\textup{I}}_{x}^{2}(\theta,\theta)\,d\theta$, $s(x)$ is scalar curvature at $x$, ${\textup{I}\!\textup{I}}_{x}$ is the second fundamental form of the embedding at $x$, and $|S^{d-1}|$ is the volume of the canonical $(d-1)$-sphere. 

\[numerator\_est\] If $f\in C^{4}(M^{d})$, then for all $x\in M^{d}$, we have $$\begin{aligned}
    \nonumber &\int_{M^{d}}\epsilon^{-d/2}K_{\epsilon}(x,y)f(y)p_{X}(y)\,dV(y) 
    \\\nonumber =&\, p_{X}(x)\left[1+\frac{\epsilon\mu_{1,2}^{(0)}}{d}\left(\frac{\Delta p_{X}(x)}{2p_{X}(x)} - \frac{1}{2}w(x)\right)f(x)\right]
    \\\nonumber &\qquad+ \epsilon\frac{\mu_{1,2}^{(0)}p_{X}(x)}{d}\left(\frac{\Delta f(x)}{2}+\frac{\nabla f(x)\cdot\nabla p_{X}(x)}{p_{X}(x)}\right) + \mathcal{O}(\epsilon^{2})\end{aligned}$$ where $w(x)$ as in lemma \[for\_numerator\_est\]. 

It is an immediate consequence of lemma \[for\_numerator\_est\] by replacing $f(y)$ by $f(y)p_{X}(y)$.

Glivenko-Cantelli class and entropy bound
-----------------------------------------

A set $\mathfrak{F}\subseteq C(M)$ is called a Glivenko-Cantelli class if $$\textup{sup}_{f\in\mathfrak{F}}|\mathbb{P}_{n}f - \mathbb{P}f|\longrightarrow 0\quad\textup{a.s.}$$

To handle the randomness, we need to control the complexity of “all possible” functions that are related to the random samples.  The following definition is the quantity we need. 

Let $(\mathcal{F},{\left\lVert\cdot\right\rVert})$ be a subset of normed space of real functions $f:M\rightarrow\mathbb{R}$. Given two functions $l$ and $u$, the bracket $[l,u]$ is the set of all functions $f$ such that $l(x)\leq f(x)\leq u(x)$ for all $x\in M$. An $\epsilon$-bracket is a bracket such that ${\left\lVertu-l\right\rVert}<\epsilon$. The bracketing number $N_{[]}(\epsilon,\mathcal{F},{\left\lVert\cdot\right\rVert})$ is the minimum number of $\epsilon$-bracket needed to cover $\mathcal{F}$. The upper and lower bounds $u$ and $l$ of the brackets need not belong to $\mathcal{F}$ but are assumed to have finite norms. 

To proceed to the spectral convergence rate, we need the following results that control the rate of convergence from finite sample points.  To this end, we need the following entropy bound [@von2008consistency Theorem 19].  Or see [@mendelson2003few] for a more systematic review of the topic. 

\[entropy\_bd\]\[Entropy bound [@mendelson2003few; @von2008consistency]\].  Let $(\mathcal{X},\mathcal{B},\mathbb{P})$ be an arbitrary probability space, $\mathcal{F}$ a class of real-valued functions on $\mathcal{X}$ with ${\left\lVertf\right\rVert}_{\infty}\leq 1$. Let $(X_{n})$ be a sequence of i.i.d.  random variables drawn from $\mathbb{P}$, and $(\mathbb{P}_{n})$ the corresponding empirical distributions.  For $\delta>0$, there exists a constant $C_E>0$ such that, for all $ n\in\mathbb{N}$, with probability higher than $1-\delta$: $$\sup_{F\in\mathcal{F}}{\left |\mathbb{P}_{n}F-\mathbb{P}F\right |}\leq\frac{C_E}{\sqrt{n}}\int_{0}^{\infty}\sqrt{\log N(\mathcal{F},r,L_{2}(\mathbb{P}_{n}))}\,dr + \sqrt{\frac{1}{n}\log\frac{2}{\delta}}.$$

The following lemma is the standard statement about the covering number when the kernel is Gaussian.  Since the proof can be found in [@wang2015spectral Lemma 4.1], we omit it here. 

\[cover\_num\_gau\](Covering number of Gaussian, [@wang2015spectral Lemma 4.1]).  Take $\epsilon>0$ and $k(x,y):=K_\epsilon(x,y)=e^{-{\left\lVertx-y\right\rVert}^2/\epsilon}$ in Definition \[Definition different classes\].  For $r>0$, we have the following bound: $$N(\mathcal{K},r,{\left\lVert\cdot\right\rVert}_{\infty})\leq\Big(\frac{24\sqrt{2d}D_{M}}{r\epsilon}\Big)^{2d},$$ where $D_{M}$ is the diameter of $M^{d}$.

Facts we need for the spectral convergence
------------------------------------------

Next we need the notion of collectively compact convergence.  Recall that $C(M)$ is a Banach space with the ${\left\lVert\cdot\right\rVert}_{\infty}$ norm. 

([@chatelin2011spectral p.122]) Let $(E, \|\cdot\|_{E})$ be an arbitrary Banach space.  A sequence of operators $T_{n}:E\rightarrow E$ converges to $T:E\rightarrow E$ collectively compactly if and only if the following conditions are satisfied:

(I) $T_{n}\rightarrow T$ pointwisely,

(II) the set $\cup_{n}(T_{n}-T)B$ is relatively compact in $E$, where $B$ is the closed unit ball centered at $0$ in $E$.

We state the following theorem for the spectral convergence. 

Let $(E,{\left\lVert\cdot\right\rVert}_{E})$ be any Banach space, $\{T_{n}\}_{n}$ and $T$ be bounded linear operators on $E$ such that $T_{n}\rightarrow T$ compactly.  Let $\lambda\in\sigma(T)$ be an isolated eigenvalue with finite multiplicity $m$, and $\mathcal N_\lambda\subset \mathbb{C}$ be an open neighborhood of $\lambda$ such that $\sigma(T)\cap \mathcal N_\lambda=\lambda$. Then:

(a) (Convergence of eigenvalues) There exists an $N\in\mathbb{N}$ such that for all $n>N$, the set $\sigma(T_{n})\cap \mathcal N_\lambda$ is an isolated part of $\sigma(T_{n})$ consisting of at most $m$ different eigenvalues, and their multiplicities sum to $m$. Moreover, the sequence of sets $\{\sigma(T_{n})\cap \mathcal N_\lambda\}\rightarrow\{\lambda\}$ in the sense that every sequence $\{\lambda_{n}\}$ with $\lambda_{n}\in\sigma(T_{n})\cap \mathcal N_\lambda$ satisfies $\lim\lambda_{n}\rightarrow\lambda$.

(b) (Convergence of spectral projections) Let $P_\lambda$ be the spectral projection of $T$ corresponding to $\lambda$. Let $P_{\lambda,n}$ be the spectral projection of $T_{n}$ corresponding to $\sigma(T_{n})\cap \mathcal N_\lambda$, then $P_{\lambda,n}\rightarrow P_\lambda$ pointwisely. 

Note that when $\lambda$ is simple with the eigenfunction $f$, this theorem can be simplified.  Indeed, there exists $N\in\mathbb{N}$ such that for all $n>N$, the sets $\sigma(T_{n})\cap \mathcal N_\lambda$ consists of a simple eigenvalue $\lambda_{n}$, and $\lim\lambda_{n}\rightarrow\lambda$. Moreover, for the corresponding eigenfunctions $f_{n}$, there exists a sequence $a_{n}\in \{1,\,-1\}$ so that ${\left\lVerta_{n}f_{n}-f\right\rVert}_{E}\rightarrow 0$.

Recall that the $L^2(M)$ space is a separable Hilbert space.  The following lemma is the key toward the spectral convergence rate [@DunsonWuWu2019].  The proof can be found in [@DunsonWuWu2019], so we omit it here.  We mention that part (c) of this lemma is also considered in [@calder2019improved] to improve the $L^2$ convergence rate. 

\[ratio of eigenvalues\] Let $A$ and $B$ be two compact self-adjoint operators from the separable Hilbert space $H$ to $H$. Let $(\cdot, \cdot)$ be the inner product of $H$. Suppose the eigenvalues of $A$, denoted as $\lambda_l(A)$, $l=1,\ldots$, are simple and positive, and the eigenvalues of $B$, denoted as $\lambda_l(B)$, $l=1,\ldots$, are simple and bounded from below so that $1=\lambda_1(A) > \lambda_2(A) > \cdots\geq 0$ and $\lambda_1(B) >  \lambda_2(B) > \cdots $.  Denote $\{u_i\}$ to be the orthonormal eigenfunctions of $A$ and $\{w_i\}$ to be orthonormal eigenfunctions of $B$. Furthermore, denote $$\gamma_i(B):=\min(\lambda_i(B)-\lambda_{i-1}(B), \,\lambda_{i+1}(B)-\lambda_i(B))\,.\label{Proof proposition 1 Definition gamma i}$$ Let $E: = A-B$. Then, for $\epsilon>0$ we have the following statements:

(a) If $|\frac{(E f,f)}{(Af,f)}| \leq \epsilon$ for all $f \in L^2$, then for all $i$, we have $$\left|\frac{1-\lambda_i(B)}{1-\lambda_i(A)}-1\right|\leq \epsilon\,.\nonumber$$

(b) If $\|Bu_i-\lambda_i(B)u_i\|_2 \leq \epsilon$, then for $a=1$ or $-1$, we have $$\|aw_i- u_i\|_2 \leq \frac{2\epsilon}{\gamma_i(B)}.  \nonumber$$ Moreover, $$|(u_i,w_i)| \geq 1-\frac{\epsilon}{\gamma_i(B)}\,.\nonumber$$

(c) The eigenvalues satisfy $$|\lambda_i(A)-\lambda_j(B)|\leq \frac{\|Ew_j\|_2}{|(u_i,w_j)|}\,.\nonumber$$

The following result describes how the spectral convergence, or more precisely, the eigenfunction convergence, happens when a sequence of operators converges.  It is a restatement of Equation (5) in Theorem 3 in [@atkinson1967numerical].  We refer the readers with interest in its derivation to [@DunsonWuWu2019 Theorem SI.1]. 

\[atkinson1967bound\]\[[@atkinson1967numerical], or [@DunsonWuWu2019 Theorem SI.1]\] Let $(E, \|\cdot\|_{E})$ be an arbitrary Banach space.  Let $\{T_n\}_{n=1}^\infty$ and $T$ be compact linear operators on $E$ such that $\{T_n\}_{n=1}^\infty$ converges to $T$ collectively compactly.  For a nonzero eigenvalue $\lambda \in \sigma(T)$, denote the corresponding spectral projection by $\textup{Pr}_{\lambda}$. Let $D\subset \mathbb{C}$ be an open neighborhood of $\lambda$ such that $\sigma(T)\cap D$ =$\{ \lambda\}$. There exists some $N\in \mathbb{N}$ such that for all $n > N$, $\sigma(T_n)\cap D=\{\lambda_n\}$. Let $\textup{Pr}_{\lambda_n}$ be the corresponding spectral projections of $T_n$ for $\lambda_n$. Let $r <|\lambda|$ and $r<\textup{dist}(\{\lambda\}, \sigma(K)\setminus \{\lambda\})$. Then, for every $x \in \textup{Pr}_\lambda(E)$, we have $$\|x-\textup{Pr}_{\lambda_n} x\|_{E} \leq \max_{z\in \Gamma_r(\lambda)} \frac{2r\|R_{z}(T)\|}{\min_{z\in \Gamma_r(\lambda)}|z|}\left(\|(T_n-T)x\|_{E} +\|R_{z}(T)x\|_{E}\|(T-T_n)T_n\|\right)\,,
\nonumber$$ where $\Gamma_r(\lambda):=\{z\in \mathbb{C}|\,|z-\lambda|=r\}$.

We need the following lemma to connect the ultimate eigenfunction convergence to Theorem \[atkinson1967bound\]. 

\[cong\_1\_dim\_proj\]\[convergence of one-dim projections [@von2008consistency Proposition 18]\] Let $\{v_{n}\}_{n=1}^\infty$ be a sequence of vectors in a Banach space $(E,{\left\lVert\cdot\right\rVert})$ with ${\left\lVertv_{n}\right\rVert}=1$. Denote $\textup{Pr}_{v_n}$ be the projection onto the one-dimensional subspaces spanned by $v_{n}$. Take $v\in E$ with ${\left\lVertv\right\rVert}=1$. Then there exists a sequence of signs $(a_{n})$ such that $${\left\lVerta_{n}v_{n}-v\right\rVert}\leq2{\left\lVertv-\textup{Pr}_{v_n}(v)\right\rVert}.$$ In particular, if ${\left\lVertv-\textup{Pr}_{v_n}(v)\right\rVert}\rightarrow 0$, then $v_{n}\rightarrow v$ in $(E,{\left\lVert\cdot\right\rVert})$ up to a change of sign. 

  ---------------------------------------------------- -- ------------------------------------------------------------------------------------------------------------
                      $K_{\textup{ref},\epsilon}(x,y)$    $\displaystyle\int_{M}K_{\epsilon}(x,z)K_{\epsilon}(z,y)p_{Y}(z)\,dV(z)$
                        $d_{\textup{ref},\epsilon}(x)$    $\displaystyle\int_{M}K_{\textup{ref},\epsilon}(x,y)p_{X}(y)\,dV(y)$
                      $M_{\textup{ref},\epsilon}(x,y)$    $\displaystyle \frac{K_{\textup{ref},\epsilon}(x,y)}{d_{\textup{ref},\epsilon}(x)}$
                       $T_{\textup{ref},\epsilon}f(x)$    $\displaystyle\int_{M}M_{\textup{ref},\epsilon}(x,y)f(y)p_{X}(y)\,dV(y)$
          $\widehat{K}_{\textup{ref},\epsilon,n}(x,y)$    $\displaystyle\frac{1}{m}\sum_{j}^{m}K_{\epsilon}(x,z_{j})K_{\epsilon}(z_{j},y)$
            $\widehat{d}_{\textup{ref},\epsilon,n}(x)$    $\displaystyle\frac{1}{n}\sum_{i=1}^{n}\widehat{K}_{\textup{ref},\epsilon,n}(x,x_{i})$
          $\widehat{M}_{\textup{ref},\epsilon,n}(x,y)$    $\displaystyle\frac{\widehat{K}_{\textup{ref},\epsilon,n}(x,y)}{\widehat{d}_{\textup{ref},\epsilon,n}(x)}$
           $\widehat{T}_{\textup{ref},\epsilon,n}f(x)$    $\displaystyle\frac{1}{n}\sum_{i=1}^{n}\widehat{M}_{\textup{ref},\epsilon,n}(x,x_{i})f(x_{i})$
                      $d_{\textup{ref},\epsilon,n}(x)$    $\displaystyle\frac{1}{n}\sum_{i=1}^n K_{\textup{ref},\epsilon}(x,x_i)$
                     $T_{\textup{ref},\epsilon,n}f(x)$    $\displaystyle\frac{1}{n}\sum_{i=1}^{n}M_{\textup{ref},\epsilon}(x,x_{i})f(x_{i})$
    $\widehat{M}^{(d)}_{\textup{ref},\epsilon,n}(x,y)$    $\displaystyle\frac{K_{\textup{ref},\epsilon}(x,y)}{\widehat{d}_{\textup{ref},\epsilon,n}(x)}$
                                         $\mathbb{P}f$    $\displaystyle\int f(x)p_X(x)\,dV(x)$
                                     $\mathbb{P}_{n}f$    $\displaystyle\frac{1}{n}\sum_{i=1}^{n}f(x_{i})$
                             $\widetilde{\mathbb{P}}g$    $\displaystyle\int g(y)p_{Y}(y)\,dV(y)$
                           $\widetilde{\mathbb{P}}_mg$    $\displaystyle\frac{1}{m}\sum_{l=1}^m g(y_l)$
  ---------------------------------------------------- -- ------------------------------------------------------------------------------------------------------------

  : Table of notation throughout the proof.  Assume $f,g\in C(M)$.

\[tab:TableOfNotationForMyResearch\]

Proof of Theorems \[biasthm\] and \[variancethm\] – pointwise convergence {#biasproof}
=========================================================================

The main theoretical contribution of this paper is handling how the overall diffusion behaves when the diffusion must goes through the landmark set.  The result will form the base of the spectral convergence proof. 

\[Proof: Theorem1: bias analysis\] Take $f\in C^{3}(M^{d})$. Then, for all $x\in M^{d}$ we have $$\begin{aligned}
    T_{\textup{ref},\epsilon}f(x) - f(x) =&\, \frac{\epsilon\mu_{1,2}^{(0)}}{d}\left(\frac{2\nabla p_{X}(x)}{p_{X}(x)}+\frac{\nabla p_{Y}(x)}{p_{Y}(x)}\right)\cdot\nabla f(x) \\
    &+ \frac{\epsilon\mu_{1,2}^{(0)}}{d}\Delta f(x) + \mathcal{O}(\epsilon^{3/2})\,,\nonumber\end{aligned}$$ where the implied constant in $\mathcal{O}(\epsilon^{3/2})$ depends on the $C^3$ norm of $f$, the $C^2$ norms of $p_X$ and $p_Y$, and the Ricci curvature of the manifold. 

The bias analysis is almost the same as those shown in [@Coifman2006; @singer2016spectral], except the extra step handling the landmark set.  By Definition \[GL\_cts\], we have $$\begin{aligned}
    \nonumber T_{\textup{ref},\epsilon}f(x) &= \frac{\int_{M}K_{\textup{ref},\epsilon}(x,y)f(y)p_{X}(y)\,dV(y)}{\int_{M}K_{\textup{ref},\epsilon}(x,y)p_{X}(y)\,dV(y)}\,.\end{aligned}$$ We first compute the numerator, which satisfies $$\begin{aligned}
    \nonumber&\int_{M}K_{\textup{ref},\epsilon}(x,y)f(y)p_{X}(y)\,dV(y)
    \\\nonumber =& \int_{M}\left(\int_{M}K_{\epsilon}(x,z)K_{\epsilon}(z,y)p_{Y}(z)\,dV(z)\right)f(y)p_{X}(y)\,dV(y)
    \\\nonumber =& \int_{M}\left(\int_{M}K_{\epsilon}(z,y)f(y)p_{X}(y)\,dV(y)\right)K_{\epsilon}(x,z)p_{Y}(z)\,dV(z)\,.\end{aligned}$$ By Lemmas \[trunc\] and \[numerator\_est\], the right hand side can be expanded and organized as $$\begin{aligned}
\nonumber &\, \epsilon^{d}\left(p_{X}(x)p_{Y}(x)+\frac{\epsilon\mu_{1,2}^{(0)}}{d}p_{Y}(x)\Delta p_{X}(x) + \frac{\epsilon\mu_{1,2}^{(0)}}{2d}p_{X}(x)\Delta p_{Y}(x)\right. \\ 
    \nonumber&\left.-\frac{\epsilon\mu_{1,2}^{(0)}}{d}w(x)p_{X}(x)p_{Y}(x)+\frac{\epsilon\mu_{1,2}^{(0)}}{d}\nabla p_{X}(x)\cdot\nabla p_{Y}(x)\right)f(x)
    \\\nonumber &+ \epsilon^{d}\left(\frac{2\epsilon\mu_{1,2}^{(0)}}{d}p_{Y}(x)\nabla p_{X}(x)+\frac{\epsilon\mu_{1,2}^{(0)}}{d}p_{X}(x)\nabla p_{Y}(x)\right)\cdot\nabla f(x) \nonumber\\
    &+ \frac{\epsilon^{d+1}\mu_{1,2}^{(0)}}{d}p_{X}(x)p_{Y}(x)\Delta f(x)+\mathcal{O}(\epsilon^{d+2})\,.\end{aligned}$$ Next, note that the denominator is just the numerator with $f(x)$ replaced by the constant function $1$. Hence, we have: $$\begin{aligned}
    \nonumber &\int_{M}K_{\textup{ref},\epsilon}(x,y)p_{X}(y)\,dV(y)
    \\\nonumber =&\, \epsilon^{d}\left(p_{X}(x)p_{Y}(x)+\frac{\epsilon\mu_{1,2}^{(0)}}{d}p_{Y}(x)\Delta p_{X}(x) + \frac{\epsilon\mu_{1,2}^{(0)}}{2d}p_{X}(x)\Delta p_{Y}(x)\right. \\ 
    \nonumber&\left.\qquad-\frac{\epsilon\mu_{1,2}^{(0)}}{d}w(x)p_{X}(x)p_{Y}(x)+\frac{\epsilon\mu_{1,2}^{(0)}}{d}\nabla p_{X}(x)\cdot\nabla p_{Y}(x)\right)+\mathcal{O}(\epsilon^{d+2})\,.\end{aligned}$$ By putting them together, we have $$\begin{aligned}
    \nonumber T_{\textup{ref},\epsilon}f(x) &= f(x) + \frac{\epsilon\mu_{1,2}^{(0)}}{d}\left(\frac{2\nabla p_{X}(x)}{p_{X}(x)}+\frac{\nabla p_{Y}(x)}{p_{Y}(x)}\right)\cdot\nabla f(x) + \frac{\epsilon\mu_{1,2}^{(0)}}{d}\Delta f(x) + \mathcal{O}(\epsilon^{2})\,.\end{aligned}$$

Note the constants in front of $\frac{\nabla p_{X}(x)}{p_{X}(x)}$ and $\frac{\nabla p_{Y}(x)}{p_{Y}(x)}$ respectively.  Intuitively, the $2$ in front of $\frac{\nabla p_{X}(x)}{p_{X}(x)}$ comes from the “2” steps diffusion. 

The variance analysis is less trivial, and we need to apply the large deviation theorem when dependence exists. 

Take $\mathcal{X}=\{x_i\}_{i=1}^n$ and $\mathcal{Y}=\{y_j\}_{i=1}^m$, where $m=[n^{\beta}]$ for some $0<\beta\leq 1$ and $[x]$ is the nearest integer of $x\in\mathbb{R}$. Take $f\in C(M^{d})$ and denote $\textbf{\textit{f}}\in\mathbb{R}^{n}$ such that $\textbf{\textit{f}}_{i}=f(x_{i})$. Let $\epsilon = \epsilon(n)$ so that $\frac{\sqrt{\log n}}{n^{\beta/2}\epsilon^{d/2+1/2}}\rightarrow 0$ and $\epsilon\rightarrow 0$ when $n\rightarrow\infty$. Then with probability higher than $1 - \mathcal{O}(1/n^{2})$, we have $$\big[(I-(D^{(\textup{R})})^{-1}W^{(\textup{R})})\textbf{\textit{f}}\;\big](i) = f(x_{i}) - T_{\textup{ref},\epsilon}f(x_{i}) + \mathcal{O}\Big(\frac{\sqrt{\log n}}{n^{\beta/2}\epsilon^{d/2-1/2}}\Big)$$ for all $i=1,2,\ldots,n$.

Define $L:=I-(D^{(\textup{R})})^{-1}W^{(\textup{R})}$. Fix some $x_{i}\in M^{d}$, by Definition \[ref\_affinity\], $$\begin{aligned}
    \nonumber(L\textbf{\textit{f}})_{i} \,&=\frac{\sum_{j=1}^{n}W_{ij}(f(x_{i})-f(x_{j}))}{\sum_{j=1}^{n}W_{ij}}=\frac{\sum_{j=1}^{n}K_{\textup{ref},\epsilon}(x_{i},x_{j})(f(x_{i})-f(x_{j}))}{\sum_{j=1}^{n}K_{\textup{ref},\epsilon}(x_{i},x_{j})}
    \\ &= \label{grid_samp}
    \frac{\frac{1}{nm}\sum_{j=1,k=1}^{n,m}\epsilon^{-d} K_{\epsilon}(x_{i},y_{k})K_{\epsilon}(y_{k},x_{j})(f(x_{i})-f(x_{j}))}{\frac{1}{nm}\sum_{j=1,k=1}^{n,m}\epsilon^{-d}K_{\epsilon}(x_{i},y_{k})K_{\epsilon}(y_{k},x_{j})}\,.\end{aligned}$$ Define two random variables $$\begin{aligned}
\label{def_F}
    F\coloneqq\epsilon^{-d}K_{\epsilon}(x_{i},Y)K_{\epsilon}(Y,X)(f(x_{i})-f(X))\end{aligned}$$ $$\begin{aligned}
\label{def_G}
    G\coloneqq\epsilon^{-d}K_{\epsilon}(x_{i},Y)K_{\epsilon}(Y,X)\end{aligned}$$ Recall that the landmark set $\mathcal{Y}=\{y_{k}\}_{k=1}^{m}$ are i.i.d.  samples from the random vector $Y$, which has the p.d.f.  $p_{Y}$. Also, the data set $\mathcal{X}=\{x_{j}\}_{j=1}^{n}$ are i.i.d.  sampled from the random vector $X$, which has the p.d.f.  $p_{X}$. Moreover, $Y$ is independent of $X$. Denote by $F_{k,j}$ one realization of $F$ when the realization of the random vector $(X,Y)$ is $(x_j,\,y_{k})$; in other words, $F_{k,j}=\epsilon^{-d}K_{\epsilon}(x_{i},y_{k})K_{\epsilon}(y_{k},x_{j})(f(x_{i})-f(x_{j}))$. Similarly for $G_{k,j}$. Hence the numerator in can be written as a random variable $$\begin{aligned}
\label{F}
    \frac{1}{mn}\textbf{F}\coloneqq\frac{1}{mn}\sum_{j=1,k=1}^{n,m}F_{k,j}\end{aligned}$$ and the denominator can be written as a random variable $$\begin{aligned}
\label{G}
    \frac{1}{mn}\textbf{G}\coloneqq\frac{1}{mn}\sum_{j=1,k=1}^{n,m}G_{k,j}\end{aligned}$$ By the law of large number, we would expect that $$(L\textbf{\textit{f}})_{i}\approx\frac{\mathbb{E}(F)}{\mathbb{E}(G)}\,.$$ We now justify this intuition, and get the convergence rates of $$\begin{aligned}
\label{top_rate}
    \frac{1}{mn}\sum_{j=1,k=1}^{n,m}F_{k,j}\longrightarrow\mathbb{E}(F)\end{aligned}$$ and $$\begin{aligned}
    \frac{1}{mn}\sum_{j=1,k=1}^{n,m}G_{k,j}\longrightarrow\mathbb{E}(G)\,,\end{aligned}$$ and hence the convergence rate of $$(L\textbf{\textit{f}})_{i}\longrightarrow\frac{\mathbb{E}(F)}{\mathbb{E}(G)}\,.$$

We start by solving \[top\_rate\] and there are three steps.**\[Step1\]**.  We know $\chi(\mathcal{A})=\mathcal{O}(\max(m,n))$ by . 

**\[Step2\]**.  We want to compute $\textup{Var}(F)=\mathbb{E}(F^{2})-(\mathbb{E}(F))^{2}$. From Lemma \[numerator\_est\], we have $$\begin{aligned}
    \nonumber\mathbb{E}(F) =&\, \epsilon^{-d}\int_{M}\left(\int_{M}K_{\epsilon}(x_{i},y)K_{\epsilon}(y,x)p_{Y}(y)\,dV(y)\right)(f(x_{i})-f(x))p_{X}(x)\,dV(x)
    \\\nonumber =&\, \epsilon^{-d}\left(\frac{\epsilon^{d+1}\mu_{1,2}^{(0)}}{d}(2p_{Y}\nabla p_{X}+p_{X}\nabla p_{Y})\cdot\nabla ((f(x_{i})-f(x))p_{X}(x))|_{x=x_{i}}\right. \\ 
    \nonumber&\left.  +\frac{\epsilon^{d}\mu_{1,2}^{(0)}}{d}p_{X}p_{Y}\Delta ((f(x_{i})-f(x))p_{X}(x))|_{x=x_{i}}+\mathcal{O}(\epsilon^{d+2})\right)
    \\\nonumber =&\, \frac{\epsilon\mu_{1,2}^{(0)}}{d}(2p_{Y}\nabla p_{X}+p_{X}\nabla p_{Y})\cdot\nabla ((f(x_{i})-f(x))p_{X}(x))|_{x=x_{i}} \\
    \nonumber& +\frac{\epsilon\mu_{1,2}^{(0)}}{d}p_{X}p_{Y}\Delta ((f(x_{i})-f(x))p_{X}(x))|_{x=x_{i}}+\mathcal{O}(\epsilon^{2});\end{aligned}$$ Similarly, by applying Lemma \[numerator\_est\] twice, we get $$\begin{aligned}
    \label{Proof Variance Theorem key variance quantity}&\mathbb{E}(F^{2}) \\
    \nonumber=&\, \epsilon^{-2d}\int_{M}\left(\int_{M}K_{\epsilon}^{2}(x_{i},y)K_{\epsilon}^{2}(y,x)p_{Y}(y)\,dV(y)\right)\\
    &\qquad\qquad\times (f(x_{i})-f(x))^{2}p_{X}(x)\,dV(x)\nonumber
    \\\nonumber =&\, \int_{M}\epsilon^{-d}\left(\int_{M}\epsilon^{-d}K_{\epsilon}^{2}(y,x)(f(x_{i})-f(x))^{2}p_{X}(x)\,dV(x)\right)\\
    &\qquad\qquad\times K_{\epsilon}^{2}(x_{i},y)p_{Y}(y)\,dV(y)\nonumber
    \\\nonumber =&\, \frac{\epsilon^{1-d}\mu_{2,0}^{(0)}\mu_{2,2}^{(0)}}{d}\Delta((f(y)-f(x_{i}))^{2}p_{X}(y)p_{Y}(y))|_{y=x_{i}} + \mathcal{O}(\epsilon^{2-d})\,.\end{aligned}$$ Without loss of generality, we assume from now that $\Delta((f(y)-f(x_{i}))^{2}p_{X}(y)p_{Y}(y))|_{y=x_{i}}$ is positive.  Therefore, when $\epsilon>0$ is sufficiently small, $\mathbb{E}(F^{2})\asymp \epsilon^{1-d}$. Since $[\mathbb{E}(F)]^2 =\mathcal{O}(\epsilon^{2})$, we have $\mathbb{E}(F^{2}) \gg  (\mathbb{E}(F))^{2}$ and hence $\textup{Var}(F)\asymp \epsilon^{1-d}$.**\[Step3\]**.  We apply Theorem \[bernstein\] to establish the large deviation bound.  From **\[Step2\]** we have $\textup{Var}(F_{k,j})\asymp \epsilon^{1-d}$, which is controlled by $|F_{k,j}|\asymp \epsilon^{-d}$ by .  Moreover, $\mathbb{E}(F_{k,j})=\mathcal{O}(\epsilon)$, so $|F_{k,j}-\mathbb{E}(F_{k,j})|\leq b$, for some $b>0$ satisfying $b\asymp \epsilon^{-d}$ when $\epsilon$ is sufficiently small.  Hence, by Theorem \[bernstein\], let $m=[n^{\beta}]$ for any $0<\beta\leq 1$, we have for all $ t>0$: $$\begin{aligned}
    \nonumber\mathbb{P}\Big(\frac{1}{mn}\textbf{F}-\mathbb{E}(F)\geq t\Big)\leq&\,\textup{exp}\left(\frac{-8(mnt)^{2}}{25\chi(\mathcal{A})(\sum_{k,j}\textup{Var}(F_{k,j})+bmnt/3)}\right)
    \\\nonumber\asymp &\,\textup{exp}\left(\frac{-8(mnt)^{2}}{25(m+n)(c_1mn\epsilon^{1-d}+c_2mnt\epsilon^{-d}/3)}\right)
    \\\nonumber=&\,\textup{exp}\left(\frac{-8mnt^{2}}{25(m+n)(c_1\epsilon^{1-d}+c_2t\epsilon^{-d}/3)}\right)\,,\end{aligned}$$ where $c_1,c_2>0$ are implied constants in $\textup{Var}(F_{k,j})\asymp \epsilon^{1-d}$ and $b\asymp \epsilon^{-d}$ respectively.  Since our goal is to estimate the Laplace-Beltrami term, which has the prefactor of order $\epsilon$, we ask $\frac{t}{\epsilon}\rightarrow 0$. As a result, the exponent becomes $$\begin{aligned}
    \label{aa}\frac{8mnt^{2}}{25(m+n)(c_1\epsilon^{1-d}+c_2t\epsilon^{-d}/3)} \geq \frac{c_3n^{1+\beta}t^{2}}{(n^{\beta}+n)\epsilon^{1-d}} \geq\frac{c_3n^{1+\beta}t^{2}}{n\epsilon^{1-d}}=\frac{c_3n^{\beta}t^{2}}{\epsilon^{1-d}}\,,\end{aligned}$$ for some constant $c_3>0$. Then, if we choose $n$ such that $\frac{c_3n^{\beta}t^{2}}{\epsilon^{1-d}}=3\log n$, we have $$\begin{aligned}
    \nonumber t\asymp \frac{\sqrt{\log n}}{n^{\beta/2}\epsilon^{d/2-1/2}}\,,\end{aligned}$$ which satisfies the request that $\frac{t}{\epsilon} \asymp \frac{\sqrt{\log n}}{n^{\beta/2}\epsilon^{d/2+1/2}}\rightarrow 0$ as $n\rightarrow\infty$ by assumption.  As a result, by the chosen $n$, we have $$\begin{aligned}
    \nonumber\mathbb{P}\Big(\frac{1}{mn}\textbf{F}-\mathbb{E}(F)\geq t\Big)\leq\textup{exp}(-3\log n) = \frac{1}{n^{3}}\,.\end{aligned}$$ Recall that we have fixed $x_{i}$ for some $i$. In order for $i=1,2,\ldots,n$, we simply use the union bound to get: $$\begin{aligned}
    \nonumber\mathbb{P}\Big(\frac{1}{mn}\textbf{F}-\mathbb{E}(F)\geq t; \textup{for}\;i=1,2,\ldots,n\Big)\leq n\times\frac{1}{n^{3}}=\frac{1}{n^{2}}\,.\end{aligned}$$

In general, when $\Delta((f(y)-f(x_{i}))^{2}p_{X}(y)p_{Y}(y))|_{y=x_{i}}=0$, the variance is smaller, and $t$ is smaller.  Hence, we conclude that with probability $1-\mathcal{O}(n^{-2})$, for all $x_i$, the numerator of $(L\textbf{\textit{f}})_{i}$ equals $$\begin{aligned}
\label{num_LF}
    \mathbb{E}(F)+\mathcal{O}\left(\frac{\sqrt{\log n}}{n^{\alpha/2}\epsilon^{d/2-1/2}}\right)\,.\end{aligned}$$ The the denominator follows the same line and we list the computation here for the convenience of the reader: $$\begin{aligned}
    \nonumber\mathbb{E}(G)\, &= \epsilon^{-d}\int_{M}\left(\int_{M}K_{\epsilon}(x_{i},y)K_{\epsilon}(y,x)p_{Y}(y)\,dV(y)\right)p_{X}(x)\,dV(x)
    \\\nonumber &= p_{X}p_{Y}(x_{i}) + \mathcal{O}(\epsilon)\end{aligned}$$ and $$\begin{aligned}
    \nonumber\mathbb{E}(G^{2}) \,&= \epsilon^{-2d}\int_{M}\left(\int_{M}K_{\epsilon}^{2}(x_{i},y)K_{\epsilon}^{2}(y,x)p_{Y}(y)\,dV(y)\right)p_{X}(x)\,dV(x)
    \\\nonumber &= \epsilon^{-d}\mu_{2,0}^{(0)}\mu_{2,0}^{(0)}p_{X}p_{Y}(x_{i}) + \mathcal{O}(\epsilon^{1-d})\,.\end{aligned}$$ By the same argument, with probability $1-\mathcal{O}(1/n^{2})$, we have for all $x_i$ that the denominator of $(L\textbf{\textit{f}})_{i}$ satisfies $$\begin{aligned}
\label{denom_LF}
    \mathbb{E}(G)+\mathcal{O}\left(\frac{\sqrt{\log n}}{n^{\alpha/2}\epsilon^{d/2}}\right)\,.\end{aligned}$$ By combining and and the binomial expansion, we conclude that for all $i=1,2,\ldots,n$, with probability $1-\mathcal{O}(1/n^{2})$, we have $$\begin{aligned}
    \nonumber(L\textbf{\textit{f}})_{i} \,=& \frac{\mathbb{E}(F)+\mathcal{O}\left(\frac{\sqrt{\log n}}{n^{\alpha/2}\epsilon^{d/2-1/2}}\right)}{\mathbb{E}(G)+\mathcal{O}\left(\frac{\sqrt{\log n}}{n^{\alpha/2}\epsilon^{d/2}}\right)}
     = \frac{\mathbb{E}(F)}{\mathbb{E}(G)} + \mathcal{O}\left(\frac{\sqrt{\log n}}{n^{\alpha/2}\epsilon^{d/2-1/2}}\right)
    \\\nonumber =&
    f(x_{i}) - T_{\textup{ref},\epsilon}f(x_{i}) + \mathcal{O}\Big(\frac{\sqrt{\log n}}{n^{\alpha/2}\epsilon^{d/2-1/2}}\Big)\,,\end{aligned}$$ and hence the proof. 

Proof of Theorems \[spec\_cong\] – spectral convergence
=======================================================

We extend the argument provided in [@belkin2007convergence; @von2008consistency], and apply tools from [@DunsonWuWu2019] to prove the spectral convergence in the $L^\infty$ sense and its corresponding rate.  We first define some notations.  Recall the definition of $K_{\textup{ref},\epsilon}$ and $d_{\textup{ref},\epsilon}$ in Definition \[GL\_cts\].  Define the normalized landmark-kernel as $$M_{\textup{ref},\epsilon}(x,y)=\frac{K_{\textup{ref},\epsilon}(x,y)}{d_{\textup{ref},\epsilon}(x)}\in C(M\times M)\,.$$ Note that by this definition we have $$T_{\textup{ref},\epsilon}f(x)=\int_{M}M_{\textup{ref},\epsilon}(x,y)f(y)p_{X}(y)\,dV(y)\,.$$ When we only have finite sample points $\{x_{i}\}_{i=1}^{n}$ and landmark set $\{z_{i}\}_{i=1}^{m}$, we need to handle various terms in Definition \[Definition various operators discrete\], and the following “intermittent” term: $$\begin{aligned}
d_{\textup{ref},\epsilon,n}(x)&\,:=\frac{1}{n}\sum_{i=1}^n K_{\textup{ref},\epsilon}(x,x_i)\,,\\
T_{\textup{ref},\epsilon,n}f(x)&\,:=\frac{1}{n}\sum_{i=1}^{n}M_{\textup{ref},\epsilon}(x,x_{i})f(x_{i})\,.\nonumber\end{aligned}$$

The following lemma says that $T_{\textup{ref},\epsilon}$, $T_{\textup{ref},\epsilon,n}$, and $\widehat{T}_{\textup{ref},\epsilon,n}$ are all “nice” integral operators.  The proof is similar to that shown in [@belkin2007convergence; @von2008consistency], so we omit it. 

The integral operators $T_{\textup{ref},\epsilon}$, $T_{\textup{ref},\epsilon,n}$, and $\widehat{T}_{\textup{ref},\epsilon,n}$ are all compact. 

The proof of Theorems \[spec\_cong\] is composed of two major parts. 

-   Part 1.  When $p_{Y}$ is well chosen, show that $\frac{T_{\textup{ref},\epsilon}-1}{\epsilon}\rightarrow \frac{\mu_{1,2}^{(0)}}{d}\Delta$ “spectrally” as $\epsilon\rightarrow 0$, and evaluate the rate that depends on $\epsilon$.

-   Part 2.  For a fixed $\epsilon>0$, show that $\widehat{T}_{\textup{ref},\epsilon,n}\rightarrow T_{\textup{ref},\epsilon}$ compactly $\textup{a.s.}$ as $n\rightarrow\infty$, and evaluate the rate that depends on $\epsilon$ and $n$.

Below, we prepare needed facts for these two major parts.  With these facts, we put them together to finish the proof in the end. 

Facts for Part 1
----------------

We need the following Proposition.  This proposition is the key step toward the spectral convergence.  Its proof is long and delicate, and can be found in [@DunsonWuWu2019], we only provide key steps and refer readers with interest to [@DunsonWuWu2019 Proposition 1] for details. 

\[T epsilon and Delta\] Assume that all eigenvalues of $\Delta$ are simple.  Denote $(\lambda_{i,\epsilon}, \phi_{i,\epsilon})$ to be the $i$-th eigenpair of $\frac{I-T_{\textup{ref},\epsilon}}{\epsilon}$ and $(\lambda_{i}, \phi_{i})$ to be the $i$-th eigenpair of $-\Delta$. Assume both $\phi_{i,\epsilon}$ and $\phi_{i}$ are normalized in the $L^2$ norm.  Fix $K\in \mathbb{N}$. Denote $$\mathsf \Gamma_K:=\min_{1 \leq i \leq K}\textup{dist}(\lambda_i, \sigma(-\Delta)\setminus \{\lambda_i\})\,.$$ Suppose $\sqrt{\epsilon} \leq \mathcal{K}_1 \min \Bigg(\bigg(\frac{\min(\mathsf\Gamma_K,1)}{\mathcal{K}_2+\lambda_K^{d/2+5}}\bigg)^2, \frac{1}{(2+\lambda_K^{d+1})^2}\Bigg)$, where $\mathcal{K}_1$ and $\mathcal{K}_2>1$ are the constants depending on $p_X$, $p_Y$, and the volume, the injectivity radius and the sectional curvature of the manifold.  Furthermore, assume $p_Y$ is properly chosen so that $\frac{2\nabla p_{X}(x)}{p_{X}(x)}+\frac{\nabla p_{Y}(x)}{p_{Y}(x)}=0$. Then, there are $a_i \in \{-1, 1\}$ such that for all $i < K$, $$\begin{aligned}
|\lambda_{i,\epsilon}-\lambda_{i}| &\, \leq \epsilon^{4/3} \,,  \label{Proposition main result statement} \\
\|a_i\phi_{i,\epsilon}-\phi_{i}\|_{\infty} &\, \leq \epsilon^{1/2} \,.  \nonumber \end{aligned}$$

Note that the kernel associated with $T_{\textup{ref},\epsilon}$ is $$\begin{aligned}
K_{\textup{ref},\epsilon}(x,y)&\coloneqq\int_{M}K_{\epsilon}(x,z)K_{\epsilon}(z,y)p_{Y}(z)\,dV(z)\,.\end{aligned}$$ While in general $K_{\textup{ref},\epsilon}$ is not Gaussian, it is smooth and decays exponentially fast.  To proceed, note that by Theorem \[Proof: Theorem1: bias analysis\], when $p_Y$ is properly chosen so that $\frac{2\nabla p_{X}(x)}{p_{X}(x)}+\frac{\nabla p_{Y}(x)}{p_{Y}(x)}=0$, we have the pointwise convergence of the eigenvalue/eigenfunction of $\frac{1-T_{\textup{ref},\epsilon}}{\epsilon}$ to those of $-\Delta$. On the other hand, if we plug the eigenfunction of $-\Delta$ into Theorem \[Proof: Theorem1: bias analysis\], the error in the pointwise convergence depends on the $C^4$ norm of the eigenfunction.  Therefore, by the standard Sobolev embedding (see [@palais1968foundations Theorem 9.2] or [@DunsonWuWu2019 Lemma SI.8]), the error of the pointwise convergence is controlled in the uniform way.  Therefore, by plugging the kernel $K_{\textup{ref},\epsilon}$ into the proof of [@DunsonWuWu2019 Proposition 1], where Lemma \[ratio of eigenvalues\] is applied to control the deviations of eigenvalues and eigenfunctions, we obtain the result.  Note that the implied constants associated with error bounds are different from those shown in [@DunsonWuWu2019 Proposition 1] due to the different kernels we choose here.  Also, note that the bandwidth used in [@DunsonWuWu2019] is $\epsilon$, while it is $\epsilon^{1/2}$ in this work. 

Facts for Part 2
----------------

This subsection is long and includes several details we need to discuss.  Overall, to link the random finite samples to the continuous and deterministic setup; that is, link $\widehat{T}_{\textup{ref},\epsilon,n}$ to $T_{\textup{ref},\epsilon}$, we consider the Glivenko-Cantelli class commonly used in the empirical processes analysis [@wellner2013weak].  We need some more definitions.  For the probability measure $d\mathbb{P}_X=p_XdV$ associated with the dataset, and a function $f\in C(M)$, introduce the abbreviation $$\mathbb{P}f\coloneqq\int f(x)\,d\mathbb{P}_X(x)\,.$$ Let $x_{1},\ldots,x_{n}$ be i.i.d.  sampled from $\mathbb{P}$, and denote by $$\mathbb{P}_{n}\coloneqq \frac{1}{n}\sum_{i=1}^{n}\delta_{x_{i}}$$ the corresponding empirical distribution, where $\delta_{x_{i}}$ is the Dirac delta measure supported at $x_i$. Note that we have $\mathbb{P}_{n}f = \frac{1}{n}\sum_{i=1}^{n}f(x_{i})$. Also, denote $$\widetilde{\mathbb{P}}f:=\int f(y)p_{Y}(y)\,dV(y)\,\,\,\mbox{ and }\,\,\,\widetilde{\mathbb{P}}_mf:=\frac{1}{m}\sum_{l=1}^m f(y_l)\,.$$

We now prepare some bounds for later proof. 

\[unif-bd\] Fix $\epsilon>0$. Set $\delta=\min K_{\epsilon}$. The following bounds hold for all $x\in M$: $$\begin{aligned}
\delta^2\leq K_{\textup{ref},\epsilon}(x,y)\leq{\left\lVertK\right\rVert}^{2}_{\infty},&\;\; \delta^2\leq\widehat{K}_{\textup{ref},\epsilon,n}(x,y)\leq{\left\lVertK\right\rVert}^{2}_{\infty}\,,\nonumber\\
C_1\epsilon^d\leq d_{\textup{ref},\epsilon}(x)\leq C_2\epsilon^d, &\;\; \delta^2\leq\widehat{d}_{\textup{ref},\epsilon,n}(x)\leq{\left\lVertK\right\rVert}^{2}_{\infty}\,,\label{Bound of d and dhat}\\
\frac{C_2\epsilon^d}{{\left\lVertK\right\rVert}_{\infty}}\leq M_{\textup{ref},\epsilon}(x,y)\leq \frac{{\left\lVertK\right\rVert}^{2}_{\infty}}{C_1\epsilon^{d}},&\;\;\frac{\delta^2}{{\left\lVertK\right\rVert}_{\infty}}\leq \widehat{M}_{\textup{ref},\epsilon,n}(x,y)\leq \frac{{\left\lVertK\right\rVert}^{2}_{\infty}}{\delta^2}\,,\nonumber\end{aligned}$$ where $C_1$ is a constant depending on the kernel, the curvature of the manifold and the minima of $p_X$ and $p_Y$. Similarly, $C_2$ is a constant depending on the kernel, the curvature of the manifold and the maxima of $p_X$ and $p_Y$.

The bounds for $K_{\textup{ref},\epsilon}(x,y)$ and $\hat K_{\textup{ref},\epsilon,n}(x,y)$ come from a trivial bound.  For $d_{\textup{ref},\epsilon}(x)$, we have $$\begin{aligned}
    \nonumber&\quad\int_{M}K_{\textup{ref},\epsilon}(x,y)P_{X}(y)\,dV(y)\,\geq \inf_{x'\in M}P_{X}(x')\int_{M}K_{\textup{ref},\epsilon}(x,y)\,dV(y)\\
    \nonumber
        &\geq \inf_{x',y'\in M}P_X(x')P_Y(y')\int_M \int_{M}K_{\epsilon}(x,z)K_{\epsilon}(z,y)\,dV(z)dV(y)\\
       \nonumber &= C_1\epsilon^{d}\,,\end{aligned}$$ Note that due to the randomness, $\widehat{d}_{\textup{ref},\epsilon,n}(x)$ can only be trivially bounded.  $M_{\textup{ref},\epsilon}(x,y)$ and $\hat M_{\textup{ref},\epsilon,n}(x,y)$ are bounded by combining the above bounds. 

Below, we list some functional spaces we need for the analysis, and show that they are Glivenko-Cantelli classes. 

\[Definition different classes\] Let $u\in C(M)$ and $k$ be the chosen Gaussian kernel stated in Theorem \[spec\_cong\].  Define $$\begin{aligned}
\mathcal{K}&\coloneqq\{k(x,\cdot);x\in M\}, \nonumber\\
\mathcal{K}\cdot\mathcal{K}&\coloneqq\{k(x,\cdot)k(\cdot,y);x,y\in M\}\,\nonumber\\
u\cdot\mathcal{M}&\coloneqq\{u(\cdot)M_{\textup{ref},\epsilon}(x,\cdot);x\in M\},\\
\mathcal{M}\cdot\mathcal{M}&\coloneqq\{M_{\textup{ref},\epsilon}(x,\cdot)M_{\textup{ref},\epsilon}(\cdot,y);x,y\in M\}\nonumber\\
\int\mathcal{K}\cdot\mathcal{K}&\coloneqq\Big\{\int k(x,z)k(z,\cdot)p_{Y}(z)\,dV(z);x\in M\Big\}\,.\nonumber\end{aligned}$$

\[G-C class\] The classes $\mathcal{K}$, $\mathcal{K}\cdot\mathcal{K}, f\cdot\mathcal{M}$, $\mathcal{M}\cdot\mathcal{M}$ and $\int\mathcal{K}\cdot\mathcal{K}$ are Glivenko-Cantelli classes. 

The proof of Lemma \[G-C class\] is standard, and can be found in, for example [@von2008consistency Proposition 11], so we omit the details. 

\[Proposition: Tepsn conv to Teps compactly\] For a fixed $\epsilon>0$, $\widehat{T}_{\textup{ref},\epsilon,n}$ converges to $T_{\textup{ref},\epsilon}$ collectively compactly $\textup{a.s.}$ as $n\rightarrow\infty$.

We verify the collectively compact convergence.  For (I), let $\widehat{M}_{\textup{ref},\epsilon,n}^{(d)}(x,y)\coloneqq\frac{K_{\textup{ref},\epsilon}(x,y)}{\widehat{d}_{\textup{ref},\epsilon,n}^{(d)}(x)}\in C(M\times M)$. Pick any $f\in C(M)$. By the triangle inequality, we have $$\begin{aligned}
    \nonumber {\left\lVert\widehat{T}_{\textup{ref},\epsilon,n}f - T_{\textup{ref},\epsilon}f\right\rVert}_{\infty} =&\, {\left\lVert\frac{1}{n}\sum_{i=1}^{n}\widehat{M}_{\textup{ref},\epsilon,n}(x,x_{i})f(x_{i}) - \int_{M}M_{\textup{ref},\epsilon}(x,y)f(y)p_{X}(y)\,dV(y)\right\rVert}_{\infty}
    \\\nonumber \leq&\, \sup_{x}{\left |\mathbb{P}_{n}\widehat{M}_{\textup{ref},\epsilon,n}(x,\cdot)f(\cdot) - \mathbb{P}M_{\textup{ref},\epsilon}(x,\cdot)f(\cdot)\right |}
    \\ \label{term1} \leq&\, \sup_{x}{\left |\mathbb{P}_{n}\widehat{M}_{\textup{ref},\epsilon,n}(x,\cdot)f(\cdot) - \mathbb{P}_{n}\widehat{M}^{(d)}_{\textup{ref},\epsilon,n}(x,\cdot)f(\cdot)\right |}
    \\ \label{term2} &\quad+\, \sup_{x}{\left |\mathbb{P}_{n}\widehat{M}^{(d)}_{\textup{ref},\epsilon,n}(x,\cdot)f(\cdot) - \mathbb{P}_{n}M_{\textup{ref},\epsilon}(x,\cdot)f(\cdot)\right |}
    \\ \label{term3} &\quad+\, \sup_{x}{\left |\mathbb{P}_{n}M_{\textup{ref},\epsilon}(x,\cdot)f(\cdot) - \mathbb{P}M_{\textup{ref},\epsilon}(x,\cdot)f(\cdot)\right |}\,.\end{aligned}$$ We bound the three terms (\[term1\]) (\[term2\]) and (\[term3\]) respectively by Lemma \[unif-bd\].  By a direct expansion, $$\begin{aligned}
    \nonumber (\ref{term1}) \,&= \sup_{x}{\left |\frac{1}{n}\sum_{i=1}^{n}\widehat{M}_{\textup{ref},\epsilon,n}(x,x_{i})f(x_{i}) - \frac{1}{n}\sum_{i=1}^{n}\widehat{M}^{(d)}_{\textup{ref},\epsilon,n}(x,x_{i})f(x_{i})\right |}
    \\\nonumber &=
    \sup_{x}{\left |\frac{1}{n}\sum_{i=1}^{n}\left[\widehat{M}_{\textup{ref},\epsilon,n}(x,x_{i})-\widehat{M}^{(d)}_{\textup{ref},\epsilon,n}(x,x_{i})\right]f(x_{i})\right |}\,,
    \end{aligned}$$ which is bounded by $$\begin{aligned}
 \nonumber &\quad
    \sup_{x}\frac{1}{n}\sum_{i=1}^{n}{\left |\widehat{M}_{\textup{ref},\epsilon,n}(x,x_{i})-\widehat{M}^{(d)}_{\textup{ref},\epsilon,n}(x,x_{i})\right |}{\left |f(x_{i})\right |}
            \\\nonumber &\leq
    {\left\lVertf\right\rVert}_{\infty}\sup_{x,y}{\left |\frac{\widehat{K}_{\textup{ref},\epsilon,n}(x,y)-K_{\textup{ref},\epsilon}(x,y)}{\widehat{d}_{\textup{ref},\epsilon,n}(x)}\right |}
    \\\nonumber &\leq \frac{{\left\lVertf\right\rVert}_{\infty}}{\delta^{2}}\sup_{x,y}{\left |\widehat{K}_{\textup{ref},\epsilon,n}(x,y)-K_{\textup{ref},\epsilon}(x,y)\right |}\,,
    \end{aligned}$$ where the last bound comes from .  Then, we spell out $K_{\textup{ref},\epsilon,n}(x,y)$ and $K_{\textup{ref},\epsilon}(x,y)$ so that (\[term1\]) is further bounded by $$\begin{aligned}
&\nonumber
    \frac{{\left\lVertf\right\rVert}_{\infty}}{\delta^{2}}\sup_{x,y}{\left |\frac{1}{m}\sum_{j}^{m}K_{\epsilon}(x,z_{j})K_{\epsilon}(z_{j},y)-\int_{M}K_{\epsilon}(x,z)K_{\epsilon}(z,y)p_{Y}(z)\,dV(z)\right |}
    \\\nonumber =&\,
    \frac{{\left\lVertf\right\rVert}_{\infty}}{\delta^{2}}\sup_{x,y}{\left |\widetilde{\mathbb{P}}_{m}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,y) - \widetilde{\mathbb{P}}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,y)\right |}
    \\\nonumber \leq &\,
    \frac{{\left\lVertf\right\rVert}_{\infty}}{\delta^{2}}\sup_{F\in\mathcal{K}\cdot\mathcal{K}}{\left |\widetilde{\mathbb{P}}_{m}F - \widetilde{\mathbb{P}}F\right |} \,,\end{aligned}$$ which goes to $0$ a.s.  as $n\rightarrow\infty$ since $\mathcal{K}\cdot\mathcal{K}$ is Glivenko-Cantelli class by Lemma \[G-C class\] and $m=n^{\beta}$ for some $\beta\in(0,1)$. For term (\[term2\]), we again expand it directly: $$\begin{aligned}
    \nonumber (\ref{term2}) =& \sup_{x}{\left |\frac{1}{n}\sum_{i=1}^{n}\widehat{M}^{(d)}_{\textup{ref},\epsilon,n}(x,x_{i})f(x_{i}) - \frac{1}{n}\sum_{i=1}^{n}M_{\textup{ref},\epsilon}(x,x_{i})f(x_{i})\right |}
    \\\nonumber \leq&\,
    \sup_{x}\frac{1}{n}\sum_{i=1}^{n}{\left |\widehat{M}^{(d)}_{\textup{ref},\epsilon,n}(x,x_{i})-M_{\textup{ref},\epsilon}(x,x_{i})\right |}{\left |f(x_{i})\right |}
            \\\nonumber \leq&\,
    {\left\lVertf\right\rVert}_{\infty}\sup_{x,y}{\left |\frac{K_{\textup{ref},\epsilon}(x,y)}{\widehat{d}_{\textup{ref},\epsilon,n}(x)}-\frac{K_{\textup{ref},\epsilon}(x,y)}{d_{\textup{ref},\epsilon}(x)}\right |} 
    \\\nonumber \leq&\,
    \frac{{\left\lVertf\right\rVert}_{\infty}{\left\lVertK\right\rVert}_{\infty}^{2}}{\delta^{4}}\sup_{x}{\left |\widehat{d}_{\textup{ref},\epsilon,n}(x) - d_{\textup{ref},\epsilon}(x)\right |}\,,
    \end{aligned}$$ where in the last bound we apply the fact that if $A,B\geq C>0$, then ${\left |A^{\beta}-B^{\beta}\right |}\leq\frac{1}{C^{1-\beta}}{\left |A-B\right |}$. Then, by spelling out $\widehat{d}_{\textup{ref},\epsilon,n}(x)$ and $d_{\textup{ref},\epsilon}(x)$, (\[term2\]) is further bounded by $$\begin{aligned}
\nonumber &\,
    \frac{{\left\lVertf\right\rVert}_{\infty}{\left\lVertK\right\rVert}_{\infty}^{2}}{\delta^{4}}\sup_{x}{\left |\mathbb{P}_{n}\widetilde{\mathbb{P}}_{m}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,\star)-\mathbb{P}\widetilde{\mathbb{P}}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,\star)\right |}
    \\ \label{term4}\leq&\,
    \frac{{\left\lVertf\right\rVert}_{\infty}{\left\lVertK\right\rVert}_{\infty}^{2}}{\delta^{4}}\sup_{x}{\left |\mathbb{P}_{n}\widetilde{\mathbb{P}}_{m}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,\star)-\mathbb{P}_{n}\widetilde{\mathbb{P}}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,\star)\right |}
    \\ \label{term5}&\quad+\,
    \frac{{\left\lVertf\right\rVert}_{\infty}{\left\lVertK\right\rVert}_{\infty}^{2}}{\delta^{4}}\sup_{x}{\left |\mathbb{P}_{n}\widetilde{\mathbb{P}}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,\star)-\mathbb{P}\widetilde{\mathbb{P}}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,\star)\right |}\end{aligned}$$ Clearly, the term $(\ref{term5})\rightarrow 0\quad\textup{a.s.}$ by lemma \[G-C class\].  For term (\[term4\]): $$\begin{aligned}
    \nonumber (\ref{term4}) =&\, \sup_{x}{\left |\mathbb{P}_{n}\widetilde{\mathbb{P}}_{m}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,\star)-\mathbb{P}_{n}\widetilde{\mathbb{P}}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,\star)\right |}
    \\\nonumber =&\,
    \sup_{x}{\left |\frac{1}{n}\sum_{i=1}^{n}\frac{1}{m}\sum_{j=1}^{m}K_{\epsilon}(x,z_{j})K_{\epsilon}(z_{j},x_{i})-\frac{1}{n}\sum_{i=1}^{n}\int K_{\epsilon}(x,z)K_{\epsilon}(z,x_{i})p_{Y}(z)\,dV(z)\right |}
    \\\nonumber \leq&\,
    \sup_{x}\frac{1}{n}\sum_{i=1}^{n}{\left |\frac{1}{m}\sum_{j=1}^{m}K_{\epsilon}(x,z_{j})K_{\epsilon}(z_{j},x_{i})-\int K_{\epsilon}(x,z)K_{\epsilon}(z,x_{i})p_{Y}(z)\,dV(z)\right |}
    \\\nonumber \leq&
    \sup_{x,y}{\left |\frac{1}{m}\sum_{j=1}^{m}K_{\epsilon}(x,z_{j})K_{\epsilon}(z_{j},y)-\int K_{\epsilon}(x,z)K_{\epsilon}(z,y)p_{Y}(z)\,dV(z)\right |}
    \\\nonumber =&\,
    \sup_{x,y}{\left |\widetilde{\mathbb{P}}_{m}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,y) - \widetilde{\mathbb{P}}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,y)\right |}\,,\end{aligned}$$ which tends to $0$ a.s.  as $n\rightarrow\infty$. Hence, the term $(\ref{term2})\rightarrow 0\;\; \textup{a.s.  as}\;n\rightarrow\infty$. Finally, the term $(\ref{term3})\rightarrow 0\;\; \textup{a.s.  as}\;n\rightarrow\infty$ by Lemma \[G-C class\].  So the condition (I) is verified. 

Next, we verify (II).  Since $T_{\textup{ref},\epsilon}$ is compact, it is enough to show the set $\cup_{n}\widehat{T}_{\textup{ref},\epsilon,n}(B)$ is relatively compact, where $B\subset (C(M),\,\|\cdot\|_\infty)$ is the unit ball centered at $0$. By the Arzela-Ascoli theorem, we need to show that the set $\cup_{n}\widehat{T}_{\textup{ref},\epsilon,n}(B)$ is:

(i) pointwisely bounded, and

(ii) equicontinuous. 

For (i), for any $x\in M$, by Lemma \[unif-bd\]: $$\begin{aligned}
    \nonumber&\sup_{f\in B,\, n\in \mathbb{N}}{\left\lVert\widehat{T}_{\textup{ref},\epsilon,n}f(x)\right\rVert}_{\infty}
    =
    \sup_{f\in B,\, n\in \mathbb{N}}{\left |\frac{1}{n}\sum_{i=1}^{n}\widehat{M}_{\textup{ref},\epsilon,n}(x,x_{i})f(x_{i})\right |}
            \\\nonumber\leq&\,
    \sup_{f\in B,\, n\in \mathbb{N}}{\left\lVertf\right\rVert}_{\infty}\frac{1}{n}\sum_{i=1}^{n}{\left |\widehat{M}_{\textup{ref},\epsilon,n}(x,x_{i})\right |}
   \leq \frac{{\left\lVertK\right\rVert}_{\infty}^{2}}{\delta^{2}} < \infty.\end{aligned}$$ For (ii), for any $x,y\in M$ that are close, a direct expansion leads to $$\begin{aligned}
    \nonumber &\sup_{f\in B,\, n\in \mathbb{N}}{\left |\widehat{T}_{\textup{ref},\epsilon,n}f(y) - \widehat T_{\textup{ref},\epsilon,n}f(x)\right |} 
    \\\nonumber=& \sup_{f\in B,\, n\in \mathbb{N}}{\left |\frac{1}{n}\sum_{i=1}^{n}\widehat{M}_{\textup{ref},\epsilon,n}(y,x_{i})f(x_{i}) - \frac{1}{n}\sum_{i=1}^{n}\widehat{M}_{\textup{ref},\epsilon,n}(x,x_{i})f(x_{i})\right |}
    \\\nonumber \leq&\,
    \sup_{f\in B,\, n\in \mathbb{N}}\frac{1}{n}\sum_{i=1}^{n}{\left |\widehat{M}_{\textup{ref},\epsilon,n}(y,x_{i})-\widehat{M}_{\textup{ref},\epsilon,n}(x,x_{i})\right |}{\left |f(x_{i})\right |}
    \\\nonumber \leq&\,
    \sup_{f\in B}{\left\lVertf\right\rVert}_{\infty}\sup_{z}{\left |\widehat{M}_{\textup{ref},\epsilon,n}(y,z)-\widehat{M}_{\textup{ref},\epsilon,n}(x,z)\right |}\,. 
    \end{aligned}$$ Clearly, since $f\in B$, ${\left\lVertf\right\rVert}_{\infty}\leq 1$. So $\sup_{f\in B,\, n\in \mathbb{N}}{\left |\widehat{T}_{\textup{ref},\epsilon,n}f(y) - \widehat T_{\textup{ref},\epsilon,n}f(x)\right |}$ is further bounded by $$\begin{aligned}
    \nonumber &\,
    \sup_{z}{\left |\frac{\widehat{K}_{\textup{ref},\epsilon,n}(y,z)}{\widehat{d}_{\textup{ref},\epsilon,n}(y)}-\frac{\widehat{K}_{\textup{ref},\epsilon,n}(x,z)}{\widehat{d}_{\textup{ref},\epsilon,n}(x)}\right |}
            \\\nonumber \leq &\,
    \frac{1}{\delta^{4}}\sup_{z}{\left |\widehat{d}_{\textup{ref},\epsilon,n}(y)\widehat{K}_{\textup{ref},\epsilon,n}(x,z)-\widehat{d}_{\textup{ref},\epsilon,n}(x)\widehat{K}_{\textup{ref},\epsilon,n}(y,z)\right |}\,,
    \end{aligned}$$ where the last inequality comes again from Lemma \[unif-bd\].  The right hand side can further be bounded by $$\begin{aligned}
    \nonumber &\,
    \frac{1}{\delta^{4}}\sup_{z}\big|\widehat{d}_{\textup{ref},\epsilon,n}(y)[\widehat{K}_{\textup{ref},\epsilon,n}(x,z)-\widehat{K}_{\textup{ref},\epsilon,n}(y,z)]\\
    &\nonumber\qquad+[\widehat{d}_{\textup{ref},\epsilon,n}(y)-\widehat{d}_{\textup{ref},\epsilon,n}(x)]\widehat{K}_{\textup{ref},\epsilon,n}(y,z)\big|
    \\\label{Bound for equicontinuity part 1} \leq&\,
    \frac{{\left\lVertK\right\rVert}_{\infty}^{2}}{\delta^{4}}\sup_{z}{\left |\widehat{K}_{\textup{ref},\epsilon,n}(x,z)-\widehat{K}_{\textup{ref},\epsilon,n}(y,z)\right |}
    \\\label{Bound for equicontinuity part 2} &\qquad+
    \frac{{\left\lVertK\right\rVert}_{\infty}^{2}}{\delta^{4}}{\left |\widehat{d}_{\textup{ref},\epsilon,n}(y)-\widehat{d}_{\textup{ref},\epsilon,n}(x)\right |}\,.\end{aligned}$$ To finish the equicontinuity argument, we bound and .  By spelling out $\widehat{K}_{\textup{ref},\epsilon,n}(x,z)$ and $\widehat{K}_{\textup{ref},\epsilon,n}(y,z)$; that is, $$\begin{aligned}
    \nonumber &\sup_{z}{\left |\widehat{K}_{\textup{ref},\epsilon,n}(x,z)-\widehat{K}_{\textup{ref},\epsilon,n}(y,z)\right |}\\
    =&\,
    \sup_{z}{\left |\frac{1}{m}\sum_{j}^{m}K_{\epsilon}(x,z_{j})K_{\epsilon}(z_{j},z)-\frac{1}{m}\sum_{j}^{m}K_{\epsilon}(y,z_{j})K_{\epsilon}(z_{j},z)\right |}\,,\nonumber\end{aligned}$$ is bounded by: $$\begin{aligned}
    \nonumber &
    \sup_{z}\frac{1}{m}\sum_{j}^{m}{\left |K_{\epsilon}(x,z_{j})-K_{\epsilon}(y,z_{j})\right |}{\left |K_{\epsilon}(z_{j},z)\right |}
    \\\nonumber \leq&\,
    {\left\lVertK_{\epsilon}\right\rVert}_{\infty}\frac{1}{m}\sum_{j}^{m}{\left |K_{\epsilon}(x,z_{j})-K_{\epsilon}(y,z_{j})\right |}
     \leq
    {\left\lVertK_{\epsilon}\right\rVert}_{\infty}\sup_{z}{\left |K_{\epsilon}(x,z)-K_{\epsilon}(y,z)\right |}\,,\end{aligned}$$ which is controlled by $Cd(x,y)$ for some constant $C>0$ due to the continuity of $K_{\epsilon}$ and the compactness of $M$. Similarly, by spelling out $\widehat{d}_{\textup{ref},\epsilon,n}(y)-\widehat{d}_{\textup{ref},\epsilon,n}(x)$; that is, $${\widehat{d}_{\textup{ref},\epsilon,n}(y)-\widehat{d}_{\textup{ref},\epsilon,n}(x)}
= {\frac{1}{n}\sum_{i=1}^{n}\widehat{K}_{\textup{ref},\epsilon,n}(y,x_{i})-\frac{1}{n}\sum_{i=1}^{n}\widehat{K}_{\textup{ref},\epsilon,n}(x,x_{i})}\,,$$ is bounded by: $$\begin{aligned}
   \nonumber &
    {\left |\frac{1}{n}\sum_{i=1}^{n}\frac{1}{m}\sum_{j}^{m}K_{\epsilon}(y,z_{j})K_{\epsilon}(z_{j},x_{i})-\frac{1}{n}\sum_{i=1}^{n}\frac{1}{m}\sum_{j}^{m}K_{\epsilon}(x,z_{j})K_{\epsilon}(z_{j},x_{i})\right |}
    \\\nonumber \leq&
    \frac{1}{nm}\sum_{i,j}{\left |K_{\epsilon}(y,z_{j})-K_{\epsilon}(x,z_{j})\right |}{\left |K_{\epsilon}(z_{j},x_{i})\right |}
    \\\nonumber \leq&
    {\left\lVertK_{\epsilon}\right\rVert}_{\infty}\sup_{z}{\left |K_{\epsilon}(x,z)-K_{\epsilon}(y,z)\right |}\,,\end{aligned}$$ which is again controlled by $Cd(x,y)$. So (ii) is verified.  Then condition (II) is verified by Arzela-Ascoli theorem.  We thus finish the proof of collectively compact convergence. 

Next, we control various functional classes that concern us with the bound shown in Lemma \[cover\_num\_gau\].  These controls are needed when we derive the convergence rate.  While the proof is standard, we provide details of how the landmark set plays a role in the bound. 

\[cover\_num\_bd\] Take $u\in C(M)$ and $\epsilon>0$. Take $k(x,y):=K_\epsilon(x,y)$ in Definition \[Definition different classes\].  For $r>0$, we have the following bounds: $$N(\mathcal{K}\cdot\mathcal{K},r,{\left\lVert\cdot\right\rVert}_{\infty})\leq N\Big(\mathcal{K},\frac{r}{2{\left\lVertK_{\epsilon}\right\rVert}},{\left\lVert\cdot\right\rVert}_{\infty}\Big)$$ $$N(\int\mathcal{K}\cdot\mathcal{K},r,{\left\lVert\cdot\right\rVert}_{\infty})\leq N\Big(\mathcal{K},\frac{r}{{\left\lVertK_{\epsilon}\right\rVert}},{\left\lVert\cdot\right\rVert}_{\infty}\Big)$$ $$N(u\cdot\mathcal{M},r,{\left\lVert\cdot\right\rVert}_{\infty})\leq N\Big(\mathcal{K},\frac{rC^{2}\epsilon^{2d}}{2{\left\lVertu\right\rVert}{\left\lVertK_{\epsilon}\right\rVert}^{3}},{\left\lVert\cdot\right\rVert}_{\infty}\Big)$$ $$N(\mathcal{M}\cdot\mathcal{M},r,{\left\lVert\cdot\right\rVert}_{\infty})\leq N\Big(\mathcal{K},\frac{rC^{3}\epsilon^{3d}}{4{\left\lVertK_{\epsilon}\right\rVert}^{5}},{\left\lVert\cdot\right\rVert}_{\infty}\Big)\,.$$

For the class $\mathcal{K}\cdot\mathcal{K}$, pick any $x_{1},x_{2},y_{1},y_{2}\in M$: $$\begin{aligned}
    \nonumber&{\left |K_{\epsilon}(x_{1},z)K_{\epsilon}(y_{1},z)-K_{\epsilon}(x_{2},z)K_{\epsilon}(y_{2},z)\right |}
    \\\nonumber=&\,
    \big|K_{\epsilon}(x_{1},z)K_{\epsilon}(y_{1},z)-K_{\epsilon}(x_{1},z)K_{\epsilon}(y_{2},z)\\
    \nonumber&\qquad\qquad+K_{\epsilon}(x_{1},z)K_{\epsilon}(y_{2},z)-K_{\epsilon}(x_{2},z)K_{\epsilon}(y_{2},z)\big|
    \\\nonumber\leq&\,
    {\left\lVertK_{\epsilon}\right\rVert}_\infty\big({\left |K_{\epsilon}(y_{1},z)-K_{\epsilon}(y_{2},z)\right |} + {\left |K_{\epsilon}(x_{1},z)-K_{\epsilon}(x_{2},z)\right |}\big)\,.\end{aligned}$$ This implies a $\frac{r}{2{\left\lVertK_{\epsilon}\right\rVert}}$-cover of $\mathcal{K}$ induces a $r$-cover of $\mathcal{K}\cdot\mathcal{K}$. Hence, $N(\mathcal{K}\cdot\mathcal{K},r,{\left\lVert\cdot\right\rVert}_{\infty})\leq N(\mathcal{K},\frac{r}{2{\left\lVertK_{\epsilon}\right\rVert}},{\left\lVert\cdot\right\rVert}_{\infty})$.

For the class $\int\mathcal{K}\cdot\mathcal{K}$, pick any $x_{1},x_{2}\in M$: $$\begin{aligned}
    \nonumber&{\left |\int K_{\epsilon}(x_{1},z)K_{\epsilon}(z,y)\,d\widetilde{\mathbb{P}}(z)-\int K_{\epsilon}(x_{2},z)K_{\epsilon}(z,y)\,d\widetilde{\mathbb{P}}(z)\right |}
   \\\nonumber \leq&\,
    {\left\lVertK_{\epsilon}\right\rVert}{\left |K_{\epsilon}(x_{1},z)-K_{\epsilon}(x_{2},z)\right |}\,,\end{aligned}$$ which implies that a $\frac{r}{{\left\lVertK_{\epsilon}\right\rVert}}$-cover of $\mathcal{K}$ induces a $r$-cover of $\int\mathcal{K}\cdot\mathcal{K}$. Hence, $N(\int\mathcal{K}\cdot\mathcal{K},r,{\left\lVert\cdot\right\rVert}_{\infty})\leq N(\mathcal{K},\frac{r}{{\left\lVertK_{\epsilon}\right\rVert}},{\left\lVert\cdot\right\rVert}_{\infty})$.

For the class $u\cdot\mathcal{M}$, pick any $x_{1},x_{2}\in M$: $$\begin{aligned}
    \nonumber&\quad{\left |u(y)M_{\textup{ref},\epsilon}(x_{1},y)-u(y)M_{\textup{ref},\epsilon}(x_{2},y)\right |}
    \\\nonumber&\leq
    {\left\lVertu\right\rVert}_\infty {\left |\frac{K_{\textup{ref},\epsilon}(x_{1},y)}{d_{\textup{ref},\epsilon}(x_{1})}-\frac{K_{\textup{ref},\epsilon}(x_{2},y)}{d_{\textup{ref},\epsilon}(x_{2})}\right |}
    \\\nonumber&\leq
    \frac{{\left\lVertu\right\rVert}_\infty}{C^{2}\epsilon^{2d}}\left({\left\lVertd_{\textup{ref},\epsilon}\right\rVert}{\left |K_{\textup{ref},\epsilon}(x_{1},y)-K_{\textup{ref},\epsilon}(x_{2},y)\right |}+{\left\lVertK_{\textup{ref},\epsilon}\right\rVert}{\left\lVertd_{\textup{ref},\epsilon}(x_{1})-d_{\textup{ref},\epsilon}(x_{2})\right\rVert}\right)
    \\\nonumber&\leq
    \frac{2{\left\lVertu\right\rVert}_\infty{\left\lVertK_{\epsilon}\right\rVert}_\infty^{2}}{C^{2}\epsilon^{2d}}{\left |K_{\textup{ref},\epsilon}(x_{1},y)-K_{\textup{ref},\epsilon}(x_{2},y)\right |}
    \\\nonumber&=
    \frac{2{\left\lVertu\right\rVert}_\infty{\left\lVertK_{\epsilon}\right\rVert}_\infty^{2}}{C^{2}\epsilon^{2d}}{\left |\int K_{\epsilon}(x_{1},z)K_{\epsilon}(z,y)\,d\widetilde{\mathbb{P}}(z)-\int K_{\epsilon}(x_{2},z)K_{\epsilon}(z,y)\,d\widetilde{\mathbb{P}}(z)\right |}
    \\\nonumber&\leq
    \frac{2{\left\lVertu\right\rVert}_\infty{\left\lVertK_{\epsilon}\right\rVert}_\infty^{3}}{C^{2}\epsilon^{2d}}{\left |K_{\epsilon}(x_{1},z)-K_{\epsilon}(x_{2},z)\right |}\,,\end{aligned}$$ which implies that a $\frac{rC^{2}\epsilon^{2d}}{2{\left\lVertu\right\rVert}_\infty{\left\lVertK_{\epsilon}\right\rVert}_\infty^{3}}$-cover of $\mathcal{K}$ induces a $r$-cover of $u\cdot\mathcal{M}$, hence $N(u\cdot\mathcal{M},r,{\left\lVert\cdot\right\rVert}_{\infty})\leq N\big(\mathcal{K},\frac{rC^{2}\epsilon^{2d}}{2{\left\lVertu\right\rVert}_\infty{\left\lVertK_{\epsilon}\right\rVert}_\infty^{3}},{\left\lVert\cdot\right\rVert}_{\infty}\big)$.

For the class $\mathcal{M}\cdot\mathcal{M}$, pick any $x_{1},x_{2},y_{1},y_{2}\in M$ and carry out a similar bound: $$\begin{aligned}
    \nonumber&\quad{\left |M_{\epsilon}(x_{1},z)M_{\epsilon}(y_{1},z)-M_{\epsilon}(x_{2},z)M_{\epsilon}(y_{2},z)\right |}
            \\\nonumber&\leq
    {\left\lVertM_{\epsilon}\right\rVert}_\infty\big({\left |M_{\epsilon}(y_{1},z)-M_{\epsilon}(y_{2},z)\right |} + {\left |M_{\epsilon}(x_{1},z)-M_{\epsilon}(x_{2},z)\right |}\big)
    \\\nonumber&\leq
    \frac{4{\left\lVertK_{\epsilon}\right\rVert}^{5}_\infty}{C^{3}\epsilon^{3d}}{\left |K_{\epsilon}(x_{1},z)-K_{\epsilon}(x_{2},z)\right |}\,,\end{aligned}$$ which implies that a $\frac{rC^{3}\epsilon^{3d}}{4{\left\lVertK_{\epsilon}\right\rVert}^{5}}$-cover of $\mathcal{K}$ induces a $r$-cover of $\mathcal{M}\cdot\mathcal{M}$. Hence, $N(\mathcal{M}\cdot\mathcal{M},r,{\left\lVert\cdot\right\rVert}_{\infty})\leq N\big(\mathcal{K},\frac{rC^{3}\epsilon^{3d}}{4{\left\lVertK_{\epsilon}\right\rVert}^{5}},{\left\lVert\cdot\right\rVert}_{\infty}\big)$. The constant C above is the C in Lemma \[unif-bd\], depends on the kernel $K$, the curvature of the manifold and the density functions $p_{X}, p_{Y}$.

With Lemma \[cover\_num\_bd\], we can derive the following technical lemma.  We provide a detailed proof to show how the landmark set impacts the bound. 

\[cant\_class\_rate\] Take $\epsilon>0$, $u\in C(M)$, and $k(x,y):=K_\epsilon(x,y)=e^{-{\left\lVertx-y\right\rVert}^2/\epsilon}$ in Definition \[Definition different classes\].  Let $\mathcal{F}_{\epsilon}\coloneqq(\mathcal{K}\cdot\mathcal{K})\cup(\int\mathcal{K}\cdot\mathcal{K})\cup(u\cdot\mathcal{M})\cup(\mathcal{M}\cdot\mathcal{M})$. With probability $1-\mathcal O(n^{-2})$, we have: $$\label{GC class bound for all classes}
\sup_{f\in\mathcal{F}_{\epsilon}}{\left |\mathbb{P}_{n}f-\mathbb{P}f\right |}=\mathcal{O}\left(\frac{\sqrt{-\log \epsilon}+\sqrt{\log n}}{\sqrt{n}}\right)\,,$$ where the implied constant depends on $d$, $\|u\|_\infty$ and the constants shown in the entropy bound shown in Theorem \[entropy\_bd\]. 

Note that the probability event space that holds depends on the chosen $u$, and the implied constant depends on $\|u\|_\infty$. This is critical when we carry out the final spectral convergence proof. 

By plugging $\delta=1/n^{2}$ into Theorem \[entropy\_bd\], we have $$\begin{aligned}
    \nonumber\sup_{F\in\mathcal{K}\cdot\mathcal{K}}{\left |\mathbb{P}_{n}F-\mathbb{P}F\right |}\,&\leq\frac{C_{E}}{\sqrt{n}}\int_{0}^{\infty}\sqrt{\log N(\mathcal{K}\cdot\mathcal{K},r,L_{2}(\mathbb{P}_{n}))}\,dr + \sqrt{\frac{1}{n}\log(2n^{2})}
    \\\nonumber&=
    \frac{C_{E}}{\sqrt{n}}\int_{0}^{1}\sqrt{\log N(\mathcal{K}\cdot\mathcal{K},r,L_{2}(\mathbb{P}_{n}))}\,dr + c_{1}\sqrt{\frac{\log n}{n}}\,,
    \end{aligned}$$ where we use the fact that $N(\mathcal{K}\cdot\mathcal{K},r,L_{2}(\mathbb{P}_{n}))=1$ when $r>1$. Indeed, for any $x,y\in M$ and $n\in\mathbb{N}$, we have $$\begin{aligned}
    \|K_\epsilon(x,\cdot)K_\epsilon(\cdot,y)\|^2_{L_{2}(\mathbb{P}_{n})}
    \leq \frac{1}{n}\sum_{l=1}^n|K_\epsilon(x,z_l)K_\epsilon(z_l,y)|^2\leq \|K_\epsilon\|_{\infty}^2=1\,.\nonumber
    \end{aligned}$$ Similarly, note that for any $x,x',y,y'\in M$ and $n\in \mathbb{N}$, we have $$\|K_\epsilon(x,\cdot)K_\epsilon(\cdot,y)-K_\epsilon(x',\cdot)K_\epsilon(\cdot,y')\|^2_{L_{2}(\mathbb{P}_{n})}\leq \|K_\epsilon(x,\cdot)K_\epsilon(\cdot,y)-K_\epsilon(x',\cdot)K_\epsilon(\cdot,y')\|^2_{\infty}\nonumber\,,$$ so we immediately have $$N(\mathcal{K}\cdot\mathcal{K},r,L_{2}(\mathbb{P}_{n}))\leq N(\mathcal{K}\cdot\mathcal{K},r,{\left\lVert\cdot\right\rVert}_{\infty})\nonumber\,$$ and hence $\int_{0}^{1}\sqrt{\log N(\mathcal{K}\cdot\mathcal{K},r,L_{2}(\mathbb{P}_{n}))}\,dr\leq
    \int_{0}^{1}\sqrt{\log N(\mathcal{K}\cdot\mathcal{K},r,{\left\lVert\cdot\right\rVert}_{\infty})}\,dr$. By Proposition \[cover\_num\_bd\], we have the bound that $$\begin{aligned}
    \nonumber &\quad   \int_{0}^{1}\sqrt{\log N(\mathcal{K}\cdot\mathcal{K},r,{\left\lVert\cdot\right\rVert}_{\infty})}\,dr 
    \leq
    \int_{0}^{1}\sqrt{\log N\Big(\mathcal{K},\frac{r}{2{\left\lVertK_{\epsilon}\right\rVert}_\infty},{\left\lVert\cdot\right\rVert}_{\infty}\Big)}\,dr 
    \\\nonumber&\leq \sqrt{2d}\int_0^1\sqrt{[48\sqrt{2d} D_M-\log \epsilon]-\log r}dr
    \\\nonumber &=\sqrt{2d}\Big(\frac{\sqrt{\pi}}{2}e^{c_\epsilon}\textup{erfc}(c_\epsilon)+\sqrt{c_\epsilon}\Big)\,,
   \end{aligned}$$ where $\textup{erfc}$ is the complementary error function, $c_\epsilon:= 48\sqrt{2d} D_M-\log \epsilon$, the second inequality comes from Theorem \[cover\_num\_gau\] and the fact that $\|K_\epsilon\|_\infty=1$. Note that since $c_\epsilon\approx -\log \epsilon$ when $\epsilon$ is small and $\textup{erfc}(x)\approx \frac{e^{-x^2}}{\sqrt{\pi}x}$, we know that $\frac{\sqrt{\pi}}{2}e^{c_\epsilon}\textup{erfc}(c_\epsilon)\to 0$ when $\epsilon$ tends to $0$. As a result, $$\begin{aligned}
    \nonumber \int_{0}^{1}\sqrt{\log N(\mathcal{K}\cdot\mathcal{K},r,{\left\lVert\cdot\right\rVert}_{\infty})}\,dr&\leq
    2\sqrt{2d}\sqrt{-\log \epsilon}
    \end{aligned}$$ when $\epsilon$ is sufficiently small.  By combining the above bounds, we have $$\begin{aligned}
    \nonumber \sup_{F\in\mathcal{K}\cdot\mathcal{K}}{\left |\mathbb{P}_{n}F-\mathbb{P}F\right |}\,&\leq  \frac{2\sqrt{2d}c\sqrt{-\log \epsilon}}{\sqrt{n}} + \frac{c_1\sqrt{\log n}}{\sqrt{n}}=\mathcal{O}\Big(\frac{\sqrt{-\log \epsilon}+\sqrt{\log n}}{\sqrt{n}}\Big)\,.\end{aligned}$$ By a similar argument we have the bound for $\sup_{F\in\int\mathcal{K}\cdot\mathcal{K}}{\left |\mathbb{P}_{n}F-\mathbb{P}F\right |}$ , $\sup_{F\in u\cdot\mathcal{M}}{\left |\mathbb{P}_{n}F-\mathbb{P}F\right |}$ and $\sup_{F\in \mathcal{M}\cdot\mathcal{M}}{\left |\mathbb{P}_{n}F-\mathbb{P}F\right |}$. Note that the implied constant of the bound for $\sup_{F\in u\cdot\mathcal{M}}{\left |\mathbb{P}_{n}F-\mathbb{P}F\right |}$ depends on $\|u\|_\infty$. The result follows. 

Below, we prepare several technical lemmas to control the spectral convergence rate. 

\[d\_hat\_bd\] Suppose $m=n^\beta$ for $\beta\in (0,1)$ and take $\epsilon=\epsilon(n)$ so that $\epsilon\to 0$ when $n\to \infty$. If we further assume that $\frac{\sqrt{\log\frac{1}{\epsilon}}+\sqrt{\log m}}{\sqrt{m}\epsilon^{d}}\rightarrow 0$, we have with probability $1-\mathcal{O}(m^{-2})$: $$\begin{aligned}
&\;\; C_1 \epsilon^{d}/2\leq\widehat{d}_{\textup{ref},\epsilon,n}(x)\leq 2C_2 \epsilon^{d} \label{Bound of d and dhat prob}\,,\end{aligned}$$ and with probability $1-\mathcal{O}(n^{-2})$: $$\begin{aligned}
&\;\; C_1 \epsilon^{d}/2\leq d_{\textup{ref},\epsilon,n}(x)\leq 2C_2 \epsilon^{d} \nonumber\,,
$$ where $C_1,C_2>0$ are constants defined in Lemma \[unif-bd\]. 

Note the difference between this lemma and Lemma \[unif-bd\].  This lemma says that while usually the kernels involved in the analysis have a wide range, with high probability, the range is well controlled. 

By the same calculation (e.g. Step 2 in the proof of Theorem \[spec\_cong\]) we have $$\begin{aligned}
    \nonumber
    \nonumber&
    \sup_{x}{\left |\widehat{d}_{\textup{ref},\epsilon,n}(x) - d_{\textup{ref},\epsilon}(x)\right |}
    \\\nonumber =&\,
    \sup_{x}{\left |\mathbb{P}_{n}\widetilde{\mathbb{P}}_{m}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,\star)-\mathbb{P}\widetilde{\mathbb{P}}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,\star)\right |}
    \\\nonumber \leq&\,
    \sup_{x}{\left |\mathbb{P}_{n}\widetilde{\mathbb{P}}_{m}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,\star)-\mathbb{P}_{n}\widetilde{\mathbb{P}}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,\star)\right |}
    \\\nonumber &\qquad+\,
    \sup_{x}{\left |\mathbb{P}_{n}\widetilde{\mathbb{P}}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,\star)-\mathbb{P}\widetilde{\mathbb{P}}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,\star)\right |}\,,
    \end{aligned}$$ which is further bounded by taking Lemma \[cant\_class\_rate\] into account: $$\begin{aligned}
    \nonumber &\,
    \sup_{f\in\mathcal{K}\cdot\mathcal{K}}{\left |\widetilde{\mathbb{P}}_{m}f - \widetilde{\mathbb{P}}f\right |} + 
    \sup_{f\in\int\mathcal{K}\cdot\mathcal{K}}{\left |\mathbb{P}_{n}f - {\mathbb{P}}f\right |}
    \\\nonumber =&\,
    \mathcal{O}\left(\frac{\sqrt{\log\frac{1}{\epsilon}}+\sqrt{\log m}}{\sqrt{m}}\right)+ \mathcal{O}\left(\frac{\sqrt{\log\frac{1}{\epsilon}}+\sqrt{\log n}}{\sqrt{n}}\right) =
    \mathcal{O}\left(\frac{\sqrt{\log\frac{1}{\epsilon}}+\sqrt{\log m}}{\sqrt{m}}\right)\end{aligned}$$ with probability $1-\mathcal{O}(m^{-2})-\mathcal{O}(n^{-2})=1-\mathcal{O}(m^{-2})$, where we use the fact that $m<n$. Then, since $d_{\textup{ref},\epsilon}(x)\geq C_1\epsilon^{d}$ by Lemma \[unif-bd\], by the assumption $\frac{\sqrt{\log\frac{1}{\epsilon}}+\sqrt{\log m}}{\sqrt{m}\epsilon^{d}}\rightarrow 0$, we have $\widehat{d}_{\textup{ref},\epsilon,n}(x)\geq C_1\epsilon^{d}/2$ with probability $1-\mathcal{O}(m^{-2})$ when $m$ is sufficiently large.  Likewise, we have $$\begin{aligned}
    \nonumber
    \nonumber&
    \sup_{x}{\left |d_{\textup{ref},\epsilon,n}(x) - d_{\textup{ref},\epsilon}(x)\right |}
    \\\nonumber =&\,
    \sup_{x}{\left |\mathbb{P}_{n}\widetilde{\mathbb{P}}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,\star)-\mathbb{P}\widetilde{\mathbb{P}}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,\star)\right |}
    \\\nonumber \leq&\,
 \sup_{f\in\int\mathcal{K}\cdot\mathcal{K}}{\left |\mathbb{P}_{n}f - \widetilde{\mathbb{P}}f\right |}=\mathcal{O}\left(\frac{\sqrt{\log\frac{1}{\epsilon}}+\sqrt{\log n}}{\sqrt{n}}\right)\,. 
    \end{aligned}$$ Again, by the assumption $\frac{\sqrt{\log\frac{1}{\epsilon}}+\sqrt{\log m}}{\sqrt{m}\epsilon^{d}}\rightarrow 0$, we have $d_{\textup{ref},\epsilon,n}(x)\geq C_1\epsilon^{d}/2$ with probability $1-\mathcal{O}(n^{-2})$ when $n$ is sufficiently large. 

Next, we control the other terms we need for the spectral convergence rate. 

\[Proposition: bound of several operators\] Suppose $m=n^\beta$ for $\beta\in (0,1]$ and take $\epsilon=\epsilon(n)$ so that $\epsilon\to 0$ when $n\to \infty$. If we further assume that $\frac{\sqrt{-\log\epsilon}+\sqrt{\log n†}}{\sqrt{n}\epsilon^{d}}\rightarrow 0$, we have with probability $1-\mathcal{O}(n^{-2})$ the following bound: $${\left\lVertT_{\textup{ref},\epsilon,n}\right\rVert}\leq \frac{2C_2}{C_1}\,,$$ where $C_1$ and $C_2$ are constants defined in Lemma \[unif-bd\].  Moreover, the following two bounds always hold: $${\left\lVertT_{\textup{ref},\epsilon}\right\rVert}\leq 1\,,\,\,{\left\lVert\widehat{T}_{\textup{ref},\epsilon,n}\right\rVert}\leq 1\,.$$

Take $f\in C(M)$ so that $\|f\|_\infty\leq 1$. Since the kernel is positive, by definition, $$\begin{aligned}
|\widehat{T}_{\textup{ref},\epsilon,n}f(x)|\leq \frac{\frac{1}{n}\sum_{i=1}^{n}\widehat{K}_{\textup{ref},\epsilon,n}(x,x_i)|f(x_{i})|}{\widehat{d}_{\textup{ref},\epsilon,n}(x)} \leq \|f\|_\infty\,,\nonumber\end{aligned}$$ $$|T_{\textup{ref},\epsilon}f(x)|\leq \frac{\int_{M}K_{\textup{ref},\epsilon}(x,y)|f(y)|p_{X}(y)\,dV(y)}{|d_{\textup{ref},\epsilon}(x)|}\leq \|f\|_\infty\,.\nonumber$$ Similarly, by Lemma \[d\_hat\_bd\], with probability $1-\mathcal{O}(n^{-2})$: $$\begin{aligned}
|T_{\textup{ref},\epsilon,n}f(x)|\leq \frac{\frac{1}{n}\sum_{i=1}^{n}K_{\textup{ref},\epsilon}(x,x_i)}{d_{\textup{ref},\epsilon}(x)}\|f\|_\infty= \frac{d_{\textup{ref},\epsilon,n}(x)}{d_{\textup{ref},\epsilon}(x)}\leq \frac{2C_2}{C_1} \,.\nonumber\end{aligned}$$ We hence finish the proof. 

\[tech\_bd\] Take $g\in C(M)$. Suppose $m=n^\beta$ for $\beta\in (0,1]$ and take $\epsilon=\epsilon(n)$ so that $\epsilon\to 0$ when $n\to \infty$. If we further assume that $\frac{\sqrt{-\log\epsilon}+\sqrt{\log m}}{\sqrt{m}\epsilon^{d}}\rightarrow 0$, when $n$ is sufficiently large, we have with probability $1-\mathcal{O}(m^{-2})$: $$\begin{aligned}
{\left\lVert\widehat{T}_{\textup{ref},\epsilon,n}-T_{\textup{ref},\epsilon,n}\right\rVert}&\,\leq \frac{2(\sqrt{-\log \epsilon}+\sqrt{\log m})}{C^2_{1}\sqrt{m}\epsilon^{2d}}\,,\label{Proposition SI3 bound Tnhat and Tn}\end{aligned}$$ and with probability $1-\mathcal{O}(n^{-2})$: $$\begin{aligned}
{\left\lVert(T_{\textup{ref},\epsilon}-T_{\textup{ref},\epsilon,n})T_{\textup{ref},\epsilon,n}\right\rVert}&\,\leq C_3\frac{\sqrt{-\log \epsilon}+\sqrt{\log n}}{\sqrt{n}}\,,\\
{\left\lVert(T_{\textup{ref},\epsilon,n}-T_{\textup{ref},\epsilon})g\right\rVert}_{\infty}&\,\leq C_4\frac{\sqrt{-\log \epsilon}+\sqrt{\log n}}{\sqrt{n}}\,,\nonumber\end{aligned}$$ where $C_3>0$ is a constant depending on the kernel, the curvature of the manifold and the minima of $p_X$ and $p_Y$, and the constants shown in the entropy bound shown in Theorem \[entropy\_bd\], and $C_4>0$ is a constant depending on the kernel, the curvature of the manifold and the minima of $p_X$ and $p_Y$, $\|g\|_\infty$ and the constants shown in the entropy bound shown in Theorem \[entropy\_bd\]. 

Note that the bound for ${\left\lVert\widehat{T}_{\textup{ref},\epsilon,n}-T_{\textup{ref},\epsilon,n}\right\rVert}$ is dominated by $\epsilon^{-2d}$ since we need to control the term $\widehat{d}_{\textup{ref},\epsilon,n}(x)$. Also, the condition $\frac{\sqrt{\log\frac{1}{\epsilon}}+\sqrt{\log m}}{\sqrt{m}\epsilon^{d}}\to 0$ does not imply that ${\left\lVert\widehat{T}_{\textup{ref},\epsilon,n}-T_{\textup{ref},\epsilon,n}\right\rVert}\to 0$.

Take $f\in C(M)$ so that $\|f\|_\infty\leq 1$. By the same calculation (e.g. Step 2 in the proof of Theorem \[spec\_cong\]) and Lemma \[d\_hat\_bd\], we have $$\begin{aligned}
        \nonumber&\, {\left\lVert(\widehat{T}_{\textup{ref},\epsilon,n}-T_{\textup{ref},\epsilon,n})f\right\rVert}_{\infty}
    =\max_{x\in M}{\left |\mathbb{P}_{n}\widehat{M}_{\textup{ref},\epsilon,n}(x,\cdot)f(\cdot) - \mathbb{P}_{n}M_{\textup{ref},\epsilon}(x,\cdot)f(\cdot)\right |}
    \\\nonumber\leq&\,\max_{x\in M}{\left |\mathbb{P}_{n}\widehat{M}_{\textup{ref},\epsilon,n}(x,\cdot)f(\cdot) -
    \mathbb{P}_{n}\widehat{M}^{(d)}_{\textup{ref},\epsilon,n}(x,\cdot)f(\cdot)\right |}
    \\\nonumber&\qquad+
    \max_{x\in M}{\left |\mathbb{P}_{n}\widehat{M}^{(d)}_{\textup{ref},\epsilon,n}(x,\cdot)f(\cdot)-\mathbb{P}_{n}M_{\textup{ref},\epsilon}(x,\cdot)f(\cdot)\right |}\,,
    \end{aligned}$$ where $\widehat{M}^{(d)}_{\textup{ref},\epsilon,n}(x,y):=\frac{K_{\textup{ref},\epsilon}(x,y)}{\widehat{d}_{\textup{ref},\epsilon,n}(x)}$. By , with probability $1-\mathcal{O}(m^{-2})$, we have $$\begin{aligned}
    \nonumber&\max_{x\in M}{\left |\mathbb{P}_{n}\widehat{M}_{\textup{ref},\epsilon,n}(x,\cdot)f(\cdot) -
    \mathbb{P}_{n}\widehat{M}^{(d)}_{\textup{ref},\epsilon,n}(x,\cdot)f(\cdot)\right |}\\
    \nonumber\leq&\, \frac{2\|f\|_\infty}{C_{1}\epsilon^{d}}\max_{x\in M}\frac{1}{n}\sum_{l=1}^n|K_{\textup{ref},\epsilon}(x,x_l)-\widehat K_{\textup{ref},\epsilon}(x,x_l)|\\
    \nonumber= &\, \frac{2\|f\|_\infty}{C_{1}\epsilon^{d}}\frac{1}{n}\sum_{l=1}^n\max_{x\in M}|\mathbb{P}_{m}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,x_l)-\widehat{\mathbb{P}}_{m}K_{\epsilon}(x,\cdot)K_{\epsilon}(\cdot,x_l) |\,. 
    \end{aligned}$$ We bound the right hand side by Lemma \[cant\_class\_rate\]; that is, with probability $1-\mathcal{O}(m^{-2})$, we have $$\begin{aligned}
    \nonumber&\, \frac{2\|f\|_\infty}{C_{1}\epsilon^{d}}\sup_{F\in\mathcal{K}\cdot\mathcal{K}}{\left |\widetilde{\mathbb{P}}_{m}F - \widetilde{\mathbb{P}}F\right |}\leq \frac{2C_3\|f\|_\infty }{C_{1}\epsilon^{d}}\frac{\sqrt{-\log \epsilon}+\sqrt{\log m}}{\sqrt{m}}\,,
    \end{aligned}$$ where $C_3$ is the implied constant in Lemma \[cant\_class\_rate\].  For the other term, with probability $1-\mathcal{O}(m^{-2})$, we have $$\begin{aligned}
    \nonumber&\sup_{x}{\left |\mathbb{P}_{n}\widehat{M}^{(d)}_{\textup{ref},\epsilon,n}(x,\cdot)f(\cdot)-\mathbb{P}_{n}M_{\textup{ref},\epsilon}(x,\cdot)f(\cdot)\right |}\\
    \nonumber\leq\,& \|f\|_\infty \max_{x\in M}\frac{1}{n}\sum_{l=1}^n\frac{K_{\textup{ref},\epsilon}(x,x_l)}{d_{\textup{ref},\epsilon}(x)\hat d_{\textup{ref},\epsilon,n}(x)}{\left |\hat d_{\textup{ref},\epsilon,n}(x)-d_{\textup{ref},\epsilon}(x)\right |}\\
    \nonumber\leq\,&    \frac{2\|f\|_\infty(\sqrt{-\log \epsilon}+\sqrt{\log m})}{C^2_{1}\sqrt{m}\epsilon^{2d}}\,,
    \end{aligned}$$ where we use the fact the $K_\epsilon$ is positive, $\|K_{\textup{ref},\epsilon}\|_\infty=1$, and Lemma \[d\_hat\_bd\].  As a result, by combining the above two bounds with a union probability bound, we conclude that when $m$ is sufficiently large, with probability $1-\mathcal{O}(m^{-2})$, $$\begin{aligned}
    {\left\lVert\widehat{T}_{\textup{ref},\epsilon,n}-T_{\textup{ref},\epsilon,n}\right\rVert} \nonumber\leq
    \frac{2(\sqrt{-\log \epsilon}+\sqrt{\log m})}{C^2_{1}\sqrt{m}\epsilon^{2d}}\,.\end{aligned}$$ The second statement follows the same argument, but with more terms to control: $$\begin{aligned}
    \nonumber&\quad        {\left\lVert(T_{\textup{ref},\epsilon}-T_{\textup{ref},\epsilon,n})T_{\textup{ref},\epsilon,n}f\right\rVert}_{\infty}
            \\\nonumber&=
    {\left\lVertT_{\textup{ref},\epsilon}\left(\frac{1}{n}\sum_{i=1}^{n}M_{\textup{ref},\epsilon}(x,x_{i})f(x_{i})\right)-T_{\textup{ref},\epsilon,n}\left(\frac{1}{n}\sum_{i=1}^{n}M_{\textup{ref},\epsilon}(x,x_{i})f(x_{i})\right)\right\rVert}_{\infty}
    \\\nonumber&=
    {\left\lVert\int M_{\textup{ref},\epsilon}(y,z)\left(\frac{1}{n}\sum_{i=1}^{n}M_{\textup{ref},\epsilon}(z,x_{i})f(x_{i})\right)d\mathbb{P}(z)\right.\\
    &\nonumber\left.\qquad\qquad\qquad-\frac{1}{n}\sum_{j=1}^{n}M_{\textup{ref},\epsilon}(y,z_{j})\left(\frac{1}{n}\sum_{i=1}^{n}M_{\textup{ref},\epsilon}(z_{j},x_{i})f(x_{i})\right)\right\rVert}_{\infty}\,,
    \end{aligned}$$ which, by noting that the term $\frac{1}{n}\sum_{i=1}^{n}f(x_{i})$ can be isolated, can be bounded by $$\begin{aligned}
   \nonumber&\quad
     {\left\lVertf\right\rVert}_\infty\sup_{y}{\left |\int M_{\textup{ref},\epsilon}(y,z)M_{\textup{ref},\epsilon}(z,x_{i})d\mathbb{P}(z)-\frac{1}{n}\sum_{j=1}^{n}M_{\textup{ref},\epsilon}(y,z_{j})M_{\textup{ref},\epsilon}(z_{j},x_{i})\right |}
    \\\nonumber&\leq
    \|f\|_\infty \sup_{F\in\mathcal{M}\cdot\mathcal{M}}{\left |\mathbb{P}F-\mathbb{P}_{n}F\right |}    \leq 
  C_3 \|f\|_\infty \frac{\sqrt{-\log \epsilon}+\sqrt{\log n}}{\sqrt{n}}$$ with probability $1-\mathcal{O}(n^{-2})$, where the first inequality comes from the fact that $\frac{1}{n}\sum_{i=1}^{n}{\left |f(x_{i})\right |}\leq 1$.

The final statement is by a direct bound: $$\begin{aligned}
    \nonumber{\left\lVert(T_{\textup{ref},\epsilon,n}-T_{\textup{ref},\epsilon})g\right\rVert}_{\infty} &= {\left\lVert\mathbb{P}_{n}M_{\textup{ref},\epsilon}(x,\cdot)f(\cdot) - \mathbb{P}M_{\textup{ref},\epsilon}(x,\cdot)f(\cdot)\right\rVert}_{\infty}
    \\\nonumber&\leq
    \sup_{F\in g\cdot\mathcal{M}}{\left |\mathbb{P}_{n}F - \mathbb{P}F\right |}\,,\end{aligned}$$ which leads to the conclusion by Lemma \[cant\_class\_rate\].  Note that due to the finite sampling, we cannot control the error simply by the $\|g\|_\infty$.

Based on the above preparation, the following proposition describes the spectral convergence of the operator $T_{\textup{ref},\epsilon,n}$ to the operator $T_{\textup{ref},\epsilon}$. Denote $\lambda_{i,\epsilon}$ is the $i$-th smallest eigenvalue of $\frac{I-T_{\textup{ref},\epsilon}}{\epsilon}$, and denote $u_{\epsilon,i}$ to be the associated eigenfunction.  Clearly, $\frac{I-T_{\textup{ref},\epsilon}}{\epsilon}$ and $T_{\textup{ref},\epsilon}$ share the same eigenfunctions.  Similarly, denote $\lambda_{\epsilon,n,i}$ to be the $i$-th smallest eigenvalue of $\frac{I-T_{\textup{ref},\epsilon,n}}{\epsilon}$, and denote $u_{\epsilon, n,i}$ to be the associated eigenfunction.  We assume that both $u_{\epsilon,i}$ and $u_{\epsilon, n,i}$ are normalized in the $L^2$ norm. 

\[Proposition: Tepsn conv to Teps with rate\] Fix $K\in \mathbb{N}$. Assume that the eigenvalues of $\Delta$ are simple.  Take $m=n^\beta$, where $\beta\in (0,1)$. Suppose $\epsilon=\epsilon(n)$ so that $\epsilon \to 0$ and $\frac{\sqrt{-\log\epsilon}+\sqrt{\log m}}{\sqrt{m}\epsilon^{d}} \rightarrow 0$, as $n \rightarrow \infty$, and $\epsilon \leq \mathcal{K}_1 \min \Bigg(\bigg(\frac{\min(\mathsf\Gamma_K,1)}{\mathcal{K}_2+\lambda_K^{d/2+5}}\bigg)^2, \frac{1}{(2+\lambda_K^{d+1})^2}\Bigg)$, where $\mathsf\Gamma_K$, $\mathcal{K}_1$ and $\mathcal{K}_2>1$ are introduced in Proposition \[T epsilon and Delta\], then there is a sequence $a_n \in \{1,-1\}$ such that with probability $1-\mathcal{O}(n^{-2})$, for all $i < K$, we have $$\begin{aligned}
& |\lambda_{\epsilon,n,i}-\lambda_{\epsilon,i}|\leq  3\mathcal{K}_3\frac{\sqrt{-\log \epsilon}+\sqrt{\log m}}{\sqrt{m}\epsilon^{2d+2}}\,, \nonumber  \\
& \|a_nu_{\epsilon,n,i}-u_{\epsilon,i}\|_{\infty} \leq \mathcal{K}_3\frac{\sqrt{-\log \epsilon}+\sqrt{\log m}}{\sqrt{m}\epsilon^{2d+2}}\,.  \nonumber \end{aligned}$$ where $\mathcal{K}_3$ is a constant depending on the kernel, the curvature of $M$, $p_X$ and $p_Y$.

Note that the imposed conditions, like $\frac{\sqrt{-\log\epsilon}+\sqrt{\log m}}{\sqrt{m}\epsilon^{d}}\to 0$, does not imply that $|\lambda_{\epsilon,n,i}-\lambda_{\epsilon,i}|\to 0$ or $\|a_nu_{\epsilon,n,i}-u_{\epsilon,i}\|_{\infty}\to 0$.

Note that due to Proposition \[T epsilon and Delta\] and the assumption that the eigenvalues of $\Delta$ are simple, we have that for any $K\in \mathbb{N}$, when $\epsilon>0$ is sufficiently small, the first smallest $K$ eigenvalues of $\frac{I-T_{\textup{ref},\epsilon}}{\epsilon}$ are simple.  Specifically, it is shown in the proof of [@DunsonWuWu2019 Proposition 1 (SI.20)] that when $\epsilon$ satisfies the assumption, for each $i<K$, we have $$\label{Spectral gap of I-Teps/eps}
\gamma_i\Big(\frac{I-T_{\textup{ref},\epsilon}}{\epsilon}\Big)\geq \frac{1}{12}\mathsf\Gamma_K\,,$$ where $\gamma_i$ is defined in . 

Fix $i<K$. Take $$r=\frac{\mathsf\Gamma_K}{24}\epsilon\,.$$ We now quantify how the sequence $\{u_{\epsilon,n,i}\}_{n=1}^\infty$ converges to $u_{\epsilon,i}$. It is clear that $\frac{I-T_{\textup{ref},\epsilon}}{\epsilon}$ and $T_{\textup{ref},\epsilon}$ share the same eigenfunctions, with the eigenvalues directly related.  Denote $\bar{\lambda}_{\epsilon,i}=1-\lambda_{\epsilon,i}\epsilon$ to be the $i$-th largest eigenvalue of $T_{\textup{ref},\epsilon}$. Similarly, this relationship holds for $\frac{I-T_{\textup{ref},\epsilon,n}}{\epsilon}$ and $T_{\textup{ref},\epsilon,n}$, and we denote $\bar{\lambda}_{\epsilon,n,i}=1-\lambda_{\epsilon,n,i}\epsilon$ to be the $i$-th largest eigenvalue of $T_{\textup{ref},\epsilon,n}$. Therefore, we can directly compare $T_{\textup{ref},\epsilon}$ and $T_{\textup{ref},\epsilon,n}$.

By Proposition \[cong\_1\_dim\_proj\], to control ${\left\lVerta_nu_{\epsilon,n,i}-u_{\epsilon,i}\right\rVert}$, we need to bound ${\left\lVertu_{\epsilon,i}-\textup{Pr}_{u_{\epsilon,n,i}}(u_{\epsilon,i})\right\rVert}$. Since $T_{\textup{ref},\epsilon,n}$ converges to $T_{\textup{ref},\epsilon}$ collectively compactly a.e.  by Lemma \[Proposition: Tepsn conv to Teps compactly\], we apply Theorem \[atkinson1967bound\] to control ${\left\lVertu_{\epsilon,i}-\textup{Pr}_{u_{\epsilon,n,i}}(u_{\epsilon,i})\right\rVert}$. To apply Theorem \[atkinson1967bound\], we need to control $\max_{z\in \Gamma_r(\bar\lambda_{\epsilon,i})} \|R_{z}(T_{\textup{ref},\epsilon})\|$, $\min_{z\in \Gamma_r(\bar\lambda_{\epsilon,i})}|z|$, ${\left\lVert(\widehat{T}_{\textup{ref},\epsilon,n}-T_{\textup{ref},\epsilon})u_{\epsilon,i}\right\rVert}$ and ${\left\lVert(T_{\textup{ref},\epsilon}-\widehat{T}_{\textup{ref},\epsilon,n})\widehat{T}_{\textup{ref},\epsilon,n}\right\rVert}$. First, according to , the spectral gap of the $i$-th largest eigenvalue of $T_{\textup{ref},\epsilon}$ is bounded from below by $\frac{\mathsf\Gamma_K}{12}\epsilon$. Therefore, by the basic bound of the resolvent (see, for example, [@DunsonWuWu2019 Lemma SI.16]) and the chosen $r$, we have $$\label{Proof spectral rate bound1}
\max_{z\in \Gamma_r(\bar\lambda_{\epsilon,i})} \|R_{z}(T_{\textup{ref},\epsilon})\|\leq \frac{1}{r}=\frac{24}{\mathsf\Gamma_K \epsilon}\,.$$ By Proposition \[T epsilon and Delta\], under the assumption, we have $\lambda_{\epsilon,i}\leq \lambda_i+\epsilon^{3/2}$, and hence if $\epsilon$ is sufficiently small, we have $$\label{Proof spectral rate bound2}
\min_{z\in \Gamma_r(\bar\lambda_{\epsilon,i})}|z|\geq 1/2\,.$$ For the remaining terms, by a direct triangular inequality, we have $$\begin{aligned}
    \nonumber{\left\lVert(\widehat{T}_{\textup{ref},\epsilon,n}-T_{\textup{ref},\epsilon})u_{\epsilon,i}\right\rVert}_\infty&\,\leq
   {\left\lVert(\widehat{T}_{\textup{ref},\epsilon,n}-T_{\textup{ref},\epsilon,n})u_{\epsilon,i}\right\rVert}_\infty+{\left\lVert(T_{\textup{ref},\epsilon,n}-T_{\textup{ref},\epsilon})u_{\epsilon,i}\right\rVert}_\infty\\
   \nonumber&\,\leq {\left\lVert\widehat{T}_{\textup{ref},\epsilon,n}-T_{\textup{ref},\epsilon,n}\right\rVert}{\left\lVertu_{\epsilon,i}\right\rVert}_\infty+{\left\lVert(T_{\textup{ref},\epsilon,n}-T_{\textup{ref},\epsilon})u_{\epsilon,i}\right\rVert}_\infty\,.\end{aligned}$$ Moreover, with probability $1-\mathcal{O}(m^{-2})$, $$\begin{aligned}
    \nonumber&{\left\lVert(T_{\textup{ref},\epsilon}-\widehat{T}_{\textup{ref},\epsilon,n})\widehat{T}_{\textup{ref},\epsilon,n}\right\rVert}
    \\\nonumber\leq&\,
    {\left\lVertT_{\textup{ref},\epsilon}\widehat{T}_{\textup{ref},\epsilon,n}-T_{\textup{ref},\epsilon}T_{\textup{ref},\epsilon,n}\right\rVert} + {\left\lVertT_{\textup{ref},\epsilon}T_{\textup{ref},\epsilon,n}-T_{\textup{ref},\epsilon,n}T_{\textup{ref},\epsilon,n}\right\rVert}
    \\\nonumber&\qquad+\,
    {\left\lVertT_{\textup{ref},\epsilon,n}T_{\textup{ref},\epsilon,n}-T_{\textup{ref},\epsilon,n}\widehat{T}_{\textup{ref},\epsilon,n}\right\rVert} + {\left\lVertT_{\textup{ref},\epsilon,n}\widehat{T}_{\textup{ref},\epsilon,n}-\widehat{T}_{\textup{ref},\epsilon,n}\widehat{T}_{\textup{ref},\epsilon,n}\right\rVert}
    \\\nonumber\leq&\,
    \left({\left\lVertT_{\textup{ref},\epsilon}\right\rVert}+{\left\lVertT_{\textup{ref},\epsilon,n}\right\rVert}+{\left\lVert\widehat{T}_{\textup{ref},\epsilon,n}\right\rVert}\right){\left\lVertT_{\textup{ref},\epsilon,n}-\widehat{T}_{\textup{ref},\epsilon,n}\right\rVert}\\
    &\nonumber\qquad+{\left\lVert(T_{\textup{ref},\epsilon}-T_{\textup{ref},\epsilon,n})T_{\textup{ref},\epsilon,n}\right\rVert}
    \\\nonumber\leq&\,
    C_{7}{\left\lVertT_{\textup{ref},\epsilon,n}-\widehat{T}_{\textup{ref},\epsilon,n}\right\rVert}+{\left\lVert(T_{\textup{ref},\epsilon}-T_{\textup{ref},\epsilon,n})T_{\textup{ref},\epsilon,n}\right\rVert}\,,\end{aligned}$$ for $C_7>0$, where the last bound comes from Lemma \[Proposition: bound of several operators\].  By Theorem \[atkinson1967bound\], we have: $$\begin{aligned}
    &\nonumber\quad{\left\lVertu_{\epsilon,i}-\textup{Pr}_{u_{\epsilon,n,i}}u_{\epsilon,i}\right\rVert}    \\\nonumber&\leq \max_{z\in \Gamma_r(\lambda_{\epsilon,i})} \frac{2r\|R_{z}(T_{\textup{ref},\epsilon})\|}{\min_{z\in \Gamma_r(\lambda_{\epsilon,i})}|z|}\Big({\left\lVert(\widehat{T}_{\textup{ref},\epsilon,n}-T_{\textup{ref},\epsilon})u_{\epsilon,i}\right\rVert} 
    \\\nonumber&\qquad\qquad\qquad\qquad+\|R_{z}(T_{\textup{ref},\epsilon})u_{\epsilon,i}\|_\infty{\left\lVert(T_{\textup{ref},\epsilon}-\widehat{T}_{\textup{ref},\epsilon,n})\widehat{T}_{\textup{ref},\epsilon,n}\right\rVert}\Big)
    \\&\leq 4\Big({\left\lVert(\widehat{T}_{\textup{ref},\epsilon,n}-T_{\textup{ref},\epsilon})\frac{u_{\epsilon,i}}{{\left\lVertu_{\epsilon,i}\right\rVert}_\infty}\right\rVert} +\frac{24}{\mathsf\Gamma_K \epsilon}{\left\lVert(T_{\textup{ref},\epsilon}-\widehat{T}_{\textup{ref},\epsilon,n})\widehat{T}_{\textup{ref},\epsilon,n}\right\rVert}\Big){\left\lVertu_{\epsilon,i}\right\rVert}_\infty  \nonumber \,,     \end{aligned}$$ where we plug in and .  Therefore, by plugging the bounds of ${\left\lVert(\widehat{T}_{\textup{ref},\epsilon,n}-T_{\textup{ref},\epsilon})u_{\epsilon,i}\right\rVert}$ and ${\left\lVert(T_{\textup{ref},\epsilon}-\widehat{T}_{\textup{ref},\epsilon,n})\widehat{T}_{\textup{ref},\epsilon,n}\right\rVert}$, with probability $1-\mathcal{O}(m^{-2})$, we have: $$\begin{aligned}
    \nonumber&{\left\lVertu_{\epsilon,i}-\textup{Pr}_{u_{\epsilon,n,i}}u_{\epsilon,i}\right\rVert}\\
    \nonumber\leq&\, 4\bigg[\Big(\frac{2}{C^2_{1}}+\frac{24C_7}{\mathsf \Gamma_K\epsilon}\Big)\frac{\sqrt{-\log \epsilon}+\sqrt{\log m}}{\sqrt{m}\epsilon^{2d}} +C_3\Big(1+\frac{24}{\mathsf\Gamma_K }\Big)\frac{\sqrt{-\log \epsilon}+\sqrt{\log n}}{\sqrt{n}\epsilon}\bigg]{\left\lVertu_{\epsilon,i}\right\rVert}_\infty\,.\end{aligned}$$ Note that as discussed after Lemma \[cant\_class\_rate\], when we apply Lemma \[cant\_class\_rate\], the bound depends on the eigenfunction.  To control ${\left\lVertu_{\epsilon,i}\right\rVert}_\infty$, note that by Proposition \[T epsilon and Delta\] and Lemma \[lemma hormander\], we have $$\|u_{\epsilon,i}\|_\infty\leq \|u_i\|_\infty+\epsilon\leq C_1\lambda_K^{(d-1)/4}+\epsilon\leq 2C_1\lambda_K^{(d-1)/4}\,,$$ where the last inequality comes from the assumption of $\epsilon$. Moreover, by the assumption of $\epsilon$, we have $$\label{controls of u infty u/Gamma and 1/Gamma}
\max\Big\{\|u_{\epsilon,i}\|_\infty,\frac{\|u_{\epsilon,i}\|_\infty}{\mathsf\Gamma_K},\frac{1}{\mathsf\Gamma_K}\Big\}\leq \epsilon^{-1}\,.$$ As a result, with probability $1-\mathcal{O}(m^{-2})$, we can find $a_n\in \{1,-1\}$ so that $$\begin{aligned}
&{\left\lVerta_nu_{\epsilon,n,i}-u_{\epsilon,i}\right\rVert}\leq 2{\left\lVertu_{\epsilon,i}-\textup{Pr}_{u_{\epsilon,n,i}}u_{\epsilon,i}\right\rVert}_\infty
\leq 192C_7\frac{\sqrt{-\log \epsilon}+\sqrt{\log m}}{\sqrt{m}\epsilon^{2d+2}}\,.\label{uepsni and uepsi difference L infty}\end{aligned}$$ By setting $\mathcal{K}_3:=192C_7$, we get the claim for the eigenvectors.  For eigenvalues, we have $$\begin{aligned}
    \nonumber&{\left |\bar\lambda_{\epsilon,i}-\bar\lambda_{\epsilon,n,i}\right |}{\left\lVertu_{\epsilon,i}\right\rVert}_\infty \\
    =&\, \nonumber{\left\lVert\bar\lambda_{\epsilon,i}u_{\epsilon,i}-\bar\lambda_{\epsilon,n,i}u_{\epsilon,i}\right\rVert}_\infty
    \\\nonumber\leq&\,
    {\left\lVert\bar\lambda_{\epsilon,i}u_{\epsilon,i}-\bar\lambda_{\epsilon,n,i}a_{n}u_{\epsilon,n,i}\right\rVert}_\infty + {\left |\bar\lambda_{\epsilon,n,i}\right |}{\left\lVerta_{n}u_{\epsilon,n,i}-u_{\epsilon,i}\right\rVert}_\infty
    \\\nonumber=&\,
    {\left\lVertT_{\textup{ref},\epsilon}u_{\epsilon,i}-a_{n}\widehat{T}_{\textup{ref},\epsilon,n}u_{\epsilon,n,i}\right\rVert}_\infty + {\left |\bar\lambda_{\epsilon,n,i}\right |}{\left\lVerta_{n}u_{\epsilon,n,i}-u_{\epsilon,i}\right\rVert}_\infty\,,
    \end{aligned}$$ which we further bound by $$\begin{aligned}
    \nonumber\quad&\,
    {\left\lVertT_{\textup{ref},\epsilon}u_{\epsilon,i}-\widehat{T}_{\textup{ref},\epsilon,n}u_{\epsilon,i}\right\rVert}+{\left\lVert\widehat{T}_{\textup{ref},\epsilon,n}u_{\epsilon,i}-a_{n}\widehat{T}_{\textup{ref},\epsilon,n}u_{\epsilon,n,i}\right\rVert}_\infty + {\left |\bar\lambda_{\epsilon,n,i}\right |}{\left\lVerta_{n}u_{\epsilon,n,i}-u_{\epsilon,i}\right\rVert}_\infty
    \\\nonumber\leq &\,
    {\left\lVert(T_{\textup{ref},\epsilon,n}-\widehat{T}_{\textup{ref},\epsilon,n})u_{\epsilon,i}\right\rVert}_\infty +{\left\lVert(T_{\textup{ref},\epsilon}-{T}_{\textup{ref},\epsilon,n})u_{\epsilon,i}\right\rVert}_\infty \\
   \nonumber &\qquad+\left({\left\lVert\widehat{T}_{\textup{ref},\epsilon,n}\right\rVert}+{\left |\bar\lambda_{\epsilon,n,i}\right |}\right){\left\lVerta_{n}u_{\epsilon,n,i}-u_{\epsilon,i}\right\rVert}_\infty
\\\nonumber\leq &\,
    \Big({\left\lVertT_{\textup{ref},\epsilon,n}-\widehat{T}_{\textup{ref},\epsilon,n}\right\rVert}+{\left\lVert(T_{\textup{ref},\epsilon}-{T}_{\textup{ref},\epsilon,n})\frac{u_{\epsilon,i}}{{\left\lVertu_{\epsilon,i}\right\rVert}_\infty}\right\rVert}_\infty\Big) {\left\lVertu_{\epsilon,i}\right\rVert}_\infty\\
   \nonumber &\qquad+2{\left\lVert\widehat{T}_{\textup{ref},\epsilon,n}\right\rVert}{\left\lVerta_{n}u_{\epsilon,n,i}-u_{\epsilon,i}\right\rVert}_\infty\,,
       \end{aligned}$$ where we use the fact that $\frac{u_{\epsilon,i}}{{\left\lVertu_{\epsilon,i}\right\rVert}_\infty}$ has the $L^\infty$ norm $1$, and ${\left |\bar\lambda_{\epsilon,n,i}\right |}\leq {\left\lVert\widehat{T}_{\textup{ref},\epsilon,n}\right\rVert}$. As a result, by Lemmas \[Proposition: bound of several operators\] and \[tech\_bd\], and , with probability $1-\mathcal{O}(m^{-2})$, ${\left |\bar\lambda_{\epsilon,i}-\bar\lambda_{\epsilon,n,i}\right |}{\left\lVertu_{\epsilon,i}\right\rVert}_\infty$ is bounded by $$\begin{aligned}
\nonumber 2C_{3}\frac{\sqrt{-\log \epsilon}+\sqrt{\log n}}{\sqrt{n}\epsilon}+ 2\mathcal{K}_3\frac{\sqrt{-\log \epsilon}+\sqrt{\log m}}{\sqrt{m}\epsilon^{2d+2}}\,.\end{aligned}$$ Finally, note that since $\beta<1$, $\frac{\sqrt{-\log \epsilon}+\sqrt{\log n}}{\sqrt{n}\epsilon}$ is always dominated by $\frac{\sqrt{-\log \epsilon}+\sqrt{\log m}}{\sqrt{m}\epsilon^{2d+2}}$. We thus finish the proof. 

Finish the proof of Theorem \[spec\_cong\]
------------------------------------------

With the above preparation, we are ready to prove the main theorem. 

By Lemma \[Proposition: Tepsn conv to Teps compactly\], $\widehat{T}_{\textup{ref},\epsilon,n}$ converges to $T_{\textup{ref},\epsilon}$ compactly $\textup{a.s.}$ as $n\rightarrow\infty$. Therefore, Proposition \[Proposition: Tepsn conv to Teps with rate\] leads to the spectral convergence of $\widehat{T}_{\textup{ref},\epsilon,n}$ to $T_{\textup{ref},\epsilon}$ with the rate.  Next, we link $T_{\textup{ref},\epsilon}$ to $-\Delta$. By Theorem \[Proof: Theorem1: bias analysis\], when $p_Y$ is properly chosen so that $\frac{2\nabla p_{X}(x)}{p_{X}(x)}+\frac{\nabla p_{Y}(x)}{p_{Y}(x)}=0$, we have the convergence of the eigenvalue/eigenfunction of $\frac{1-T_{\textup{ref},\epsilon}}{\epsilon}$ to those of $-\Delta$. Thus, by Proposition \[T epsilon and Delta\], we have the spectral convergence of $\frac{1-T_{\textup{ref},\epsilon}}{\epsilon}$ to $-\Delta$ with the rate.  Finally, we put all the above together and finish the spectral convergence proof. 

More relevant work {#Appendix:MoreRelatedWork}
==================

There are some works focusing on accelerating the spectral clustering.  In Pham and Chen’s work [@pham2018large], the authors propose to use a smaller subset as landmarks, and take the co-clustering idea into account for the spectral clustering purpose.  The landmarks are chosen by uniform sampling or by the k-means clustering.  Although this landmark idea is the same to our proposed algorithm at the first glance, they are different since the affinity matrix between data and landmarks are normalized differently.  Specifically, the authors normalize the affinity matrix by both row and column sums.  This approach can thus be classified as the third category.  In Yan, Huang and Jordan [@yan2009fast], the authors apply the k-means or random projection tree first to pre-group the data sets and obtain $m$ centroids and hence $m$ groups.  Then, apply the spectral clustering on these $m$ centroids to obtain $k$ new centroids.  Finally, merge the original $m$ groups according to these $k$ centroids.  This is thus also classified as the third category.  Although these approaches do not focus on recovering the intrinsic geometry, like the geodesic distance or almost isometric embedding, we still consider them as relative studies since the spectral clustering is directly related to the GL.  Indeed, under the manifold setup, the spectral clustering can be understood as finding the connected components, which is related to recovering the null space of the Laplace-Beltrami operator.  Moreover, the behavior of those algorithms on the manifold setup is not clear. 

Salhov, et.al. [ @salhov2015approximately] propose a dictionary-based method, where they incrementally construct an approximative map by using a single scan of the data.  This algorithm is greedy, inefficient and sensitive to the scan order, as is commented in their followup paper [@bermanis2016incomplete].  Bermanis, et.al. [ @bermanis2016incomplete] provide an incomplete pivoted QR-based deterministic method for dimensionality reduction [*after*]{} running the DM.  In general, given a generic matrix $A$, the authors select a subset of columns of $A$ called pivots, and then use them to perform an incomplete QR factorization to approximate $A$.  This method is applied to reduce the dimension of the DM, and is claimed to preserve the data geometry up to a user-specified distortion rate.  Based on the nature of this algorithm, while it looks relevant at the first glance, this approach is not directly related to our acceleration mission. 

Another useful and closely related algorithms are the subspace sampling [@drineas2006subspace] and the CUR decomposition [@mahoney2009cur; @wang2013improving].  Those algorithms aim to find interpretable $k$ rank approximation of the original data matrix with respect to the matrix Frobenius norm.  For example, the SVD would give us [@mahoney2009cur]: $(1/2)\textup{age}- (1/\sqrt{2})\textup{height}+(1/2)\textup{income}$. This kind of linear combination of uncorrelated features are difficult to interpret in some situations.  The CUR on the other hand, is interpretable as it decomposes the original data matrix into a small number of actual columns and rows.  In general, the variable selection algorithms, like LASSO, can also be considered as dimension reduction algorithms, although the purpose of variable selection is totally different and the computational complexity might not be suitable for our purpose.  However, from the aspect of reducing the dimension of the dataset, in the broad sense they are closely related to algorithms in the first category. 

In [@cheng2005compression], the authors give a new decomposition form to compress a rank-deficient matrices when the SVD cannot be used efficiently.  Specifically, a matrix $A$ of rank $k$ is represented as $A=UBV$, where $B$ is a $k\times k$ submatrix of $A$, and $U,V$ are well-conditioned matrices that each contain a $k \times k$ identity submatrix.  The geometric interpretation of this decomposition is: columns (rows) of A are expressed as linear combinations of $k$ selected columns (rows) of A, the selection induces the matrix $B$, such that in the new coordinate system, the action of $A$ is represented by the action of $B$.

More numerical results {#More numerical results}
======================

Visualization of torus {#More numerical results torus}
----------------------

We consider the torus $\mathbb{T}^{2}$ embedded in $\mathbb{R}^3$ via $(\theta,\phi)\mapsto \begin{bmatrix}(0.8 + 0.4  \cos\theta) \cos\phi & (0.8 + 0.4  \cos\theta) \sin\phi & 0.6 \sin\theta\end{bmatrix}^\top$, where $\theta,\phi\in [0,2\pi)$.

We uniformly sample (with respect to the metric on $T^{2}$) $5,000$ points from $T^{2}$ to be our dataset and $200$ points to be the landmark set.  No noise is added to the dataset.  To enhance the visualization, in addition to plotting the embedding, we also display the embedding by cutting them vertically. 

![ First row from left to right: the clean Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the clean Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](cleanT2.pdf "fig:"){width="2.4cm"} ![ First row from left to right: the clean Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the clean Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](cleanT2_DM.pdf "fig:"){width="2.4cm"} ![ First row from left to right: the clean Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the clean Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](cleanT2_nys.pdf "fig:"){width="2.4cm"} ![ First row from left to right: the clean Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the clean Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](cleanT2_Haddad.pdf "fig:"){width="2.4cm"}![First row from left to right: the clean Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the clean Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](cleanT2_ref.pdf "fig:"){width="2.4cm"}\
![ First row from left to right: the clean Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the clean Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](cleanT2_Vcut.pdf "fig:"){width="2.4cm"} ![ First row from left to right: the clean Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the clean Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](cleanT2_DM_Vcut.pdf "fig:"){width="2.4cm"} ![ First row from left to right: the clean Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the clean Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](cleanT2_nys_Vcut.pdf "fig:"){width="2.4cm"} ![ First row from left to right: the clean Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the clean Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](cleanT2_Haddad_Vcut.pdf "fig:"){width="2.4cm"}![First row from left to right: the clean Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the clean Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](cleanT2_ref_Vcut.pdf "fig:"){width="2.4cm"}

Next, to show the robustness of the Roseland for the visualization purpose, we compare these methods when the dataset is corrupted by noises.  We embed the torus to the first three coordinates of $\mathbb{R}^{p}$, where $p=100$. That is, we first uniformly sample $5,000$ points from $\mathbb{T}^{2}$ to be our dataset and $200$ points to be the landmark set in $\mathbb{R}^{3}$ just as the above.  Then, we embed these clean datasets to $\mathbb{R}^{p}$ by appending zeros to both the dataset and the landmark set.  Next we keep the landmark set clean and add Gaussian noise $\epsilon_{i}$ to the dataset, where $\epsilon_{i}$ are iid from $\mathcal{N}(0, \frac{1}{\sqrt{p}} I_{p\times p})$. To enhance the visualization, in addition to plotting the embedding, we also display the embedding by cutting them vertically. 

![ First row from left to right: the noisy Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the noisy Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](noisyT2.jpg "fig:"){width="2.4cm"} ![ First row from left to right: the noisy Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the noisy Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](noisyT2_DM.pdf "fig:"){width="2.4cm"} ![ First row from left to right: the noisy Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the noisy Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](noisyT2_nys.pdf "fig:"){width="2.4cm"} ![ First row from left to right: the noisy Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the noisy Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](noisyT2_Haddad.pdf "fig:"){width="2.4cm"}![First row from left to right: the noisy Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the noisy Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](noisyT2_ref.pdf "fig:"){width="2.4cm"}\
![ First row from left to right: the noisy Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the noisy Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](noisyT2_Vcut.pdf "fig:"){width="2.4cm"} ![ First row from left to right: the noisy Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the noisy Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](noisyT2_DM_Vcut.pdf "fig:"){width="2.4cm"} ![ First row from left to right: the noisy Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the noisy Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](noisyT2_nys_Vcut.pdf "fig:"){width="2.4cm"} ![ First row from left to right: the noisy Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the noisy Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](noisyT2_Haddad_Vcut.pdf "fig:"){width="2.4cm"} ![ First row from left to right: the noisy Torus dataset and clean subset (in red), where the data size is 5,000 and the landmark set size is $200$, the DM embedding, the embedding by the Nyström extension, the HKC embedding, and the Roseland embedding.  Second row from left to right: the noisy Torus dataset cut vertically, the DM embedding cut vertically, the Nyström extension cut vertically, the HKC embedding cut vertically, and the Roseland embedding cut vertically.](noisyT2_ref_Vcut.pdf "fig:"){width="2.4cm"}

Clean data set and clean subsets with various sizes
---------------------------------------------------

We investigate the influence of the size of the subset.  Again, we sample 10,000 points uniformly from $S^{1}$ to be the dataset and vary the subset size from 100 to 1,000.  The result is shown in Figure \[Fig: Various size comparison\]. 

![ Clean data set and clean subsets with various sizes.  Left: relative $L^{\infty}$ and $L^{2}$ error of the $1^{\textup{st}}$ non-trivial eigenvector by Nyström and Roseland with ground truth using different subset sizes.  Right: relative $L^{\infty}$ and $L^{2}$ error of the $10^{\textup{th}}$ non-trivial eigenvector by Nyström and Roseland with ground truth using different subset sizes.\[Fig: Various size comparison\]](cleanS1_clean_sub_vary_err_eigenvec1.pdf "fig:"){width="6cm"} ![ Clean data set and clean subsets with various sizes.  Left: relative $L^{\infty}$ and $L^{2}$ error of the $1^{\textup{st}}$ non-trivial eigenvector by Nyström and Roseland with ground truth using different subset sizes.  Right: relative $L^{\infty}$ and $L^{2}$ error of the $10^{\textup{th}}$ non-trivial eigenvector by Nyström and Roseland with ground truth using different subset sizes.\[Fig: Various size comparison\]](cleanS1_clean_sub_vary_err_eigenvec10.pdf "fig:"){width="6cm"}

Clean data set and clean subset
-------------------------------

We sample $n=10,000$ points from the $S^1$ model without noise to be the dataset, and independently sample another $m=500$ points uniformly to be the landmark set for the Roseland.  For a fair comparison, we follow the same rule and take the landmark set to be the input to the Nyström extension and as the reference set of the HKC.  The results are shown in Figures \[S1clean\_eigenvalue\] and \[S1clean\_eigenvector\].  Clearly, when the data is clean, the Nyström extension, the HKC and the Roseland all can successfully recover first few eigenfunctions of the Laplace-Beltrami operator, while the recovery by the Nyström extension deviates a bit.  This comparison is quantified by the relative error shown in Figure \[S1clean\_eigenvector\].  Note that since the HKC is not designed to recover the Laplace-Beltrami operator, it is not clear how to estimate the eigenvalues, so we do not compare the eigenvalues.  But simply compare the eigenfunctions. 

![ Illustration of the Nyström method and the Roseland on the clean dataset and landmark set.  We superimpose the top 12 non-trivial eigenvectors determined by the Nyström extension and the Roseland (in blue) superimposed on the ground truth (in red).  Top two rows: the first 6 eigenvectors; bottom two rows: the $7^{\textup{th}}$ to the $12^{\textup{th}}$ eigenvectors.\[S1clean\_eigenvector\]](cleanS1_clean_sub_eigvec1-6.jpg "fig:"){width="12cm"} ![ Illustration of the Nyström method and the Roseland on the clean dataset and landmark set.  We superimpose the top 12 non-trivial eigenvectors determined by the Nyström extension and the Roseland (in blue) superimposed on the ground truth (in red).  Top two rows: the first 6 eigenvectors; bottom two rows: the $7^{\textup{th}}$ to the $12^{\textup{th}}$ eigenvectors.\[S1clean\_eigenvector\]](cleanS1_clean_sub_eigvec7-12.jpg "fig:"){width="12cm"}

![ Illustration of the Nyström method and the Roseland on the clean dataset and clean landmark set.  Left: superimpose the top 18 non-trivial eigenvalues by the Nyström extension and the Roseland, along with the ground truth.  Middle: relative error of eigenvalues.  Right: Relative $L^{\infty}$ and $L^{2}$ error of the top 12 non-trivial eigenvectors by the Nyström extension and the Roseland.  \[S1clean\_eigenvalue\]](cleanS1_clean_sub_eigval.pdf "fig:"){width="3.6cm"} ![ Illustration of the Nyström method and the Roseland on the clean dataset and clean landmark set.  Left: superimpose the top 18 non-trivial eigenvalues by the Nyström extension and the Roseland, along with the ground truth.  Middle: relative error of eigenvalues.  Right: Relative $L^{\infty}$ and $L^{2}$ error of the top 12 non-trivial eigenvectors by the Nyström extension and the Roseland.  \[S1clean\_eigenvalue\]](cleanS1_clean_sub_eigval_err.pdf "fig:"){width="3.6cm"} ![ Illustration of the Nyström method and the Roseland on the clean dataset and clean landmark set.  Left: superimpose the top 18 non-trivial eigenvalues by the Nyström extension and the Roseland, along with the ground truth.  Middle: relative error of eigenvalues.  Right: Relative $L^{\infty}$ and $L^{2}$ error of the top 12 non-trivial eigenvectors by the Nyström extension and the Roseland.  \[S1clean\_eigenvalue\]](cleanS1_clean_sub_err.pdf "fig:"){width="5cm"}

Noisy data set and clean subset
-------------------------------

We now show that we can use few clean data points as the landmark set to “denoise” the data set by the Roseland.  We uniformly sample $n=90,000$ points from the $S^1$ to be the dataset and independently sample another $m=300$ points uniformly to be the landmark set for the Roseland; that is, $\beta=0.5$. To make it more challenging, we first embed all points to $\mathbb{R}^{p}$, where $p=100$. Then, we keep the landmark set clean and add independent Gaussian noise $\epsilon_{i}$ to the dataset, where $\epsilon_{i}$ are i.i.d.  sampled from $\mathcal{N}(0, \frac{1}{\sqrt{p}} I_{p\times p})$.

The results are shown in Figures \[S1noisy\_embedding\] , \[S1noisy\_eigenvalue\] and \[S1noisy\_eigenvector\].  Clearly, while the Nyström method embeds $S^1$ successfully, the embedding by the HKC is cleaner, and the embedding by the Roseland is the cleanest.  Again, the Nyström method can only successfully recover first few eigenfunctions of the Laplace-Beltrami operator, but they are noisier.  The reason we have “seemingly” similar eigenvectors but different embeddings is the same as the previous section. 

![ Left: the noisy data and clean landmark set.  Middle left: the Nyström embedding.  Middle right: the HKC embedding.  Right: the Roseland embedding.\[S1noisy\_embedding\] ](noisy_S1_clean_subset.jpg "fig:"){width="3cm"} ![ Left: the noisy data and clean landmark set.  Middle left: the Nyström embedding.  Middle right: the HKC embedding.  Right: the Roseland embedding.\[S1noisy\_embedding\] ](noisy_S1_clean_subset_nys.jpg "fig:"){width="3cm"} ![ Left: the noisy data and clean landmark set.  Middle left: the Nyström embedding.  Middle right: the HKC embedding.  Right: the Roseland embedding.\[S1noisy\_embedding\] ](noisy_S1_clean_subset_haddad.jpg "fig:"){width="3cm"}![Left: the noisy data and clean landmark set.  Middle left: the Nyström embedding.  Middle right: the HKC embedding.  Right: the Roseland embedding.\[S1noisy\_embedding\] ](noisy_S1_clean_subset_ref.pdf "fig:"){width="3cm"}

![ Left: superimpose the top 18 non-trivial eigevalues by the Nyström extension and the Roseland with the ground truth.  Middle: relative errors of eigenvalues.  Right: relative $L^{\infty}$ and $L^{2}$ errors of the top 12 non-trivial eigenvectors by the Nyström extension, the HKC and the Roseland with the ground truth.  \[S1noisy\_eigenvalue\]](noisyS1_clean_sub_eigval.pdf "fig:"){width="4cm"} ![ Left: superimpose the top 18 non-trivial eigevalues by the Nyström extension and the Roseland with the ground truth.  Middle: relative errors of eigenvalues.  Right: relative $L^{\infty}$ and $L^{2}$ errors of the top 12 non-trivial eigenvectors by the Nyström extension, the HKC and the Roseland with the ground truth.  \[S1noisy\_eigenvalue\]](noisyS1_clean_sub_eigval_err.pdf "fig:"){width="4cm"} ![ Left: superimpose the top 18 non-trivial eigevalues by the Nyström extension and the Roseland with the ground truth.  Middle: relative errors of eigenvalues.  Right: relative $L^{\infty}$ and $L^{2}$ errors of the top 12 non-trivial eigenvectors by the Nyström extension, the HKC and the Roseland with the ground truth.  \[S1noisy\_eigenvalue\]](noisyS1_clean_sub_err.pdf "fig:"){width="4.2cm"}

![ Illustration of the Nyström method and the Roseland on the noisy dataset and clean landmark set.  We superimpose the top 12 non-trivial eigenvectors determined by the Nyström extension, the HKC and the Roseland (in blue) superimposed on the ground truth (in red).  Top three rows: the first 6 eigenvectors; bottom three rows: the $7^{\textup{th}}$ to the $12^{\textup{th}}$ eigenvectors.\[S1noisy\_eigenvector\]](noisyS1_clean_sub_eigvec1-6.jpg "fig:"){width="12cm"} ![ Illustration of the Nyström method and the Roseland on the noisy dataset and clean landmark set.  We superimpose the top 12 non-trivial eigenvectors determined by the Nyström extension, the HKC and the Roseland (in blue) superimposed on the ground truth (in red).  Top three rows: the first 6 eigenvectors; bottom three rows: the $7^{\textup{th}}$ to the $12^{\textup{th}}$ eigenvectors.\[S1noisy\_eigenvector\]](noisyS1_clean_sub_eigvec7-12.jpg "fig:"){width="12cm"}

From the plots we can see that while Nyström extension also did a good job in top eigenvectors or low frequency information, the Roseland is better on extracting higher frequency information. 

[^1]: Theoretically, it can reach ${O}(N^{\omega+\eta})$, where the $N^{\omega}$ part comes from the algorithm of matrix multiplication, and an arbitrary $\eta>0$ [@Demmel2007].  Note that when $M$ is dense, $\omega=\omega_0\approx 2.376$ [@Coppersmith1990].  However, the implied constant in these asymptotic is too large and cannot be practical [@le2012faster].