Diffusion $K$-means clustering on manifolds: provable exact recovery via semidefinite relaxations

TL;DR

This paper introduces a diffusion K-means clustering method on manifolds, utilizing semidefinite relaxations to achieve provable exact recovery under certain geometric conditions, suitable for complex non-linear data.

Contribution

It proposes a novel diffusion K-means algorithm on manifolds with convex relaxation via SDP, including an adaptive version, and provides theoretical guarantees for exact recovery.

Findings

Exact recovery under suitable manifold separability conditions

Convex SDP relaxations efficiently solve the clustering problem

Adaptive localized diffusion K-means adapts to local data structures

Abstract

We introduce the {\it diffusion $K$-means} clustering method on Riemannian submanifolds, which maximizes the within-cluster connectedness based on the diffusion distance. The diffusion $K$-means constructs a random walk on the similarity graph with vertices as data points randomly sampled on the manifolds and edges as similarities given by a kernel that captures the local geometry of manifolds. The diffusion $K$-means is a multi-scale clustering tool that is suitable for data with non-linear and non-Euclidean geometric features in mixed dimensions. Given the number of clusters, we propose a polynomial-time convex relaxation algorithm via the semidefinite programming (SDP) to solve the diffusion $K$-means. In addition, we also propose a nuclear norm regularized SDP that is adaptive to the number of clusters. In both cases, we show that exact recovery of the SDPs for diffusion $K$-means can be achieved under suitable between-cluster separability and within-cluster connectedness of the submanifolds, which together quantify the hardness of the manifold clustering problem. We further propose the {\it localized diffusion $K$-means} by using the local adaptive bandwidth estimated from the nearest neighbors. We show that exact recovery of the localized diffusion $K$-means is fully adaptive to the local probability density and geometric structures of the underlying submanifolds.

Tables

Table 1. Table 1. ℓ 1 superscript ℓ 1 \ell^{1} estimation errors of the SDP solutions of various diffusion K 𝐾 K -means clustering methods.

Method	Estimation error
	DGP=1	DGP=2	DGP=3
DKM	$4.7642\times {10}^{-6}$	$3.3258\times {10}^{-4}$	–
LDKM	$5.2835\times {10}^{-5}$	0.0049	0.0451

Table 2. Table 2. Classification errors of various diffusion K 𝐾 K -means and spectral clustering methods.

Method	Classification error
	DGP=1	DGP=2	DGP=3	DGP=3’
DKM	0	$1.3021\times {10}^{-4}$	–	–
SC-UN	0	0	–	–
SC-RWN	0	0	–	–
SC-NJW	0	0	–	–
LDKM	0	0.0018	0.0086	0.0594
LSC-UN	0	$4.7917\times {10}^{-4}$	0.0098	0.0801
LSC-RWN	0	0.0016	0.0105	0.0802
LSC-NJW	0	0.0399	0.0084	0.0884

Table 3. Table 3. Percentages of correctly estimated number of clusters by the regularized diffusion K 𝐾 K -means and its local scaling version using Algorithm  1 among 1,000 simulations.

Method	DGP=1	DGP=2	DGP=3
RDKM	94.30%	95.10%	–
LRDKM	99.20%	83.10%	97.70%

Equations

$$\kappa (x,y)=\exp \left(-\frac{\Vert x-y{\Vert }^2}{2h^2}\right),$$

$$d(x)={\int }_S\kappa (x,y)\mathrm{d}\mu (y)$$

$$p(x,y)=\frac{\kappa (x,y)}{d(x)},$$

$${\int }_{S}p(x,y)\mathrm{d}\mu (y)=1.$$

$$\mathcal{P}f(x)={\int }_{S}p(x,y)f(y)\mathrm{d}\mu (y).$$

$${\mathcal{P}}^{t}f(x)={\int }_S{p}_t(x,y)f(y)\mathrm{d}\mu (y),t=1,2,\dots ,$$

$$\pi (x)=\frac{\mathrm{d}\Pi }{\mathrm{d}\mu }(x)=\frac{d(x)}{{\int }_{S}d(y)\mathrm{d}\mu (y)}.$$

$$\Pi ({\mathcal{P}}^{t}f)=\Pi (f)$$

$$\pi (x)p(x,y)=\pi (y)p(y,x),\forall x,y\in S.$$

$$\mathcal{R}(x,y)=\frac{\kappa (x,y)}{\sqrt{d(x)}\sqrt{d(y)}},\forall x,y\in S.$$

$${\int }_S{\int }_S\mathcal{R}(x,y{)}^2\mathrm{d}\mu (x)\mathrm{d}\mu (y)<\infty ,$$

$$\mathcal{R}(x,y)=\sum _{j=0}^{\infty }{\lambda }_j{\varphi }_j(x){\varphi }_j(y),$$

$$p(x,y)=\sum _{j=0}^{\infty }{\lambda }_j{\psi }_j(x){\tilde{\psi }}_j(y),$$

$$\mathcal{P}{\psi }_j={\lambda }_j{\psi }_j,j=0,1,\dots .$$

$${\mathcal{D}}_t(x,y):=\Vert p_t(x,\cdot )-p_t(y,\cdot ){\Vert }_{L^2(\mathrm{d}\mu /d)}={\left\{{\int }_S[p_t(x,z)-p_t(y,z){]}^2\frac{\mathrm{d}\mu (z)}{d(z)}\right\}}^{1/2}.$$

$${\mathcal{D}}_t^2(x,y)=\sum _{j=0}^{\infty }{\lambda }_j^{2t}[{\psi }_j(x)-{\psi }_j(y){]}^2$$

$${\mathcal{D}}_t(x,y)={\left\{\sum _{j=1}^{\infty }{\lambda }_j^{2t}[{\psi }_j(x)-{\psi }_j(y){]}^2\right\}}^{1/2}.$$

$${\Psi }_t^{(q)}(x)=({\lambda }_0^t{\psi }_0(x),{\lambda }_1^t{\psi }_1(x),\dots ,{\lambda }_q^t{\psi }_q(x){)}^T,$$

$${\mathcal{D}}_t(x,y)=\underset{q\to \infty }{\lim }\Vert {\Psi }_t^{(q)}(x)-{\Psi }_t^{(q)}(y){\Vert }_2.$$

$$d_n(x)=n{\int }_{S_n}\kappa (x,y)\mathrm{d}{\mu }_n(y)=\sum _{i=1}^n\kappa (x,X_i),\forall x\in S_n,$$

$$P_n(x,y)=\frac{\kappa (x,y)}{d_n(x)}\forall x,y\in S_n.$$

$${\pi }_n(x)=\frac{d_n(x)}{{\sum }_{i=1}^n{d}_n(X_i)}\forall x\in S_n.$$

$${\mathcal{D}}_{n,t}(x,y)=\Vert P_n^t(x,\cdot )-P_n^t(y,\cdot ){\Vert }_{L^2(\text{mbox}diag(D_n^{-1}))}={\left\{\sum _{i=1}^n[P_n^t(x,X_i)-P_n^t(y,X_i){]}^2\frac{1}{d_n(X_i)}\right\}}^{1/2},$$

$${\mathcal{D}}_{n,t}(x,y)$$

$$\displaystyle\Psi_{n,t}:\,S_{n}\to\mathbb{R}^{n},\quad x\mapsto\big{(}\lambda_{n,0}^{t}\psi_{n,0}(x),\lambda_{n,1}^{t}\psi_{n,1}(x),\ldots,\lambda_{n,n-1}^{t}\psi_{n,n-1}(x)\big{)}^{T}.$$

$${\Delta }_{\mathcal{M}}f=-\text{mbox}div({\nabla }_{\mathcal{M}}f),$$

$${\int }_{\mathcal{M}}g(x){\Delta }_{\mathcal{M}}f\mathrm{d}\text{mbox}{Vol}_{\mathcal{M}}(x)={\int }_{\mathcal{M}}⟨{\nabla }_{\mathcal{M}}g(x),{\nabla }_{\mathcal{M}}f(x)⟩\mathrm{d}\text{mbox}{Vol}_{\mathcal{M}}(x),$$

$$\Delta f=-\sum _{j=1}^q\frac{{\partial }^{2}f}{\partial x_j^2}.$$

$${\Delta }_{\mathcal{M}}f=-\frac{1}{\sqrt{\det (\mathbf{G})}}\sum _{j=1}^q\frac{\partial }{\partial e^j}\left(\sqrt{\det (\mathbf{G})}\sum _{i=1}^q{g}^{ij}\frac{\partial f}{\partial e^i}\right),$$

$${\Delta }_{\mathcal{M}}f:={\Delta }_{\mathcal{M},\rho }f=-\frac{1}{{\rho }^2}\text{mbox}div({\rho }^2\nabla f)={\Delta }_{\mathcal{M}}f-2⟨\nabla \log \rho ,\nabla f⟩.$$

Full text

X. Chen’s research is supported in part by NSF DMS-1404891, NSF CAREER Award DMS-1752614, UIUC Research Board Awards (RB17092, RB18099), and a Simons Fellowship (Award Number: 663673). Y. Yang’s research is supported in part by NSF DMS-1810831. This work is completed in part with the high-performance computing resource provided by the Illinois Campus Cluster Program at UIUC. Authors are listed in alphabetical order

1. Introduction

This article studies the clustering problem of partitioning $n$ data points to $K$ disjoint (smooth) Riemannian submanifolds with $1\leqslant K\leqslant n$.

1.1. Problem formulation

Let $\mathcal{M}_{k},k=1,\dots,K$ be compact and connected Riemannian manifolds of dimension $q_{k}$. Suppose that $\mathcal{M}_{k}$ can be embedded as a submanifold of an ambient Euclidean space $\mathbb{R}^{p}$ equipped with the Euclidean metric $\|\cdot\|$ (i.e., there is an immersion $\varphi_{k}:\mathcal{M}_{k}\to\mathbb{R}^{p}$ such that the differential $\mathrm{d}\varphi_{k}$ is injective and $\varphi_{k}$ is a homeomorphism onto $\varphi_{k}(\mathcal{M}_{k})\subset\mathbb{R}^{p}$; cf. [10]). In our clustering setting, we work with disjoint submanifolds $\mathcal{M}_{1},\dots,\mathcal{M}_{K}$ in $\mathbb{R}^{p}$ and denote $S=\bigsqcup_{k=1}^{K}\mathcal{M}_{k}$ as their disjoint union. Each smooth submanifold $\mathcal{M}_{k}$ is endowed with the Riemannian metric $\rho_{k}$ induced from $\|\cdot\|$, and we denote $\mathscr{D}_{k}$ as the Borel $\sigma$-algebra on $\mathcal{M}_{k}$ (i.e., the $\sigma$-algebra generated by the open balls in $\mathcal{M}_{k}$ with respect to $\rho_{k}$). Let $X_{1}^{n}:=\{X_{1},\dots,X_{n}\}$ be a sequence of independent random variables taking values in $S$. Suppose that there exists a clustering structure $G_{1}^{*},\dots,G_{K}^{*}$ (i.e., a partition on $[n]:=\{1,\dots,n\}$ satisfying $\bigsqcup_{k=1}^{K}G_{k}^{*}=[n]$) such that each of the $n$ data points belongs to one of the $K$ clusters: if $i\in G_{k}^{*}$, then $X_{i}\sim\mu_{k}$ for some probability distribution $\mu_{k}$ supported on $\mathcal{M}_{k}$. Given the observations $X_{1}^{n}$, the task of this paper is to develop computationally tractable algorithms with strong theoretical guarantees for recovering the true clustering structure $G_{1}^{*},\dots,G_{K}^{*}$.

Classical clustering methods such as $K$-means [21] and mixture models [13] assume that data points from each cluster are sampled in the neighborhood (with the same dimension) of a centroid, where $\mathcal{M}_{k}$ contains only one point in $\mathbb{R}^{p}$. Such methods are effective for partitioning data with ellipsoidal contours, which implicitly implies that the similarity (or affinity) criteria of centroid-based clustering methods target on some notions of “compactness”. In modern applications such as image processing and computer vision [31, 33, 11], structured data with geometric features are commonly seen as clusters without necessarily being close together and having the same dimension. Figure 1 is an illustration for such observation on a synthetic data sampled from a noisy version of three clusters with one disk and two annuli. In this example, it is visually clear to distinguish the three clusters, however the $K$-means method fails to correctly cluster the data points. There are two main reasons for the failure of $K$-means. First, the north pole and south pole in the outer annulus have the largest Euclidean distance among all data points, even though they belong to the same cluster. Second, the annuli and the disk live in different dimensions. In particular, the annulus is a one-dimensional circle in $\mathbb{R}^{2}$ that is locally isometric to the real line and the disk has dimension two. Thus these geometric concerns motivate us to seek a more natural and flexible notion of closeness for clustering analysis. In this paper, we shall focus on the clustering criterion based on the connectedness, which is suitable for simultaneously addressing the two issues. First, connectedness is a graph property that does not rely on the physical distance: two vertices are connected if there is a path joining them. This extends the closeness from the local neighborhood to the global sense. Second, connectivity is a viable notion for clustering components of mixed dimensions, as long as all clusters live in the same ambient space where the graph connectivity weights can be computed.

In the population version, a clustering component can be viewed as a smooth submanifold, embedded in $\mathbb{R}^{p}$. In Riemannian geometry, a Riemannian submanifold $\mathcal{M}$ in $\mathbb{R}^{p}$ is said to be connected if for any $x,y\in\mathcal{M}$, there is a parameterized regular curve joining $x$ and $y$. Thus an appealing notion of $\mathcal{M}$ for being a cluster is that $\mathcal{M}$ is a compact and connected component in $\mathbb{R}^{p}$. In our setting, a clustering model is the union of $K$ disjoint submanifolds $\mathcal{M}_{1},\dots,\mathcal{M}_{K}$, and each $\mathcal{M}_{k}$ is equipped with a probability distribution $\mu_{k}$. Thus for the underlying true clustering model, there are $K$ connected graphs that do not overlap. In the sample version, data points are randomly generated from $\{(\mathcal{M}_{k},\mu_{k})\}_{i=1}^{k}$ with a clustering structure $\{G_{k}^{*}\}_{k=1}^{K}$. Typically, weighted graphs computed from the observed data points are fully connected (e.g., based on the Gaussian kernel). Thus a fundamental challenge of clustering analysis is to recover the true clustering structure from a noisy and fully connected weighted graph on the data.

1.2. Our contributions

In this paper, we propose a new clustering method, termed as the diffusion $K$-means, for manifold clustering. The diffusion $K$-means contains two key ingredients. First, it constructs a random walk (i.e., a Markov chain) on the weighted random graph with data points $X_{1}^{n}$ as vertices and edge weights computed from a kernel representing the similarity of data in a local neighborhood. By running the Markov chain forward in time, the local geometry (specified by the kernel bandwidth) will be integrated at multiple time scales to reveal global (topological) structures such as the connectedness properties of the graph. In the limiting case as the sample size tends to infinity and the bandwidth tends to zero, the random walk becomes a diffusion process over the manifold. By looking at the spectral decomposition of this limiting diffusion process, the evaluations of the eigenfunctions at vertices $X_{1}^{n}$ can be viewed as a continuous embedding of the data, called diffusion map, into a higher-dimensional Euclidean space. Second, once the diffusion map is obtained, we can compute the diffusion distance [9] and the $K$-means algorithm (with the Euclidean metric) can be naturally extended with the diffusion affinity as the similarity measure. Since the diffusion distance/affinity captures the connectedness among vertices on the weighted random graph, the diffusion $K$-means aims to maximize the within-cluster connectedness, which can be recast as an assignment problem via a 0-1 integer program.

Because 0-1 integer programming problems with a non-linear objective function is generally $\mathsf{NP}$-hard, solving the diffusion $K$-means is computationally intractable, i.e., polynomial-time algorithms with exact solutions only exists in special cases. This motivates us to consider semidefinite programming (SDP) relaxations. We propose two versions of SDP relaxations of the diffusion $K$-means. The first one requires the knowledge of the number of clusters, and it can be viewed as an extension from Peng and Wei’s SDP relaxation [25] for the $K$-means (as well as Chen and Yang’s SDP relaxation [7] for the generalized $K$-means for non-Euclidean data in an inner product space) to the manifold clustering setting with diffusion distances. Figure 1 (on the right) shows that the SDP relaxed diffusion $K$-means can correctly identify the three clusters in the previous example. The second SDP relaxation does not require the number of clusters as an input. The idea is to drop the constraint on the trace of the clustering membership matrix (which involves number of clusters $K$), and to add a penalization term on the diffusion $K$-means objective function. Thus it can be seen as a nuclear norm regularized version of the SDP for diffusion $K$-means that is adaptive to the number of clusters. For both SDP relaxations of the diffusion $K$-means, we show that exact recovery can be achieved when the underlying submanifolds are well separated and subsamples within each submanifold are well connected.

Since the diffusion $K$-means and its regularized version have only one (non-adaptive) bandwidth parameter to control the local geometry, they may fail for clustering problems with unbalanced sizes, mixed dimensions, and different densities. In such situations, a random walk on the vertices sampled from regions of low density mixes slower than that from regions of high density. This motivates us to consider a variant of diffusion $K$-means, termed as the localized diffusion $K$-means, by using data-dependent local bandwidth. We adopt the self-tuning procedure from [42] where local adaptive bandwidth is estimated from the nearest neighbors and we show that the localized diffusion $K$-means is adaptive to the local geometry and the local sampling density for the purpose of exact recovery of the true clustering structure.

To summarize, our contributions are listed as below.

1. (1)

We introduce the diffusion $K$-means clustering method for manifold clustering, which integrates the nonlinear embedding via the diffusion maps and the $K$-means clustering. 2. (2)

We propose two versions of the SDP relaxations of the diffusion $K$-means: one requires to know the number of clusters, and the other one does not require such knowledge as an input (and thus it is adaptive to the unknown number of clusters). 3. (3)

We derive the exact recovery property of the SDP relaxed diffusion $K$-means in terms of two hardness parameters of the clustering problem: one reflects the separation of the submanifolds, and the other one quantifies the degree of connectedness of the submanifolds. 4. (4)

We combine the local scaling procedures with the diffusion $K$-means and its regularized version, and derive their adaptivity when the clustering problems have unbalanced sizes, mixed dimensions, and different densities.

1.3. Related work

There is a large collection of clustering methods and algorithms in literature, which can be broadly classified into two categories: hierarchical clustering and partition-based clustering. Hierarchical clustering recursively divides data points into groups in either a top-down or bottom-up way. Such algorithms are greedy and they often get stuck into local optimal solutions.

Partition-based clustering methods such as $K$-means clustering [21] and spectral clustering [38] directly assign each data point with a group membership. Perhaps one of the most widely used clustering methods is the $K$-means method, due to the existence of algorithms with linear sample complexity (such as Lloyd’s algorithm [19]). However, the $K$-means clustering converges locally to a stationary point that depends on the initial partition. Recent theoretical studies in [20] show that, given a proper initialization (such as spectral clustering), Lloyd’s algorithm for optimizing the $K$-means objective function can consistently recover the clustering structures. Exact and approximate recovery of various convex relaxations for the $K$-means and mixture models are studied in literature [25, 16, 12, 29, 3]. To the best of our knowledge, existing theoretical guarantees developed for the convex relaxed $K$-means clustering assumes that the clusters are sampled in a neighborhood of a centroid. Thus results derived for $K$-means in literature cannot be directly compared with our results.

On the other hand, spectral clustering methods [31, 24] take the similarity matrix as the input and solve the clustering problem by applying $K$-means to top eigenvectors of the graph Laplacian matrix or its normalized versions [8]. In essence, spectral clustering contains two steps: (i) the Laplacian eigenmaps embed data into feature spaces, and (ii) $K$-means on top eigenvectors serves as a rounding procedure to obtain the true clustering structure [38]. Conventional intuition for spectral clustering is that the embedding step (i) often “magnifies” the cluster structure from the dataset to the feature space such that it can be revealed by a relatively simple algorithm (such as $K$-means) in step (ii). However, theoretical guarantees (such as exact recovery) for the spectral clustering is rather vague in literature, partially due to its two-step nature which complicates its theoretical analysis. For instance, [39] study the convergence of spectral properties of random graph Laplacian matrices constructed from sample points and they establish the consistency of the spectral clustering in terms of eigenvectors. However, they do not address the problem of the exact recovery property of the clustering structure. Similar results along this direction can be found in [27, 30, 35]. [17] propose similar SDP relaxations for the spectral clustering as in the present paper with the diffusion distance metric replaced with the graph Laplacian. Specifically, it is shown in [17] that those SDP relaxations can exactly recover the true clustering structure under a spectral proximity condition. Such condition is deterministic and difficult to check for general data generation models. (A particular checkable random model is the stochastic ball model [17].) During the preparation of this work, we notice a recent work [22] which proposes a similar idea of applying the diffusion distance as the similarity metric for clustering based on fast search and find of density peaks clustering (FSFDPS) [26]. To prove exact recovery, [22] requires strong deterministic assumptions on the Markov transition matrix associated with the diffusion process that could be difficult to check under their stochastic clustering model.

Literature on theoretical guarantees for manifold clustering is rather scarce, with a few exceptions [1, 18]. Near-optimal recovery of some emblematic clustering methods based on pairwise distances of data is derived under a condition that the minimal signal separation strength over all pairs of submanifolds is larger than a threshold. Compared with our diffusion $K$-means with local scaling, results established in [1] are non-adaptive to the local density and (geometric) structures of the submanifolds (cf. Theorem 4.10 and 4.12 ahead). [18] derive recovery guarantees for manifold clusters using a data-dependent metric called the longest-leg path distance (LLPD) that adapts to the geometry of data, where the data points are drawn from a mixture of uniform distributions on disjoint low-dimensional geometric objects.

1.4. Notation

For a matrix $A\in\mathbb{R}^{n\times n}$ and index subsets $G,G^{\prime}\subset[n]$, we use notation $A_{GG^{\prime}}$ to denote the submatrix of $A$ with rows being selected by $G$ and columns by $G^{\prime}$, and $\mbox{diag}(A)$ the $n$-dimensional vector composed of all diagonal entries of $A$. Let $\|A\|_{\infty}=\max_{1\leqslant i,j\leqslant n}|A_{ij}|$ and $\|A\|_{1}=\sum_{i,j=1}^{n}|A_{ij}|$ denote the $\ell_{\infty}$ and the $\ell_{1}$ norm of the vectorization $\mbox{vec}(A)$ of matrix $A$. Let ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|A\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\ast}$ and ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|A\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mbox{\scriptsize op}}$ denote the nuclear norm and the operator norm of matrix $A$. For two matrices $A$ and $B$ in $\mathbb{R}^{n\times n}$, we use $\langle A,B\rangle=\operatorname{tr}(A^{T}B)$ to denote a matrix inner product. We shall use $c,c^{\prime},c_{1},c_{2},\dots,C,C^{\prime},C_{1},C_{2},\dots$ to denote positive and finite (non-random) constants whose values may depend on the submanifolds $\{\mathcal{M}_{k}\}_{k=1}^{K}$ and the probability distributions $\{\mu_{k}\}_{k=1}^{K}$ supported on $\{\mathcal{M}_{k}\}_{k=1}^{K}$ and whose values may vary from place to place.

The rest of the paper is organized as follows. In Section 2, we discuss some related background on diffusion distances and nonlinear embeddings in Euclidean spaces, as well as the Laplace-Beltrami operator for the heat diffusion process on Riemannian manifolds. In Section 3, we introduce the diffusion $K$-means and its SDP relaxations. Regularized and localized diffusion $K$-means clustering methods are also proposed in this section. In Section 4, we present our main results on the exact recovery property of the SDP relaxed diffusion $K$-means. Simulation studies are presented in Section 5. Proofs are relegated to Section 6.

2. Preliminaries

Let $S\subset\mathbb{R}^{p}$ and $\mu$ be a probability distribution on $S$. Let $S_{n}=\{X_{1},X_{2},\ldots,X_{n}\}$ be $n$ i.i.d. random variables in $S$ sampled from $\mu$, and $\mu_{n}=n^{-1}\sum_{i=1}^{n}\delta_{X_{i}}$ the empirical distribution. In Section 2, we discuss the Euclidean embedding via the diffusion distance in the general setting. Thus in this section, we do not assume $S$ to be a disjoint union of $K$ submanifolds in $\mathbb{R}^{p}$ and the sample $S_{n}$ does not necessarily have a clustering structure.

2.1. Euclidean embedding and diffusion distances

Let $\kappa:S\times S\to\mathbb{R}$ be a positive semidefinite kernel that satisfies:

1. (i)

symmetry: $\kappa(x,y)=\kappa(y,x)$, 2. (ii)

positivity preserving: $\kappa(x,y)\geqslant 0$.

A kernel is a similarity measure between points of $S$. A widely used example is the Gaussian kernel:

[TABLE]

where $h>0$ is the bandwidth parameter that captures the local similarity of points in $S$. Given a kernel $\kappa$ with property (i) and (ii), we can define a reversible Markov chain on $S$ via the normalized graph Laplacian constructed as follows. Specifically, for any $x\in S$, let

[TABLE]

be the degree function of the graph on $S$. For simplicity, we assume $d(x)>0$ for all $x\in S$. Define

[TABLE]

which satisfies the positivity preserving property (ii) and the conservation property

[TABLE]

Thus $p(x,y)$ can be viewed as the one-step transition probability of a (stationary) Markov chain on $S$ from $x$ to $y$. We shall write this Markov chain (i.e., random walk) as $\mathcal{W}=(S,\mu,p)$, where $p(\cdot,\cdot)$ is called the transition kernel of $\mathcal{W}$. Equivalently, we can describe $\mathcal{W}$ by the bounded linear operator $\mathscr{P}:L^{2}(\mathrm{d}\mu)\to L^{2}(\mathrm{d}\mu)$ defined as

[TABLE]

Here $L^{2}(\mathrm{d}\mu):=L^{2}(S,\mathrm{d}\mu)$ is the class of squared integrable functions on $S$ with respect to $\mu$. In literature, $\mathscr{P}$ is often called the diffusion operator for the following reason. If we denote $p_{t}(x,y)$ as the $t$-step transition probability of the Markov chain $\mathcal{W}$ from $x$ to $y$ in $S$, then

[TABLE]

form a semi-group of bounded linear operators on $L^{2}(\mathrm{d}\mu)$. Let $\Pi$ be a stationary distribution of the Markov chain $\mathcal{W}$ over $S$. Then $\Pi$ is absolutely continuous with respect to $\mu$, and the probability density function $\pi$ of $\Pi$ with respect to $\mu$ is given by the Radon-Nikodym derivative

[TABLE]

Since $\Pi$ is the stationary measure of the Markov chain $\mathcal{W}$ with transition $P$, we have

[TABLE]

for all bounded measurable functions $f$, where $\Pi(f):=\int_{S}f(x)\,\mathrm{d}\Pi(x)$. Note that, since the kernel $\kappa$ is symmetric, $\mathcal{W}$ is reversible and satisfies the detailed balance condition:

[TABLE]

Lemma 2.1* (Spectral decomposition of the Markov chain $\mathcal{W}$).*

Let

[TABLE]

If

[TABLE]

then the following statements hold.

1. (1)

There exists a sequence of nonnegative eigenvalues $\lambda_{0}\geqslant\lambda_{1}\geqslant\cdots\geqslant 0$ such that

[TABLE]

where $\{\phi_{j}\}_{j=0}^{\infty}$ is the set of associated eigenfunctions to $\{\lambda_{j}\}_{j=0}^{\infty}$, and $\{\phi_{j}\}_{j=0}^{\infty}$ forms an orthonormal basis of $L^{2}(\mathrm{d}\mu)$. In particular, $\lambda_{0}=1$ and $\phi_{0}(x)=\sqrt{\pi(x)}$. 2. (2)

The transition probability $p(x,y)$ admits the following decomposition

[TABLE]

where $\psi_{j}(x)=\phi_{j}(x)/\sqrt{d(x)}$ and $\widetilde{\psi}_{j}(y)=\phi_{j}(y)\sqrt{d(y)}$. 3. (3)

The diffusion operator $\mathscr{P}$ satisfies

[TABLE]

In addition, $\lambda_{0}=1$ and $\psi_{0}\equiv(\int_{S}d(y)\,\mathrm{d}\mu(y))^{-1/2}$.

The proof of Lemma 2.1 is given in Appendix A, and our argument is similar to Lemma 12.2 in [15] in the finite-dimensional setting. If $\mathcal{W}$ is irreducible (i.e., the graph on $S$ is connected in that for all $x,y\in S$, there is some $t>0$ such that $p_{t}(x,y)>0$), then the stationary distribution $\pi$ is unique. Thus if we run this Markov chain $\mathcal{W}$ forward in time, then the local geometry (captured by the kernel $\kappa$ which is parameterized by the bandwidth $h$) will be integrated to reveal global structures of $S$ at multiple (time) scales. In particular, we can define a class $\{\mathscr{D}_{t}\}_{t\in\mathbb{N}_{+}}$ of diffusion distances [9] on $S$ by

[TABLE]

Roughly speaking, for each $t\in\mathbb{N}_{+}$ and $x,y\in S$, the diffusion distance $\mathscr{D}_{t}(x,y)$ quantifies the total number of paths with length $t$ connecting $x$ and $y$ (see Figure 2), thereby reflecting the local connectivity at the time scale $t$.

Lemma 2.2* (Spectral representation of diffusion distances).*

If the Markov chain $\mathcal{W}=(S,\mu,p)$ is irreducible, then we have

[TABLE]

for all $t\in\mathbb{N}_{+}$ and $x,y\in S$.

The proof of Lemma 2.2 is given in Appendix A. For an irreducible Markov chain, the spectral gap is strictly positive (i.e., $|\lambda_{j}|<1$ for all $j>0$). Based on the spectral decomposition in Lemma 2.2 and noting that $\psi_{0}\equiv(\int_{S}d(y)\,\mathrm{d}\mu(y))^{-1/2}$ is constant, we see that the diffusion distance can be written as

[TABLE]

In this case, the diffusion distance $\mathscr{D}_{t}(x,y)$ decays to zero as $t$ increases, provided that $x$ and $y$ belong to a connected component of the graph on $S$. In particular, the decay rate of the spectrum quantifies the connectivity of points in the graph on $S$. Given a positive integer $q\in\mathbb{N}_{+}$, the diffusion maps $\{\Psi_{t}\}_{t\in\mathbb{N}}$ are defined as

[TABLE]

where the $\ell$-th component $\Psi_{t\ell}^{(q)}(x)$ is the $\ell$-th diffusion coordinate in $\mathbb{R}^{q}$. Thus we obtain an embedding of $(S,\mu)$ into the Euclidean space $\mathbb{R}^{q}$ in the limiting sense that

[TABLE]

2.2. Empirical diffusion embedding

Recall that $S_{n}=\{X_{1},X_{2},\ldots,X_{n}\}$ are $n$ i.i.d. random variables in $S$ sampled from $\mu$, and $\mu_{n}=n^{-1}\sum_{i=1}^{n}\delta_{X_{i}}$ is the empirical distribution. Given $S_{n}$, we can consider finite sample approximations $\{\mathscr{D}_{n,t}\}_{t\in\mathbb{N}_{+}}$ to the underlying population level quantities $\{\mathscr{D}_{t}\}_{t\in\mathbb{N}_{+}}$. More precisely, consider a weighted graph with nodes corresponding to the elements in $S_{n}$, where the weight between a pair $(X_{i},X_{j})$ of nodes is $\kappa(X_{i},X_{j})$, for $i,j\in[n]$. Define the (rescaled) empirical degree function $d_{n}:\,S_{n}\to\mathbb{R}_{+}$ by

[TABLE]

where we added an extra $n$-factor so that $d_{n}(X_{i})$ is also the degree of node $X_{i}$ in the weighted graph. Let $D_{n}$ denote the $n$-by-$n$ diagonal matrix whose $i$-th diagonal entry is $d_{n}(X_{i})$. Consider the (empirical) random walk $\mathcal{W}_{n}=(S_{n},\mu_{n},P_{n})$ over $S_{n}$ with transition probability

[TABLE]

The (empirical) stationary distribution $\pi_{n}$ of the random walk $\mathcal{W}_{n}$ over $S_{n}$ becomes

[TABLE]

For any vector $v\in\mathbb{R}^{n}_{+}$, let $L^{2}(v)=\{u=(u_{1},\ldots,u_{n})\in\mathbb{R}^{n}:\,\|u\|_{L^{2}(v)}=\sum_{i=1}^{n}v_{i}\,u_{i}^{2}\}$ denote a weighted $L^{2}$ space over $S_{n}$. We define the empirical diffusion distances $\{\mathscr{D}_{n,t}\}_{t\in\mathbb{N}_{+}}$ as

[TABLE]

for all $x,y\in S_{n}$ and $t\in\mathbb{N}_{+}$. Roughly speaking, $\mathscr{D}_{n,t}$ provides an empirical estimate to $\mathscr{D}_{t}$. Similar to the spectral representation (5) for $\mathscr{D}_{t}$, we also have the following spectral representation of $\mathscr{D}_{n,t}$ (see Appendix B),

[TABLE]

where $1=\lambda_{n,0}\geqslant\lambda_{n,1}\geqslant\cdots\geqslant\lambda_{n,n-1}\geqslant 0$ are the nonnegative eigenvalues (due to the positive semidefiniteness of the kernel $\kappa$) of the transition probability operator $P_{n}$, which can be identified with a matrix in $\mathbb{R}^{n\times n}$ with $[P_{n}]_{ij}=P_{n}(X_{i},X_{j})$ as its $(i,j)$-th element, and $\psi_{n,0},\psi_{n,1},\ldots,\psi_{n,n-1}:\,S_{n}\to\mathbb{R}$ are the associated eigen-functions on $S_{n}$ with unit $L^{2}(\mbox{diag}(D_{n}))$ norm, which can be identified with vectors in $\mathbb{R}^{n}$ with $[\psi_{n,j}]_{i}=\psi_{n,j}(X_{i})$ as the $i$-th element of $\psi_{n,j}$ for $i\in[n]$ and $j=0,1\ldots,n-1$. The empirical diffusion distance $\mathscr{D}_{n,t}(X_{i},X_{j})$ between two nodes $X_{i}$ and $X_{j}$ is also the Euclidean distance between their embeddings $\Psi_{n,t}(X_{i})$ and $\Psi_{n,t}(X_{j})$ via the empirical diffusion map

[TABLE]

2.3. The Laplace-Beltrami operator on Riemannian manifolds

The Laplace-Beltrami operator on Riemannian manifolds is a generalization of the Laplace operator on Euclidean spaces. Let $f:\mathcal{M}\to\mathbb{R}$ be an (infinitely) differentiable function with continuous derivatives on a $q$-dimensional compact and smooth Riemannian manifold and $\nabla_{\mathcal{M}}f$ be the gradient vector field on $\mathcal{M}$ (i.e., $\nabla_{\mathcal{M}}f(x)$ is the deepest direction of ascent for $f$ at the point $x\in\mathcal{M}$). The Laplace-Beltrami operator $\Delta_{\mathcal{M}}$ is defined as the divergence of the gradient vector

[TABLE]

where the div operator is relative to the volume form $\mbox{Vol}_{\mathcal{M}}$ of $\mathcal{M}$. Here we adopt the convention with the minus sign of the divergence such that $\Delta_{\mathcal{M}}$ is a positive-definite operator. With integration-by-parts, we have for any two differentiable functions $f$ and $g$,

[TABLE]

where the inner product is taken in the $q$-dimensional tangent space of $\mathcal{M}$ (at the point $x$). In a Euclidean space (i.e., $\mathcal{M}=\mathbb{R}^{q}$), the Laplace-Beltrami operator is the usual Laplace operator

[TABLE]

On a general $q$-dimensional Riemannian manifold $\mathcal{M}$, the Laplace-Beltrami operator in a local coordinate system $(e^{1},\dots,e^{q})$ with a metric tensor $\mathbf{G}=(g_{ij})_{i,j=1}^{q}$ is given by

[TABLE]

where $g^{ij}$ are the entries of $\mathbf{G}^{-1}$. In the special case $\mathcal{M}=\mathbb{R}^{q}$, $\mathbf{G}$ is the $q\times q$ identity matrix.

More generally, if there is a distribution with a positive and Lipschitz continuous density $\rho>0$ on $\mathcal{M}$, then the drift Laplace-Beltrami operator is defined as

[TABLE]

Note that $\Delta_{\mathcal{M}}$ is a self-adjoint positive-definite compact operator, its spectrum contains a sequence of nonnegative eigenvalues $0\leqslant\lambda_{0}\leqslant\lambda_{1}\leqslant\cdots$. If in addition $\mathcal{M}$ is connected, then the second smallest eigenvalue $\lambda_{1}>0$. As we will show, $\lambda_{1}$ depends on the connectivity of the manifold (Figure 2), thus characterizing the limiting mixing time of the empirical random walk $\mathcal{W}_{n}$ over the $S_{n}$ as $n\to\infty$ and $h\to 0^{+}$, when $S_{n}$ is sampled from the manifold $\mathcal{M}$.

3. Diffusion $K$-means

Recall that in our clustering model, $S_{n}=\{X_{1},X_{2},\ldots,X_{n}\}$ is a sample of independent random variables taking values in $S$, where $S$ is the union of $K$ disjoint Riemannian submanifolds $\mathcal{M}_{1},\dots,\mathcal{M}_{K}$ embedded in the ambient space $\mathbb{R}^{p}$. The clustering problem is to divide these $n$ data points into $K$ clusters, so that points in the same cluster belongs to the same connected component in $S$, based on certain similarity measures between the points. In particular, the (classical) $K$-means clustering method minimizes the total intra-cluster squared Euclidean distances in $\mathbb{R}^{p}$

[TABLE]

over all possible partitions on $[n]$, where $|G_{k}|$ is the cardinality of $G_{k}$. Dropping the sum of squared norms $\sum_{i=1}^{n}\|X_{i}\|^{2}$, we see that the $K$-means clustering is equivalent to the maximization of the total within-cluster covariances

[TABLE]

Here, $a_{ij}=X_{i}^{T}X_{j}$ can be viewed as a similarity measure specified by the Euclidean space inner product $\langle X_{i},X_{j}\rangle_{\mathbb{R}^{p}}$. In general, we can replace the Euclidean inner product with any other inner product over $S_{n}$ [7]. For manifold clustering, we replace it with the inner product induced from the empirical diffusion distance, that is,

[TABLE]

Henceforth, we will refer to $\langle\cdot,\,\cdot\rangle_{\mathscr{D}_{n,t}}$ as the diffusion affinity. Interestingly, we can obtain this diffusion affinity value without explicitly conducting eigen-decomposition (spectral decomposition) to the transition probability matrix $P_{n}=D_{n}^{-1}\mathcal{K}_{n}$ (or the symmetrized matrix $D_{n}^{-1/2}\mathcal{K}_{n}D_{n}^{-1/2}$), where recall that D_{n}=\mbox{diag}\big{(}d_{n}(X_{1}),\ldots,d_{n}(X_{n})\big{)}\in\mathbb{R}^{n} is the degree diagonal matrix, and \mathcal{K}_{n}=\big{[}\kappa(X_{i},X_{j})\big{]}_{n\times n}\in\mathbb{R}^{n\times n} is the empirical kernel matrix. In fact, we may use the following relation that links the empirical diffusion affinity with entries in matrix $P_{n}$ raising to power $2t$ (see Appendix B for details),

[TABLE]

This motivates a $K$-means clustering method via diffusion distances, referred to as the diffusion $K$-means as

[TABLE]

for the tuning parameter $t$ interpreted as the number of steps in the empirical random walk $\mathcal{W}_{n}$. Note that here the affinity matrix $P_{n}^{2t}D_{n}^{-1}=D_{n}^{-1/2}(D_{n}^{-1/2}\mathcal{K}_{n}D_{n}^{-1/2})^{2t}D_{n}^{-1/2}\in\mathbb{R}^{n\times n}$ is symmetric. In light of the connections between the diffusion distance and the random walk $\mathcal{W}_{n}$ over $S_{n}$ in Section 2.2, the diffusion $K$-means attempts to maximize the total within-cluster connectedness.

Remark 3.1* (Intuition of diffusion $K$-means).*

In Section 2.1, we see that, on a connected submanifold $\mathcal{M}_{k}$, the (population) diffusion process (2) converges to the stationary distribution (3):

[TABLE]

In fact, since the kernel $\kappa$ is positive semidefinite, this convergence holds at a geometric rate governed by the spectral gap of the Laplace-Beltrami operator on $\mathcal{M}_{k}$ (cf. (43) in the proof of Lemma 6.6). Thus the empirical version of the diffusion (i.e., the Markov chain on the random graph generated by $X_{1}^{n}$ and $\kappa$) obeys

[TABLE]

for any two data points $X_{i},X_{j}\in\mathcal{M}_{k}$. On the other hand, if the separation between the submanifolds is large enough and $t$ is not so large, then the probability to diffuse from one cluster to another one is small (cf. Lemma 6.3). Thus we expect that

[TABLE]

for any two data points $X_{i}\in\mathcal{M}_{k}$ and $X_{j}\in\mathcal{M}_{m}$ such that $k\neq m$. This means that the within-cluster entries of the empirical diffusion affinity are larger than the between-cluster entries. In particular, for suitably large $t\in\mathbb{N}_{+}$, the empirical diffusion affinity matrix $A_{n}:=A_{n,t}=P_{n}^{2t}D_{n}^{-1}$ tends to become close to a block-diagonal matrix

[TABLE]

where $N_{k}=\sum_{\ell,\ell^{\prime}\in G_{k}^{*}}\kappa(X_{\ell},X_{\ell^{\prime}})$. Since each diagonal block of $A_{n}$ tends to be a constant matrix, if we run the Markov chain for a suitably long time, then this block-diagonal structure in the limit precisely conveys the true clustering structure in that $i,j\in G_{k}^{*}$ if and only if $\lim_{t\to\infty}[A_{n,t}]_{i,j}=N_{k}^{-1}>0$. The trade-off regime of $t$ (cf. (17) in Theorem 4.1) is determined by the non-asymptotic bounds on the convergence of the empirical diffusion maps to its population version (cf. Lemma 6.2 and 6.3), as well as the submanifolds separation. $\blacksquare$

Note that, for every partition $G_{1},\dots,G_{K}$, there is a one-to-one $n\times K$ assignment matrix $H=(h_{ik})\in\{0,1\}^{n\times K}$ such that $h_{ij}=1$ if $i\in G_{k}$ and $h_{ij}=0$ if $i\notin G_{k}$. Thus the diffusion $K$-means clustering problem can be recast as a 0-1 integer program:

[TABLE]

where $\mathbf{1}_{n}$ denotes the $n\times 1$ vector of all ones and $B=\text{diag}(n_{1}^{-1},\dots,n_{K}^{-1})$, where $n_{k}=|G_{k}|$ for $k=1,\ldots,K$ is the size of the $k$-th cluster.

The diffusion $K$-means clustering problem (9) is often computationally intractable, namely, polynomial-time algorithms with exact solutions only exist in special cases [32]. For instances, the (classical) $K$-means clustering is an $\mathsf{NP}$-hard integer programming problem with a non-linear objective function [25]. Exact and approximate recovery of various SDP relaxations for the $K$-means [25, 16, 12, 29, 14] are studied in literature. However, it remains a challenging task to provide statistical guarantees for clustering methods that can capture non-linear features of data taking values on manifolds.

3.1. Semidefinite programming relaxations

We consider the SDP relaxations for the diffusion $K$-means clustering. Note that every partition $G_{1},\dots,G_{K}$ of $[n]$ can be represented by a partition function $\sigma:[n]\to[K]$ via $G_{k}=\sigma^{-1}(k),k=1,\dots,n$. If we change the variable $Z=HBH^{T}$ in the 0-1 integer program formulation (9) of the diffusion $K$-means, then $Z$ satisfies the following properties:

[TABLE]

For the diffusion $K$-means $B=\text{diag}(n_{1}^{-1},\dots,n_{K}^{-1})$, the last constraint in (10) reduces to $Z\mathbf{1}_{n}=\mathbf{1}_{n}$, which does not depend on the partition function $\sigma$. Thus we can relax the diffusion $K$-means clustering to the SDP problem:

[TABLE]

where $Z\succeq 0$ means that $Z$ is positive semidefinite and $Z\geqslant 0$ means that all entries of $Z$ are non-negative, and matrix $A=[a_{ij}]:=A_{n}=P_{n}^{2t}D_{n}^{-1}\in\mathbb{R}^{n\times n}$. We shall use $\hat{Z}$ to estimate the true “membership matrix” $Z^{*}$, where

[TABLE]

Note that $Z^{*}$ is a block diagonal matrix (up to a permutation) of rank $K$. If $X_{1},\dots,X_{n}\in\mathbb{H}$ (i.e., $\mathbb{S}=\mathbb{H}$) for some Hilbert space $\mathbb{H}$ and $a_{ij}=\langle X_{i},X_{j}\rangle_{\mathbb{H}}$ is the inner product between $X_{i}$ and $X_{j}$, then (11) is the SDP for kernel $K$-means proposed in [7]. In particular, [25] consider the special case for the (Euclidean) $K$-means, where $\mathbb{H}=\mathbb{R}^{p}$ and $a_{ij}=X_{i}^{T}X_{j}$. Observe that the SDP relaxation (11) does not require the knowledge of the cluster sizes other than the number of clusters $K$. Thus it can handle the general case for unequal cluster sizes.

3.2. Regularized diffusion $K$-means

In practice, the number $K$ of clusters is rarely known. Note that the SDP problem (11) depends on $K$ only through the constraint $\operatorname{tr}(Z)=K$. Therefore we propose a regularized diffusion $K$-means estimator by dropping the constraint on the trace and penalizing $\operatorname{tr}(Z)$ as follows,

[TABLE]

where $\rho>0$ is the regularization parameter. Recall that ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|Z\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\ast}$ denotes the nuclear norm of a matrix $Z$ (i.e., ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|Z\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\ast}$ is the sum of the singular values of $Z$). Since ${\left|\kern-1.07639pt\left|\kern-1.07639pt\left|Z\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\ast}=\operatorname{tr}(Z)$ for $Z\in\mathscr{C}$, it is interesting to note that (13) is the same as the nuclear norm penalized diffusion $K$-means

[TABLE]

Recall that the true membership matrix $Z^{*}$ has a block diagonal structure with rank $K$, the nuclear norm penalized diffusion $K$-means (14) can be thought as an $\ell^{1}$ norm convex relaxation of the $\mathsf{NP}$-hard rank minimization problem (i.e., minimizing the number of non-zero eigenvalues). Thus the parameter $K$ in the clustering problem plays a similar role as the sparsity (or low-rankness) parameter in the matrix completion context [5]. Hence the SDP problem (14) can be viewed as a (further) convex relaxation of the infeasible SDP problem (11) when $K$ is unknown. Note that similar regularizations have been considered in [3] for the $G$-latent clustering models and in [40] for stochastic block models.

It remains a question to choose the value of $\rho$. Larger values of $\rho$ will lead to solutions containing less number of clusters (with larger sizes). In particular, when matrix $A$ is positive-definite, the following Lemma 3.2 shows that the solution reduces to a rank one matrix that assigns all points into a giant cluster when $\rho$ is large enough, and becomes the identity matrix that assigns $n$ points into $n$ distinct clusters when $\rho$ is small enough. In addition, the trace $\operatorname{tr}(\tilde{Z}_{\rho})$ of the solution is nonincreasing in $\rho>0$.

Lemma 3.2* (Monotonicity of tuning path).*

Suppose $A$ is a positive definite matrix, and let $\lambda_{\max}(A)$ and $\lambda_{\min}(A)$ denote its respective largest and smallest eigenvalues.

(1) If $n\rho>\lambda_{\max}(A)$, then $Z^{\diamond}=n^{-1}J_{n}$, where $J_{n}$ is the $n\times n$ matrix of all ones, is the unique solution of the SDP (13).

(2) If $n\rho<\lambda_{\min}(A)$, then $Z^{\dagger}=I_{n}$, the $n\times n$ identity matrix, is the unique solution of the SDP (13).

(3) If $\tilde{Z}_{\rho_{1}}$ and $\tilde{Z}_{\rho_{2}}$ are two solutions of the SDP (13) with the regularization parameter taking values $\rho_{1}$ and $\rho_{2}$, respectively. If $\rho_{1}<\rho_{2}$, then $\operatorname{tr}(\tilde{Z}_{\rho_{1}})\geqslant\operatorname{tr}(\tilde{Z}_{\rho_{2}})$. Furthermore, if at least one of the two SDPs has a unique solution, then $\operatorname{tr}(\tilde{Z}_{\rho_{1}})>\operatorname{tr}(\tilde{Z}_{\rho_{2}})$.

According to the interpretation of the SDP (11), the trace $\operatorname{tr}(Z)$ of the solution can be viewed as the fitted number of clusters. Consequently, Lemma 3.2 implies that smaller values of $\rho$ will result in more clusters (with smaller sizes) in $\tilde{Z}_{\rho}$. In practice, we need to properly select the tuning parameter $\rho$. We propose the following decision rule for this purpose. For each $\rho$, we run the SDP problem (13) and extract the value of $\operatorname{tr}(\tilde{Z}_{\rho})$. Then we plot the solution path $\operatorname{tr}(\tilde{Z}_{\rho})$ versus $\log(\rho)$ and pick the values of $\rho$ which spend the longest time (on the logarithmic scale) with a flat value of $\operatorname{tr}(\tilde{Z}_{\rho})$. Here we recommend using the logarithmic scale since values of $\rho$ with non-trivial solutions $\operatorname{tr}(\tilde{Z}_{\rho})$ tends to be close to zero. Algorithm 1 below summarizes this decision rule.

Figure 3 shows the empirical result on the three clusters (one disk and two annuli) example in Section 1. According to Lemma 3.2, the estimated number of clusters, proxied by $\operatorname{tr}(\tilde{Z}_{\rho})$, is a non-increasing function of $\rho$. In particular, the trace $\operatorname{tr}(\tilde{Z}_{\rho})$ in the solution path in the upper left panel of Figure 3 stabilizes around $2$ and $3$, indicating that both $2$ and $3$ are candidate values for the number of clusters. In particular, the interval of $\rho$ (on the logarithmic scale) corresponding to value $3$ is much larger than that to value $2$, indicating that $3$ is more likely to be the true number of clusters (cf. the (correct) case in Figure 3). In Section 4, we will use our theory to explain this stabilization phenomenon, which partially justifies our $\rho$ selection rule.

In addition, as we can see from the rest three panels in Figure 3, by gradually increasing $\rho$, the adaptive diffusion $K$-means method produces a hierarchical clustering structure. Unlike the top-down or bottom-up clustering procedures which are based on certain greedy rule and can incur inconsistency, the hierarchical clustering structure produced by our approach is consistent — it does not depend on the order of partitioning or merging due to the uniqueness of the global solution from the convex optimization via the SDP. Similar observations can be drawn on another example shown in Figure 4 containing a uniform sample on three rectangles (see DGP 2 in our simulation studies Section 5 for details).

Further, it is interesting to observe in Figure 3 that the regularized diffusion $K$-means tuned with two clusters yields a merge between the outer annulus and the disk, which gives the largest total diffusion affinity in the objective function (7) among the three possible combinations of the true clusters. Since diffusion affinity decays exponentially fast to zero in the squared Euclidean distance (for the Gaussian kernel), the diffusion affinity matrix $A=P_{n}^{2t}D_{n}^{-1}$ tends to have a block diagonal structure, as weights between points belonging to different clusters are exponentially small (cf. (8) and Lemma 6.4). Thus running SDP for examples with relatively well separated clusters, such as the one in Figure 3, tends to merge two clusters with largest within-cluster diffusion affinities that is irrespective of the between-cluster Euclidean distances. This may lead to a visually less appealing merge as in the Euclidean distance case (cf. the (under) case in Figure 3). On the other hand, the regularized diffusion $K$-means is able to produce more reasonable partition in splitting the clusters (cf. the bottom-left panel in Figure 3). In particular, if the regularization parameter $\rho$ is chosen such that the corresponding number $\hat{K}$ of clusters in the SDP solution is greater than $K$, then this will cause a split in one of the true clustering structures that minimizes the between-cluster diffusion affinities after the splitting. Moreover, in our simulation studies (setup DGP 3 in Section 5), we observe that the SDP relaxed regularized diffusion $K$-means performs much better in harder cases than the spectral clustering methods when the true clusters are not well separated.

3.3. Localized diffusion $K$-means

For clustering problems with different sizes, dimensions and densities, the diffusion $K$-means may have limitations since only one bandwidth parameter $h$ is used to control the local geometry on the domain. More precisely, according to our theory (for example, Theorem 4.1 below), the optimal choice of the bandwidth parameter for a $q$-dimensional submanifold as one connected component in our clustering model is $h\asymp(\log\tilde{n}/\tilde{n})^{1/q}$, where $\tilde{n}$ corresponds to the sample size within this cluster and thus depends on the local cluster size or density level. Figure 5 demonstrates such an example for a mixture of three bivariate Gaussians, which consist of one larger Gaussian component with low density and two smaller Gaussian components with high density. Empirically, the diffusion $K$-means fails on this example (even after tuning) for the reason that the larger and smaller clusters have very different local densities. This motivates us to consider a variant of diffusion $K$-means, termed as localized diffusion $K$-means, by using local adaptive bandwidth $h_{i}=h(X_{i})$ for each $X_{i}$, $i\in[n]$. In particular, we adopt the self-tuning procedure from [42] by setting $h_{i}$ to be $\|X_{i}-X_{i}^{(k_{0})}\|$, where $X_{i}^{(k_{0})}$ denotes the $k_{0}$-th nearest neighbor to $X_{i}$, and replacing $\mathcal{K}_{n}$ with $\mathcal{K}_{n}^{\dagger}=[\kappa^{\dagger}(X_{i},X_{j})]_{n\times n}$ (and accordingly replacing $D_{n}$ with $D_{n}^{\dagger}$ corresponding to the diagonal degree matrix associated with $\mathcal{K}^{\dagger}_{n}$) given by

[TABLE]

in the SDP (11). Note that $\mathcal{K}^{\dagger}_{n}$ is generally no longer a positive semidefinite matrix. Intuitively for $i\in G_{k}^{*}$, the local scaling $h_{i}$ automatically adapts to the local density $p_{k}(X_{i})$ about $X_{i}$, the cluster size $n_{k}$ and the dimension $q_{k}$ for the $k$-th Riemannian submanifold. Specifically, for each cluster $k=1,\ldots,K$, the $n_{k}$-by-$n_{k}$ submatrix $[\mathcal{K}_{n}^{\dagger}]_{G^{\ast}_{k}G^{\ast}_{k}}$ resembles the Gaussian kernel matrix with a homogeneous bandwidth $h_{k}\asymp(\log n_{k}/n_{k})^{1/q_{k}}$ that adapts to the local geometry in $\mathcal{M}_{k}$. For points $X_{i}$ and $X_{j}$ belonging to distinct clusters that are properly separated, $\kappa^{\dagger}(X_{i},X_{j})$ tends to be close to zero and is less affected by the choice of $h_{i}$ and $h_{j}$. Heuristically, $h_{i}$ is larger for lower density regions where the degree function of $X_{i}$ is smaller so that the random walk can speed up mixing at such lower density regions. Overall, such a locally adaptive choice of bandwidth improves the mixing time of the random walk within each cluster, while leaves the between cluster jumping probabilities remaining small. As a consequence, the pairwise diffusion affinity matrix $A$ in our SDP formulation (11) tends to exhibit a clearer block form reflecting the clustering structure. To compute $h_{i}$, we only need to specify $k_{0}$ to replace the (non-adaptive) bandwidth parameter $h$ whose value depend on the unknown cluster sizes $n_{k}$, dimensions $q_{k}$ of submanifolds $\mathcal{M}_{k}$, and the underlying probability density functions $p_{k}$ on $\mathcal{M}_{k}$. In contrast, the simple choice $k_{0}=\lfloor C\log{n}\rfloor$ guarantees that the local scaling $h_{i}$ adapts to the local density (cf. Theorem 4.10 in Section 4).

4. Main results

In this section, we assume each $\mathcal{M}_{k}$ is a compact connected $q_{k}$-dimensional Riemannian submanifold embedded in $\mathbb{R}^{p}$ with bounded diameter ${\rm Diam_{k}}$, absolute sectional curvature value ${\rm Sec_{k}}$, injectivity radius ${\rm Inj_{k}}$, and inverse reach ${\rm Rch_{k}^{-1}}$, where \rm Rch_{k}:\,=\sup\big{\{}t>0:\,\forall x\in\mathbb{R}^{p},\,\mbox{dist}(x,\mathcal{M}_{k})\leqslant t,\,\exists!y\in\mathcal{M}_{k}\,\mbox{s.t.~{}dist}(x,\mathcal{M}_{k})=\|x-y\|\big{\}}. Throughout the rest of the paper, we assume that each $\mu_{k}$ has a ${\rm Lip_{k}}$-Lipschitz continuous density function $p_{k}$ with respect to the Riemannian volume measure on $\mathcal{M}_{k}$, such that

[TABLE]

for some constant $\eta_{k}>0$. For simplicity, we fix the Gaussian kernel in (1) for the diffusion $K$-means and its variants, which can be relaxed to more general isotropic kernel with exponential decay in the squared Euclidean distance (cf. Remark 4.7).

4.1. Exact recovery of diffusion $K$-means

Let $\delta=\min_{1\leqslant k\neq k^{\prime}\leqslant K}\|\mathcal{M}_{k}-\mathcal{M}_{k^{\prime}}\|$, where $\|A-B\|=\inf_{x\in A,y\in B}\|x-y\|$ denote the (Euclidean) distance between two disjoint sets $A$ and $B$. Recall that the size of the true cluster $n_{k}=|G_{k}^{*}|$ for $k=1,\ldots,K$ and let $\underline{n}=\min_{1\leq k\leq K}n_{k}$ denote the minimal cluster size.

Theorem 4.1* (Exact recovery of SDP for diffusion $K$-means).*

Fix any $\xi\geqslant 1$. Let $c_{1},c_{2},c_{3},M_{0}$ be some positive constants only depending on $\{q_{k},{\rm Sec_{k}},{\rm Inj_{k}},{\rm Rch_{k}},\eta_{k},{\rm Lip_{k}}\}_{k=1}^{K}$, and $C_{1},C_{2}$ only depending on $\{q_{k},{\rm Sec_{k}},{\rm Inj_{k}},{\rm Rch_{k}},\eta_{k},{\rm Lip_{k}}\}_{k=1}^{K}$ and $\xi$. If $\underline{n}\geqslant M_{0}$ and

[TABLE]

where $\mathbf{1}(q_{k}=2)$ equals to $1$ if $q_{k}=2$ and [math] otherwise, then we can achieve exact recovery, that is $\hat{Z}=Z^{\ast}$, with probability at least $1-c_{3}K\,\underline{n}^{-\xi}$ as long as

[TABLE]

where $\nu_{d}=\frac{\pi^{d/2}}{\Gamma(d/2+1)}$ denotes the volume of the $d$-dim unit ball, and $\lambda_{1}(\mathcal{M}_{k})>0$ denotes the second smallest eigenvalue of the drift Laplace-Beltrami operator (with the minus sign) on $\mathcal{M}_{k}$ defined in (6).

A proof of this theorem is provided in Section 6.1. A crucial step in the proof is bounding from below the absolute spectral gap of the transition matrix associated with the restricted random walk $\mathcal{W}_{n}$ onto each submanifold $\mathcal{M}_{k}$ for $k=1,\ldots,K$. We do so by applying a comparison theorem of Markov chains to connect this spectral gap with the eigensystem of the drift Laplace-Beltrami operator over the submanifold $\mathcal{M}_{k}$ that captures the diffusion geometry. In particular, our proof borrows existing results [4, 34] on error estimates of using the spectrum of a random geometric graph to approximate the eigensystem of the drift Laplace-Beltrami operator in the numerical analysis literature.

Remark 4.2* (Comments on the exact recovery condition (17)).*

We begin with comments on the second term on the left hand side of (17). We point out that this is essentially a mixing condition on random walks over the $K$ submanifolds. In particular, it is due to a relation between the mixing times of the heat diffusion process on each submanifold and its discretized random walk $\mathcal{W}_{n}$ over vertices sampled from the submanifold.

For simplicity, we illustrate this relation for $K=1$ and uniform density $p_{1}$ on $\mathcal{M}_{1}=\mathbb{S}^{1}$, where $\mathbb{S}^{1}$ is the unit circle in $\mathbb{R}^{2}$. As a one-dimensional compact smooth manifold (without boundary), $\mathbb{S}^{1}$ can be (globally) parametrized by the angle $\theta\in[-\pi,\pi)$. Under this parametrization, the density function $u(\tau,\theta)$ of the heat diffusion process on $\mathbb{S}^{1}$, as a function of time $\tau$ and location $\theta$, is determined by the corresponding heat equation,

[TABLE]

where $\Delta$ is the Laplace-Beltrami operator on $\mathbb{S}^{1}$. Under the same parametrization, the Laplace operator $\Delta=-{\mathrm{d}^{2}\over\mathrm{d}\theta^{2}}$ on $\mathbb{S}^{1}$ (with the minus sign) admits the following eigen-decomposition

[TABLE]

where $\iota=\sqrt{-1}$. That is, $(\lambda_{n}(\mathbb{S}^{1}),e_{n}):=(n^{2},e^{\iota n\theta})$ is an eigen-pair of $\Delta$, which implies that $\Delta$ is a positive semidefinite and unbounded operator on $L^{2}(\mathbb{S}^{1})$ functions.

Now we can solve the heat equation (18) by expanding it with respect to this orthonormal basis. More precisely, for any $f\in L^{2}(\mathbb{S}^{1})$, the Fourier transform of $f$ is given by

[TABLE]

where $\langle f,g\rangle=(2\pi)^{-1}\int_{\mathbb{S}^{1}}f(\theta)\overline{g(\theta)}\,\mathrm{d}\theta$ is the standard inner product on $\mathbb{S}^{1}$. Then the solution to heat equation (18) with the initial distribution $u(0,\theta)=f(\theta)$ is given by

[TABLE]

So if $\mathbb{S}^{1}$ is insulated, then as $\tau\to\infty$ the heat flow has a constant equilibrium state with the value equal to the average of the initial heat distribution, namely $\lim_{\tau\to\infty}u(\tau,\theta)=(2\pi)^{-1}\int_{\mathbb{S}^{1}}f(\theta)\,\mathrm{d}\theta$. In particular, the second smallest eigenvalue $\lambda_{1}(\mathbb{S}^{1})=1$ characterizes the mixing rate of the heat diffusion process.

Now we consider the “inverse Fourier transform” by expressing the solution $u$ in terms of Green’s function, also called the heat kernel, on $\mathbb{S}^{1}$ as

[TABLE]

where ${\tilde{e}}^{\iota n\theta}=e_{n}/\sqrt{2\pi}$ is the rescaled orthonormal basis of $L^{2}(\mathbb{S}^{1})$. Then we obtain that

[TABLE]

defines an (integral) heat diffusion operator on $\mathbb{S}^{1}$. Then the Laplace operator on $\mathbb{S}^{1}$ can be seen as the generator of the heat diffusion process:

[TABLE]

Similarly, the normalized graph Laplacian $L_{n}=I_{n}-P_{n}=I_{n}-D_{n}^{-1}\mathcal{H}_{n}$ corresponding to the random walk $\mathcal{W}_{n}$ over a random sample of $S_{n}$ on $\mathbb{S}^{1}$ can also be seen as a discrete generator of $\mathcal{W}_{n}$.

Recall that in the heat diffusion process, the second smallest eigenvalue $\lambda_{1}(\mathbb{S}^{1})$ of its generator, i.e., the Laplace-Beltrami operator, characterizes the mixing rate of the heat diffusion process. Similarly, in the random walk $\mathcal{W}_{n}$ (as a discretization of the heat process), the second smallest eigenvalue $\lambda_{j}(L_{n})$ of its discrete generator, i.e., the normalized graph Laplacian operator $L_{n}$, characterizes its mixing rate. From Lemma 6.8, the spectrum of these two operators are related in the sense that for each $j=1,2\dots,$ with probability at least $1-c_{1}n^{-c_{2}}$,

[TABLE]

where $\lambda_{j}(\mathbb{S}^{1})$ is the $j$-th eigenvalue of $\Delta$ and $\lambda_{j}(L_{n})$ is the $j$-th eigenvalue of the normalized graph Laplacian. This means that we must change the time clock unit of the random walk on the graph by multiplying a factor of $h^{2}$ to approximate its underlying heat diffusion process. Thus the term $h^{2}t$ in the second term of (17) is the right time scale $\tau$ for running the heat diffusion process on the manifold, and we need $h^{2}t\to\infty$ for the heat diffusion process converges to an equilibrium distribution. On finite data, this means that the random walk converges to its stationary distribution over the points $S_{n}$ sampled from the submanifold. Using this correspondence, the second term on the left hand side of (17) is a mixing condition on random walks over the $K$ submanifolds.

The first term on the left hand side of (17) can be seen as a separation requirement of the $K$ disjoint submanifolds $\mathcal{M}_{1},\dots,\mathcal{M}_{K}$. In particular, if $t=n^{\epsilon}$ for some $\epsilon>0$ (i.e., we run the random walk in polynomial times/steps), then the minimal separation should obey

[TABLE]

in order to achieve the exact recovery for the manifold clustering problem.

Combining the two terms of (17), we see that steps of the random walk must be properly balanced: we would like the random walk on the similarity graph to sufficiently mix within each cluster (second term of (17)), while it does not overly mix to merge the true clusters (first term of (17)). This reflects the multi-scale property of the diffusion $K$-means. $\blacksquare$

Remark 4.3* (Comparisons with simple thresholding).*

In view of the approximation property (8) of the empirical diffusion affinity matrix $A_{n}=P_{n}^{2t}D_{n}^{-1}$ to a block-diagonal matrix in Remark 3.1, one can show that a simple thresholding of the matrix $A_{n}$ also yields the exact recovery for a properly chosen threshold. Indeed, by the triangle inequality, we have for any $i,j\in G_{k}^{\ast}$ in the same cluster,

[TABLE]

Choosing a threshold value $\gamma$ such that

[TABLE]

one can completely separate the block-diagonal entries from the off-diagonal ones, thus achieving exact recovery. This argument leads to the following lemma.

Lemma 4.4* (Master condition for thresholding to achieve exact recovery).*

If

[TABLE]

then thresholded estimator on the empirical diffusion affinity matrix $A_{n}$ with the threshold value satisfying (21) yields exact recovery.

Note that condition (22) in Lemma 4.4 is slightly weaker than the master condition (34) of the SDP relaxed diffusion $K$-means in Lemma 6.1 (up to a factor of $K^{-1}$ for balanced clusters, say). On the other hand, thresholding has a tuning parameter $\gamma$ and it is unclear how to develop a principled procedure to choose the threshold value satisfying (21). On the contrary, our SDP relaxed diffusion $K$-means only requires the knowledge of the number of cluster $K$ and it is tuning-free in that sense. $\blacksquare$

It is worthy to note that the second smallest eigenvalue of the Laplace-Beltrami operator in condition (17) of Theorem 4.1 can be regarded as characterizations of the connectedness of the submanifolds, where the latter can be formally quantified by the Cheeger isoperimetric constant defined as follows.

Definition 4.5* (Cheeger isoperimetric constant).*

Let $\mathcal{M}$ be a $q$-dimensional compact Riemannian manifold. Let $\mbox{Vol}(\mathcal{A})$ denote the volume of a $q$-dimensional submanifold $\mathcal{A}\subset\mathcal{M}$ and $\mbox{Area}(\mathcal{E})$ denote the $(q-1)$-dimensional area of a submanifold $\mathcal{E}$. The Cheeger isoperimetric constant of $\mathcal{M}$ is defined to be

[TABLE]

where the infimum of the normalized manifold cut (in the curly brackets) is taken over all smooth $(q-1)$-dimensional submanifolds $\mathcal{E}$ of $\mathcal{M}$ that cut $\mathcal{M}$ into two disjoint submanifolds $\mathcal{M}_{1}$ and $\mathcal{M}_{2}$ such that $\mathcal{M}=\mathcal{M}_{1}\bigsqcup\mathcal{M}_{2}$.

In words, $\mathfrak{h}(\mathcal{M})$ quantifies the minimal area of a hypersurface that bisects $\mathcal{M}$ into two disjoint pieces (cf. Figure 2). Smaller values of $\mathfrak{h}(\mathcal{M})$ mean that $\mathcal{M}$ is less connected – in particular, $\mathfrak{h}(\mathcal{M})=0$ implies that there are two disconnected components in $\mathcal{M}$. The Cheeger isoperimetric constant may also be analogously defined for a graph and its value (i.e., the conductance of the graph) is closely related to the normalized graph cut problem. Suppose we have an i.i.d. sample $X_{1},\dots,X_{n}$ drawn from the uniform distribution on $\mathcal{M}$ and $\mathcal{G}_{n}$ is the neighborhood random graph with an edge between $X_{i}$ and $X_{j}$ if $\|X_{i}-X_{j}\|\leq h$. It is shown in [2] that the normalized graph cut (after a suitable normalization) converges to the normalized manifold cut, yielding an asymptotic upper bound on the conductance of $\mathcal{G}_{n}$ based on $\mathfrak{h}(\mathcal{M})$. See also [36] for improved results.

Corollary 4.6*.*

Under the setting of Theorem 4.1, if

[TABLE]

where $\underline{\mathfrak{h}}=\min_{1\leqslant k\leqslant K}\frac{\nu_{q_{k}}\,\mathfrak{h}(\mathcal{M}_{k})}{4\sqrt{q_{k}+2}\,(2\pi)^{q_{k}/4}}$, then $\hat{Z}=Z^{\ast}$ with probability at least $1-c_{3}K\,\underline{n}^{-\xi}$.

Since $\delta$ reflects the separation of the submanifolds of $S=\bigsqcup_{k=1}^{K}\mathcal{M}_{k}$ and $\mathfrak{h}(\mathcal{M}_{k})$ reflects the degree of connectedness of the submanifold $\mathcal{M}_{k}$, the (overall) hardness of the manifold clustering problem is determined by $(\delta,\underline{\mathfrak{h}})$. In particular, if $t=n^{\epsilon}$ for some $\epsilon>0$, then we require that

[TABLE]

in addition to (20). Our results in the rest subsections can also be stated via this geometric quantity of the Cheeger isoperimetric constant.

Remark 4.7* (Comments on the Gaussian kernel and bandwidth parameter).*

In Theorem 4.1 and Corollary 4.6, the kernel $k$ is assumed to be the Gaussian kernel in (1). However, these exact recovery results do not rely on the particular choice of the Gaussian kernel. Specifically, Theorem 4.1 and Corollary 4.6 still hold, as long as the kernel is isotropic and satisfies the exponential decay in the squared Euclidean distance. For the heat kernel on $\mathbb{R}$ (i.e., Green’s function associated with the heat equation on $\mathbb{R}$)

[TABLE]

it can be viewed as an approximation to the short time dynamics of the heat kernel (19) on $\mathbb{S}^{1}$ in Remark 4.2 as $\tau\to 0^{+}$ (cf. Chapter 1 in [28]). Hence, we can approximate the short time behavior of the heat flow on the compact manifold $\mathbb{S}^{1}$ by that of the non-compact manifold $\mathbb{R}$, where the latter is governed by the Gaussian heat kernel on $\mathbb{R}$. Setting $h^{2}=2\tau$ in (24) and noticing that the normalization $(4\pi\tau)^{-1/2}$ does not affect the results in Theorem 4.1 and Corollary 4.6 since the SDP solution in (11) is invariant under scaling. Thus the bandwidth parameter in the Gaussian kernel in (1) has the time scale interpretation in terms of the heat flow dynamics, in addition to capturing the local neighborhood geometry of the submanifolds. $\blacksquare$

4.2. Exact recovery of the regularized and localized diffusion $K$-means

In this subsection, we extend the exact recovery results to the two variants of the diffusion $K$-means. First, we consider the regularized diffusion $K$-means (13) that does not require knowledge of the true number of clusters $K$.

Theorem 4.8* (Exact recovery of SDP for regularized diffusion $K$-means).*

Suppose all conditions in Theorem 4.1 are true. In addition, if the regularization parameter satisfies

[TABLE]

then we can achieve exact recovery for $\tilde{Z}$ from the regularized diffusion $K$-means with probability at least $1-c_{3}K\,\underline{n}^{-\xi}$.

Remark 4.9* (Implication on our $\rho$ selection Algorithm 1).*

Condition (25), as a sufficient condition for the exact recovery, provides some justification of our $\rho$ selection Algorithm 1, in particular, the reason of why using the logarithmic scale. More precisely, we observe from Lemma 3.2 that an upper bound for $\rho$ to produce non-trivial clustering is $n^{-1}$ times the largest eigenvalue of the affinity matrix $A$ in the SDP, which is of order $n^{-1}\max_{k}\{n_{k}h^{q_{k}}\}^{-1}=\mathcal{O}((n\log n)^{-1})$ (from the approximating form (8) and Lemma 6.4 in the proof of Theorem 4.1). As a consequence, in the original scale, the interval length of those $\rho$ that underestimates $K$ is of order $(n\log n)^{-1}$, which is comparable to the range of $\rho$ corresponding to exact recovery (correct $K$) implied by (25) as $(n\log n)^{-1}$. On the other hand, the range of $\rho$ corresponding to exact recovery will dominate if we instead consider the logarithmic scale. Precisely, on the logarithmic scale, the interval length for underestimating $K$ is of order $\log n$, while the range of $\log\rho$ implied by (25) becomes \log n-C^{\prime}\,\big{(}\log n-\min\big{\{}\delta^{2}/(2h^{2}),\min_{1\leqslant k\leqslant K}\{\lambda_{1}(\mathcal{M}_{k})\}\,h^{2}t\big{\}}\big{)}, which is of order $\Omega(n^{\iota})$ for some constant $C^{\prime}>0$ and $\iota>0$ as long as both $\delta^{2}/h^{2}$ and $h^{2}t$ are bounded below by $n^{\iota}$. (The latter requirement can be easily satisfied, for example, in our simulations the choice of $t=n^{1.2}$ is good enough to produce robust results.) In particular, this suggests that the interval length of $\log\rho$ for exact recovery is proportional to $\delta^{2}/h^{2}$, which can be viewed as a signal-to-noise ratio characteristic. $\blacksquare$

Now let us turn to the localized diffusion $K$-means that locally selects the node-wise bandwidth adapting to the local geometric structure.

Theorem 4.10* (Exact recovery of SDP for localized diffusion $K$-means).*

Let $\delta_{kk^{\prime}}=\|\mathcal{M}_{k}-\mathcal{M}_{k^{\prime}}\|$. Fix any $\xi\geqslant 1$. Let $c_{1},c_{2},c_{3},c_{4},M_{0}$ be some positive constants only depending on $\{q_{k},{\rm Sec_{k}},{\rm Inj_{k}},{\rm Rch_{k}},\eta_{k},{\rm Lip_{k}}\}_{k=1}^{K}$, and $C,C^{\prime},C_{1},C_{2}$ only depending on $\{q_{k},{\rm Sec_{k}},{\rm Inj_{k}},{\rm Rch_{k}},\eta_{k}$, ${\rm Lip_{k}}\}_{k=1}^{K}$ and $\xi$. If $\underline{n}\geqslant M_{0}$, the number of neighbor parameter $k_{0}$ satisfies $k_{0}=\lfloor C\log n\rfloor$, and

[TABLE]

for each distinct $k,k^{\prime}\in[K]$, then with probability at least $1-c_{3}K\,\underline{n}^{-\xi}$, the followings are true.

(1) For each $i\in G^{\ast}_{k}$, its local bandwidth parameter $h_{i}$ satisfies

[TABLE]

(2) We can achieve exact recovery for $\tilde{Z}$ from the localized diffusion $K$-means as long as

[TABLE]

Remark 4.11* (Advantages of the localization).*

The first result  Part (1) in Theorem 4.10 shows that our localized selection scheme via nearest neighbors truly leads to bandwidth adaptation to the unknown submanifold dimension $q_{k}$ and the unknown true cluster size $n_{k}$, by only sacrificing a $\log n$ term (from $\log n_{k}$ to $\log n$ as compared with the optimal bandwidth choice $(\log n_{k}/n_{k})^{1/q_{k}}$ from Theorem 4.1). The second result Part (2) in Theorem 4.10 indicates the advantages of using the localized node-wise bandwidth, by comparing the condition  (27) with those in Theorem 4.1. In particular, in order for the lower bound condition on the global bandwidth $h$ in Theorem 4.1 to hold, the smallest $h$ would be $\max_{k}h_{k}$, where $h_{k}=(\log n_{k}/n_{k})^{q_{k}}$ denotes the optimal bandwidth in the $k$-th cluster $\mathcal{M}_{k}$. Note that this lower bound on $h$ is uniformly larger than the magnitudes of localized bandwidth provided in (26). As a consequence, this large $h$ would require the same separation condition as $\delta_{kk^{\prime}}\geqslant\max_{k}h_{k}$ for each pair $(\mathcal{M}_{k},\mathcal{M}_{k^{\prime}})$ of distinct clusters. In comparison, the new sufficient condition (27) for exact recovery only needs a cluster-dependent separation condition as $\delta_{kk^{\prime}}\geqslant\max\{h_{k},h_{k^{\prime}}\}$, which can be substantially weaker than $\delta_{kk^{\prime}}\geqslant\max_{k}h_{k}$ if clusters are highly unbalanced with unequal sizes and mixed dimensions. $\blacksquare$

Finally, we can further combine the regularized diffusion $K$-means with local adaptive bandwidths into the localized and regularized diffusion $K$-means. The following result is an immediate consequence by combining the proofs of Theorem 4.8 and Theorem 4.10, and thus its proof is omitted.

Theorem 4.12* (Exact recovery of SDP for localized and regularized diffusion $K$-means).*

Suppose all conditions in Theorem 4.1 and Theorem 4.10 are true. In addition, if the regularization parameter satisfies

[TABLE]

then we can achieve exact recovery for $\tilde{Z}$ from the localized and regularized diffusion $K$-means with probability at least $1-c_{3}K\,\underline{n}^{-\xi}$.

Remark 4.13* (Comments on bandwidth adaptivity).*

It is interesting to note the signal separation and random walk mixing of the localized diffusion $K$-means and its regularized version are adaptive to local probability density and local geometric structures of the Riemannian submanifolds. In [1], nearly-optimal exact recovery of a collection of clustering methods based on pairwise distances of data is derived under a condition that the minimal signal separation strength over all pairs of submanifolds is larger than a threshold (even for their local scaling version, cf. Proposition 3 therein). Thus results established in [1] are non-adaptive to the local density and (geometric) structures of the submanifolds. In addition, the $\max_{k}\{n_{k}h_{k}^{q_{k}}\}$ on the right hand side of (17) now reduces to $\log(n)$ as in (28). This means that the localized diffusion $K$-means tends to increase the signal-to-noise ratio (cf. Remark 4.11) as well as the upper bound on $\rho$ for exact recovery, thereby widening the interval length of $\log(\rho)$ corresponding to the true clustering structure and improves the performance of the $\rho$ selection Algorithm 1. $\blacksquare$

5. Simulations

In this section, we assess the empirical performance of the diffusion $K$-means on some simulation examples. We generate $n=768$ data points from the following three data generation mechanisms (DGPs).

The clustering structure contains three disjoint submanifolds:

• •

$\mathcal{M}_{1}=\mbox{unit disk}$ and $n/4$ data points are uniformly sampled on $\mathcal{M}_{1}$,

• •

$\mathcal{M}_{2}=\mbox{annulus with radius }2.5$ and $n/4$ data points are uniformly sampled on $\mathcal{M}_{2}$,

• •

$\mathcal{M}_{3}=\mbox{annulus with radius }4$ and $n/2$ data points are uniformly sampled on $\mathcal{M}_{3}$,

where all $\mathcal{M}_{1},\mathcal{M}_{2},\mathcal{M}_{3}$ are centered at the origin $(0,0)$. 2. 2.

The clustering structure contains three disjoint rectangles:

• •

$\mathcal{M}_{1}=\{(-15,-8),(-15,8),(-8,8),(8,8)\}$,

• •

$\mathcal{M}_{2}=\{(10,3),(10,8),(15,3),(15,8)\}$,

• •

$\mathcal{M}_{3}=\{(10,-8),(10,-3),(15,-8),(15,-3)\}$,

where data points are uniformly distributed on $\mathcal{M}_{1}\bigsqcup\mathcal{M}_{2}\bigsqcup\mathcal{M}_{3}$. 3. 3.

The clustering structure is a mixture of three bivariate Gaussians:

[TABLE]

where $(\alpha_{1},\alpha_{2},\alpha_{3})=(1/3,1/3,1/3)$, $\mu_{1}=(-6,0),\mu_{2}=(0,0),\mu_{3}=(2.5,0)$, $\sigma_{1}=2$, and $\sigma_{2}=\sigma_{3}=0.5$.

Our simulation setups are similar to [42, 23]. Note that the sampling density in DGP 1 and 2 is uniform on the disjoint submanifolds, the hardness of the problems is mainly determined by the geometry, and we thus expect the diffusion $K$-means and its localized version can both succeed in these two cases. In addition, since DGP 1 contains two annuli that are less connected than the rectangles and ellipsoids, we expect that, for the localized diffusion $K$-means (with self-tuned bandwidths), more random walk steps are needed for DGP 1 to correctly identify the clusters than those for DGP 2 and DGP 3. In our simulation studies, we use $t=n^{2}$ for DPG 1 and $t=n^{1.2}$ for both DGP 2 and DGP 3 (all with local scaling). Further, DGP 3 has a mixture of Gaussian densities, the local scaling is expected to improve the performance of the diffusion $K$-means. In fact, we have observed in Figure 5 that the diffusion $K$-means without local scaling does not work for DGP 2. It is also known that spectral clustering methods fail on such setup [23]. Thus we do not report results on DGP 2 without local scaling for all competing methods since it does not provide meaningful comparisons with other setups.

For the SDP relaxed diffusion $K$-means clustering methods, we report the $\ell^{1}$ estimation error of for estimating the true clustering membership $Z^{*}$ and the (normalized) Hamming distance error for classifying the clustering labels. For an SDP estimator $\hat{Z}$, the $\ell^{1}$ estimation error is defined as $n^{-1}\|\hat{Z}-Z^{*}\|_{1}$. To get the clustering labels from the SDP estimator $\hat{Z}$, we extract its top $K$ eigenvectors and perform a $K$-means rounding algorithm to get the estimated partition $\hat{G}_{1},\dots,\hat{G}_{K}$. Then the classification error is defined as $n^{-1}\sum_{k=1}^{K}\|\mathbf{1}_{\hat{G}_{k}}-\mathbf{1}_{G^{*}_{k}}\|_{1}$, where $\mathbf{1}_{G_{k}}=(1_{(X_{1}\in G_{k})},\dots,1_{(X_{n}\in G_{k})})^{T}$.

In each setup, our results are reported on 1,000 simulations. For brevity, DKM stands for the diffusion $K$-means, RDKM for the (nuclear norm) regularized diffusion $K$-means, LDKM for the localized diffusion $K$-means, and LRDKM for the localized and regularized diffusion $K$-means. In the cases of no local scaling, the steps of random walks is fixed as $t=n^{1.2}$ in all setups. In the cases of local scaling, the nearest neighborhood size is chosen as $\lfloor\log{n}\rfloor$ for DPG 1 and DGP 3, and as $\lfloor 0.5\log{n}\rfloor$ for DGP 2.

For the comparison purpose, we also include three spectral clustering methods: the unnormalized spectral clustering (SC-UN), the random walk normalized spectral clustering (SC-RWN) [31], a symmetrically normalized spectral clustering (SC-NJW) proposed in [24]. For each spectral clustering method, we also consider their localized versions (LSC-UN, LSC-RWN, LSC-NJW) by replacing the kernel matrix $\mathcal{K}_{n}$ with $\mathcal{K}^{\dagger}_{n}$.

We can draw several observations from the simulation studies. First, the estimation error agrees well with our exact recovery theory for SDP relaxed DKM and LDKM, given the number of clusters (cf. Table 1). Second, all methods works relatively better for DGP 1 and DGP 2 since the separation signal strength is stronger than DGP 3 (cf. Table 2). Third, the RDKM and LRDKM perform well in selecting the true number of clusters (cf. Table 3) using Algorithm 1.

We also modify DGP 3 to make the problem harder. We consider the mixture of three Gaussian with parameters $(\alpha_{1},\alpha_{2},\alpha_{3})=(1/4,1/4,1/2)$, $\mu_{1}=(-6,0),\mu_{2}=(0,0),\mu_{3}=(1.45,0)$, $\sigma_{1}=2$, and $\sigma_{2}=\sigma_{3}=0.5$. This setup is denoted as DGP 3’. For DGP 3’, LDKM has much smaller classification errors than all spectral methods with local scaling (i.e., LSC-UN, LSC-RWN, LSC-NJW); see last column of Table 2.

6. Proofs

Recall that $n_{k}=|G_{k}^{\ast}|$ is the size of $k$-th true cluster index set $G_{k}^{\ast}$, and let $N_{k}=\sum_{i,j\in G_{k}^{\ast}}\kappa(X_{i},\,X_{j})$ denote the total within-weight in $G_{k}^{\ast}$. For any subset $G\subset[n]$, we use $\mathbf{1}_{G}$ to denote the all-one vector whose size equal to the size of $G$.

6.1. Proof of Theorem 4.1

For simplicity of notation, we use $A:=A_{n}$ to denote the empirical diffusion affinity matrix $P_{n}^{2t}D_{n}^{-1}$ in the proof, and recall

[TABLE]

At a high level, our strategy is to show that for suitably large $t\in\mathbb{N}_{+}$, the matrix $A_{n}$ tends to become close to a block-diagonal matrix, where each diagonal block tends to be a constant matrix (cf. equation (8)). Based on this approximation, we expect the global optimum $\hat{Z}$ to share a similar block-diagonal structure, thereby recovers the true membership matrix $Z^{\ast}$ in (12) which takes the form of

[TABLE]

To put this intuition in a technical form, since $Z^{\ast}$ defined in (12) is also a feasible solution belonging to the convex set $\mathscr{C}_{K}$, we have by the optimality of $\hat{Z}$ that

[TABLE]

We analyze the two sums separately as follows.

The first sum: By noticing that $Z^{\ast}_{G_{k}^{\ast}G_{m}^{\ast}}$ is a zero matrix for each pair $k\neq m\in[K]$, we have the following bound

[TABLE]

where the leading factor $\max_{k\neq m}\|[A_{n}]_{G_{k}^{\ast}G_{m}^{\ast}}\|_{\infty}$ is expected to be small due to the approximating structure (8).

The second sum: Since we expect the $k$th block $[A_{n}]_{G_{k}^{\ast}G_{k}^{\ast}}$ of $A_{n}$ in the diagonal to be close to $N_{k}^{-1}\mathbf{1}_{G_{k}^{\ast}}\mathbf{1}^{T}_{G_{k}^{\ast}}$ (cf. (8)), we can subtract and add the same term to decompose it into

[TABLE]

The first term on the right hand side can be bounded by applying Hölder’s inequality,

[TABLE]

where again the leading factor \max_{k}\big{\|}\,[A_{n}]_{G_{k}^{\ast}G_{k}^{\ast}}-N_{k}^{-1}\mathbf{1}_{G_{k}^{\ast}}\mathbf{1}_{G_{k}^{\ast}}^{T}\,\big{\|}_{\infty} is expected to be small due to the approximating structure (8). Now consider the second term. By the definition of $\mathscr{C}_{K}$, the sums of entries in each row of $\hat{Z}$ and $Z^{\ast}$ are equal (to one), we have for fixed $k\in[K]$ and each $i\in G^{\ast}_{k}$,

[TABLE]

Since $Z^{\ast}_{ij}=0$ for each pair $(i,j)$ with $i\in G_{k}^{\ast}$ and $j\notin G_{k}^{\ast}$, and $\hat{Z}$ has nonnegative entries, we can sum up the preceding display over all $i\in G_{k}^{\ast}$ to obtain

[TABLE]

Putting pieces together, we obtain

[TABLE]

Now by combining inequalities (29), (30) and (31), we can reach the following inequality

[TABLE]

Since $\hat{Z}\in\mathscr{C}_{K}\subset\mathscr{C}$, according to inequalities (51)-(53) in Lemma C.1 in Appendix C,

[TABLE]

The last two displays (32) and (33) imply the exact recovery $\hat{Z}=Z^{\ast}$ as long as

[TABLE]

Lemma 6.1* (Master condition for the diffusion $K$-means to achieve exact recovery).*

If (34) holds, then we can achieve exact recovery $\hat{Z}=Z^{*}$.

Proof of Lemma 6.1.

Suppose $\hat{Z}\neq Z^{*}$. Then $\|\hat{Z}-Z^{*}\|_{1}>0$. From (32) and (34), we have

[TABLE]

where the last equality follows from $Z^{*}_{G_{k}^{\ast}G_{m}^{\ast}}=0$ for $m\neq k$. From (33),

[TABLE]

Combining the last two inequalities, we get a contradiction. So we must have $\hat{Z}=Z^{*}$. $\blacksquare$

To further proceed, we will make use of following two lemmas to provide high probability bounds for the empirical diffusion affinity entries deviating from their expectations. Proofs of Lemma 6.2 and 6.3 are deferred to the following subsections.

Lemma 6.2*.*

Let $\kappa=\max_{1\leqslant k\neq k^{\prime}\leqslant K}\sup_{x\in\mathcal{M}_{k},\,x^{\prime}\in\mathcal{M}_{k^{\prime}}}\kappa(x,\,x^{\prime})$ and $\tau=\inf_{x,y\in S:\,\|x-y\|\leqslant h}\kappa(x,y)$. If $c\,(\log n_{k}/n_{k})^{1/q_{k}}\,(\log n_{k})^{\mathbf{1}(q_{k}=2)/4}\leqslant h$, then for any $\xi\geqslant 1$ and any $n_{k}\geqslant M_{0}$,

[TABLE]

holds with probability at least $1-c_{1}\,n_{k}^{-\xi}$, where the spectral gap $\gamma(P_{n,k})$ satisfies

[TABLE]

Here, $\nu_{d}=\frac{\pi^{d/2}}{\Gamma(d/2+1)}$ is the volume of the $d$-dim unit ball, $\lambda_{1}(\mathcal{M}_{k})>0$ is the second smallest eigenvalue of the Laplace-Beltrami operator on $\mathcal{M}_{k}$ (cf. Section 2.3), constants $(c_{1},M_{0})$ depend only on the submanifold $\mathcal{M}_{k}$ through $\{q_{k},{\rm Diam_{k}},{\rm Sec_{k}}$, ${\rm Inj_{k}},{\rm Rch_{k}},{\rm Rch_{k}}\}$ and the density $p_{k}$ through $\{\eta_{k},{\rm Lip_{k}}\}$, and constants $(C,C^{\prime})$ depend only on $\{q_{k},{\rm Sec_{k}},{\rm Inj_{k}},{\rm Rch_{k}}$, $\eta_{k},{\rm Lip_{k}},\xi\}$.

Lemma 6.3*.*

Let $\kappa=\max_{1\leqslant k\neq k^{\prime}\leqslant K}\sup_{x\in\mathcal{M}_{k},\,x^{\prime}\in\mathcal{M}_{k^{\prime}}}\kappa(x,\,x^{\prime})$. Suppose conditions (16) and (17) in Theorem 4.1 are satisfied. Then for any $\xi\geqslant 1$, it holds with probability at least $1-c_{1}K\,\underline{n}^{-\xi}$ that

[TABLE]

where recall $\underline{n}=\min_{k\in[K]}n_{k}$ and $C$ is a constant only depending on $\{q_{k},{\rm Sec_{k}},{\rm Inj_{k}},{\rm Rch_{k}}$, $\eta_{k},{\rm Lip_{k}}\}_{k=1}^{K}$ and $\xi$.

Recall that $\kappa(x,y)=\exp\{-\|x-y\|^{2}/(2h^{2})\}$ is the Gaussian kernel. Consequently, we may choose $\kappa=\exp\{-\delta^{2}/(2h^{2})\}$ and $\tau=e^{-1/2}$ in these two lemmas. By Weyl’s law for the growth of eigenvalues of the Laplace-Beltrami operator (cf. Remark 6 in [34]), we have $\lambda_{j}(\mathcal{M}_{k})\sim j^{2/q_{k}}$. Thus $\lambda_{1}(\mathcal{M}_{k})$ is bounded. Since the bandwidth parameter satisfies $c_{1}(\log n_{k}/n_{k})^{1/q_{k}}\,(\log n_{k})^{\mathbf{1}(q_{k}=2)/4}\leqslant h\leqslant c_{2}$ for some sufficient large constant $c_{1}$ and sufficiently small constant $c_{2}$ for all $k\in[K]$, inequality (35) in Lemma 6.2 becomes

[TABLE]

Therefore, the following two bounds hold with probability at least $1-c_{3}K\,\underline{n}^{-\xi}$,

[TABLE]

for some constants $c_{1},c_{2},c_{3}$ only depending on $\{q_{k},{\rm Sec_{k}},{\rm Inj_{k}},{\rm Rch_{k}},\eta_{k},{\rm Lip_{k}}\}_{k=1}^{K}$ and $C$ only depending on $\{q_{k},{\rm Sec_{k}},{\rm Inj_{k}},{\rm Rch_{k}},\eta_{k},{\rm Lip_{k}}\}_{k=1}^{K}$ and $\xi$. The last two inequalities combining with Lemma 6.1 imply the exact recovery $\hat{Z}=Z^{\ast}$ as long as

[TABLE]

Finally, the claimed result follows from the preceding display and the following lemma of a high probability bound on $N_{k}$, whose proof is postponed to Section 6.9.

Lemma 6.4* (Concentration for total within-weight in $G_{k}^{*}$).*

Suppose the density $p_{k}$ satisfies (15) and the bandwidth $h\leqslant c$ for some constant $c>0$. Then for any $\xi\geqslant 1$, it holds with probability at least $1-c_{2}\,n_{k}^{-\xi}$ that

[TABLE]

where $\beta_{k}=\mathbb{E}_{X\sim p_{k}}[p_{k}(X)]$ and $C$ is a constant only depending on $\{q_{k},{\rm Sec_{k}},{\rm Inj_{k}},{\rm Rch_{k}},\eta_{k},{\rm Lip_{k}}\}_{k=1}^{K}$ and $\xi$.

6.2. Proof of Corollary 4.6

Corollary 4.6 follows from Theorem 4.1 and the Cheeger inequality [6]:

[TABLE]

6.3. Proof of Theorem 4.8

Similar to the proof of Theorem 4.1 in Section 6.1, we can use the optimality of $\tilde{Z}$ and the feasibility of $Z^{\ast}$ for the SDP program (13) to obtain the following basic inequality,

[TABLE]

Since the only place where we used the constraint $\operatorname{tr}(Z)=K$ in Section 6.1 is Lemma C.1 in Appendix C, the analysis of the first two sums in (36) still apply, leading to

[TABLE]

which holds with probability at least $1-c_{3}K\,\underline{n}^{-c_{3}}$.

Now we apply Lemma C.2 in Appendix C to obtain

[TABLE]

Combining this inequality with (37), we obtain

[TABLE]

This implies the exact recovering $\tilde{Z}=Z^{\ast}$ provided that

[TABLE]

Finally, the claimed result follows by combining the above with Lemma 6.4.

6.4. Proof of Lemma 6.2

We consider a fixed $k\in[K]$ throughout this proof. Recall that $P_{n}=D_{n}^{-1}\mathcal{K}_{n}$ defines the random walk $\mathcal{W}_{n}$ over $S_{n}=\{X_{1},X_{2},\ldots,X_{n}\}$, and $A_{n}=P_{n}^{2t}D_{n}^{-1}$. Now consider a new random walk $\mathcal{W}_{n,k}$ over the $k$th cluster $G_{k}^{\ast}$ defined in the following way. For simplicity of notation, we may rearrange the nodes order so that $G_{k}^{\ast}=\{1,2,\ldots,n_{k}\}$. Then the transition probability matrix $P_{n,k}\in\mathbb{R}^{n_{k}\times n_{k}}$ of $\mathcal{W}_{n,k}$ is defined as

[TABLE]

is the induced degree function within cluster $G_{k}^{\ast}$. Similar to the diagonal degree matrix $D_{n}$, we denote by $D_{n,k}\in\mathbb{R}^{n_{k}\times n_{k}}$ the diagonal matrix whose $i$th diagonal entry is $d_{n,k}(X_{i})$ for $i\in[n_{k}]$. Note that $N_{k}=\sum_{i\in G_{k}^{\ast}}d_{n,k}(X_{i})$ the total degrees within $G_{k}^{\ast}$ so that $N_{k}\geqslant n_{k}\min_{i\in G^{\ast}_{k}}d_{n,k}(X_{i})$. It is easy to see that the limiting distribution of $\mathcal{W}_{n,k}$ is $\pi_{n,k}=N_{k}^{-1}\mbox{diag}(D_{n,k})\in\mathbb{R}^{n_{k}}$. Under the separation condition on $\delta$ in the lemma, the probability of moving out from $G_{k}^{\ast}$ is exponentially small, suggesting that we may approximate the sub-matrix $[P_{n}^{2t}]_{G_{k}^{\ast}G_{k}^{\ast}}$ of $P_{n}^{2t}$ with $P_{n,k}^{2t}$. We will formalize this statement in the rest of the proof.

First, we apply the triangle inequality to obtain

[TABLE]

where $T_{1}$ captures the difference between $[P_{n}^{2t}]_{G_{k}^{\ast}G_{k}^{\ast}}$ and $P_{n,k}^{2t}$, and $T_{2}$ characterizes the convergence of the Markov chain $\mathcal{W}_{n,k}$ to its limiting distribution $\pi_{n,k}$ after $2t$ steps.

Recall that $\kappa=\max_{1\leqslant k\neq k^{\prime}\leqslant K}\sup_{x\in\mathcal{M}_{k},\,x^{\prime}\in\mathcal{M}_{k^{\prime}}}\kappa(x,\,x^{\prime})$ is the minimal between-cluster affinity. The first term $T_{1}$ and the second term $T_{2}$ can be bounded via two lemmas below.

Lemma 6.5* (Term $T_{1}$).*

If $n\kappa\leqslant\min_{i\in G_{k}^{*}}d_{n,k}(X_{i})$, then

[TABLE]

Lemma 6.6* (Term $T_{2}$).*

Let $\tau=\inf_{x,y\in S:\,\|x-y\|\leqslant h}\kappa(x,y)$. For each $i\in G_{k}^{\ast}$, let $d_{k}^{\dagger}(X_{i})$ denote total number of points in $\{X_{j}\}_{j\in G_{k}^{\ast}}$ inside the $d$-ball centered at $X_{i}$ with radius $h$. For any $\xi\geqslant 1$, if $c\,(\log n_{k}/n_{k})^{1/q_{k}}\,(\log n_{k})^{\mathbf{1}(q_{k}=2)/4}\leqslant h$, then it holds with probability at least $1-c_{1}\,n_{k}^{-\xi}$ that

[TABLE]

with

[TABLE]

where recall that $\lambda_{1}(\mathcal{M}_{k})>0$ denotes the second smallest eigenvalue of the Laplace-Beltrami operator on $\mathcal{M}_{k}$, and $\nu_{d}=\frac{\pi^{d/2}}{\Gamma(d/2+1)}$ denotes the volume of the $d$-dim unit ball. Here, constant $c_{1}$ depends only on the submanifold $\mathcal{M}_{k}$ through $\{q_{k},{\rm Diam_{k}},{\rm Sec_{k}}$, ${\rm Inj_{k}},{\rm Rch_{k}},{\rm Rch_{k}}\}$ and the density $p_{k}$ through $\{\eta_{k},{\rm Lip_{k}}\}$, and constant $C$ depends on $\{q_{k},{\rm Sec_{k}},{\rm Inj_{k}},{\rm Rch_{k}},\eta_{k},{\rm Lip_{k}},\xi\}$.

Proofs of these two lemmas are provided in Sections 6.6 and 6.7.

Combining upper bounds for the two terms in inequality (38) together, we can reach

[TABLE]

where the spectral gap $\gamma(P_{n,k})$ satisfies  (39). It remains to prove some high probability bounds for $d_{n,k}(X_{i})$ and $d_{k}^{\dagger}(X_{i})$, which are summarized in the following lemma.

Lemma 6.7* (Concentrations for node degrees).*

Suppose the density $p_{k}$ satisfies (15) and the bandwidth $h\leqslant c_{1}$ for some constant $c_{1}>0$. Then for any $\xi\geqslant 1$, it holds with probability at least $1-c_{2}\,n_{k}^{-\xi}$ that

[TABLE]

where $\nu_{d}=\frac{\pi^{d/2}}{\Gamma(d/2+1)}$ denotes the volume of an unit ball in $\mathbb{R}^{d}$. Here constants $c_{1},c_{2}$ only depend on the manifold $\mathcal{M}_{k}$ through $\{q_{k},{\rm Sec_{k}},{\rm Inj_{k}}\}$ and the density function $p_{k}$ through $\{\eta_{k},{\rm Lip_{k}}\}$, and constant $C$ depends on $\{q_{k},{\rm Sec_{k}},{\rm Inj_{k}},\eta_{k},{\rm Lip_{k}},\xi\}$.

A proof of this lemma is deferred to Section 6.8.

Finally, by combining this lemma with the last display, and applying the uniform boundedness condition (15) on $p_{k}$, we obtain that for any $\xi\geqslant 1$, if $c\,(\log n_{k}/n_{k})^{1/q_{k}}\,(\log n_{k})^{\mathbf{1}(q_{k}=2)/4}\leqslant h$, then for any $n_{k}\geqslant M_{0}$,

[TABLE]

holds with probability at least $1-c_{1}\,n_{k}^{-\xi}$, where the spectral gap $\gamma(P_{n,k})$ satisfies

[TABLE]

Here, constants $(c_{1},M_{0})$ depend only on the submanifold $\mathcal{M}_{k}$ through $\{q_{k},{\rm Diam_{k}},{\rm Sec_{k}}$, ${\rm Inj_{k}},{\rm Rch_{k}},{\rm Rch_{k}}\}$ and the density $p_{k}$ through $\{\eta_{k},{\rm Lip_{k}}\}$, and constant $C,C^{\prime}$ depends on $\{q_{k},{\rm Sec_{k}},{\rm Inj_{k}},{\rm Rch_{k}}$, $\eta_{k},{\rm Lip_{k}},\xi\}$.

6.5. Proof of Lemma 6.3

For each indices $i$ and $j$ that belong to two difference clusters $G_{k}^{\ast}$ and $G_{m}^{\ast}$ with $k\neq m$, we have

[TABLE]

Let $\mathcal{W}_{n}=\{Y_{t}:\,t\geqslant 0\}$, with $Y_{t}$ denote the state of the Markov chain $\mathcal{W}_{n}$ at time $t$. Define T_{k}(i)=\min\big{\{}t\in\mathbb{N}_{+}:\,Y_{t}\not\in G_{k}^{\ast},\,Y_{0}=i\big{\}} denote the first exit time from $G_{k}^{\ast}$ of $\mathcal{W}_{n}$ starting from $Y_{0}=i$. Then it is easy to see that

[TABLE]

Since for each $i\in G_{k}^{\ast}$, the one-step transition probability of moving out from $G_{k}^{\ast}$ is bounded by $n\,\kappa\,\max_{i\in G_{k}^{\ast}}d_{n}^{-1}(X_{i})$, we have

[TABLE]

provided that $n\,\kappa\leqslant\min_{i\in G_{k}^{\ast}}d_{n}(X_{i})$. Therefore, for each $i\in G_{k}^{\ast}$ and $j\in G_{m}^{\ast}$ with $k\neq m$, we have

[TABLE]

By Lemma 6.7 and (15), we have that with probability at least $1-c_{1}n^{-\xi}$, $C_{1}n_{k}h^{q_{k}}\leqslant d_{n,k}(X_{i})\leqslant C_{2}n_{k}h^{q_{k}}$ for some constants $C_{1}$ and $C_{2}$ only depending on $\mathcal{M}_{k}$ and $p_{k}$. Note that condition (17) in Theorem 4.1 yields that $n\kappa=ne^{-\delta^{2}/(2h^{2})}\leqslant Cn_{k}h^{q_{k}}$. Then the claimed bound is implied by the above combined with the union bound.

6.6. Proof of Lemma 6.5

By adding and subtracting the same term we obtain

[TABLE]

where the second inequality is due to the fact that each row sum of $[P_{n}^{2t}]_{G_{k}^{\ast}G_{k}^{\ast}}$ is at most one.

Now we consider the first term in (40). Recall that $\mathcal{W}_{n}=\{Y_{t}:\,t\geqslant 0\}$, where $Y_{t}$ is the state of the Markov chain $\mathcal{W}_{n}$ at time $t$. Note that for each $i,j\in G_{k}^{\ast}$, we have

[TABLE]

As a consequence, we have by the law of total probability and the Markov property of $\mathcal{W}_{n}$ that for each $i,j\in G_{k}^{\ast}$ and $s\in\mathbb{N}_{+}$,

[TABLE]

For each pair $(j,\ell)$ belonging to different clusters, noting that $d_{n,k}(X_{i})\leqslant d_{n}(X_{i})$, we have

[TABLE]

which implies that for each $\ell\in G_{k}^{\ast}$,

[TABLE]

Combining the last three displays, we obtain for each $i,j\in G_{k}^{\ast}$ and $s\in\mathbb{N}^{+}$,

[TABLE]

which can be further simplified into

[TABLE]

Now we can recursively apply this two-sided inequality to get

[TABLE]

By taking $s=2t$ and applying the inequality $(1-x)^{s}\geqslant 1-xs$ for $s\in\mathbb{N}_{+}$ and $x\in[0,1]$, the above can be further reduced into

[TABLE]

Then we get

[TABLE]

which is an upper bound to the first term in inequality (40) for $T_{1}$.

The second term in inequality (40) can be bounded as

[TABLE]

Finally, by combining (40), (41) and (42), we obtain

[TABLE]

6.7. Proof of Lemma 6.6

Recall that $\gamma(P_{n,k})=1-\lambda_{1}(P_{n,k})$ denote spectral gap of the transition matrix $P_{n,k}$, where $\lambda_{1}(P_{n,k})$ denotes the second largest eigenvalue of $P_{n,k}\in\mathbb{R}^{n_{k}\times n_{k}}$ (due to similar arguments as in Appendix B, $P_{n,k}$ has $n_{k}$ real eigenvalues with $1$ as the largest one). In addition, since the kernel function $k$ is positive semidefinite, all eigenvalues of $P_{n,k}$ are nonnegative, meaning that $\gamma(P_{n,k})$ is equal to the absolute spectral gap $1-\max\{\lambda_{1}(P_{n,k}),\,\lambda_{n_{k}-1}(P_{n,k})\}$ where $\lambda_{n_{k}-1}(P_{n,k})$ is the $n_{k}$th (smallest) eigenvalue of $P_{n,k}$.

Therefore, according to the relationship between the mixing time of a Markov chain and its absolute spectral gap (see, for example, equation (12.11) in [15]), we have for each $i,j\in G_{k}^{\ast}$,

[TABLE]

where \pi_{n,k}=\big{(}d_{n,k}(X_{1})/N_{k},\ldots,d_{n,k}(X_{n_{k}})/N_{k}\big{)}^{T}\in\mathbb{R}^{n_{k}} is the limiting distribution of induced Markov chain $\mathcal{W}_{n,k}$ over $G_{k}^{\ast}$ with transition probability matrix $P_{n,k}$. This leads to a bound on $T_{2}$ as

[TABLE]

Therefore, it remains to provide a lower bound on the spectral gap $\gamma(P_{n,k})=1-\lambda_{1}(P_{n,k})$. We do so by applying a comparison theorem of Markov chains (Lemma 13.22 in [15]), where we compare the spectral gap of $P_{n,k}$ with that of a standard random walk on a random geometric graph over $\{X_{i}\}_{i\in G_{k}^{\ast}}$, where each pair of nodes are connected if and only if they are at most $h$ far away from each other. The spectrum of the normalized graph Laplacian of the latter is known to behave like the eigensystem of the Laplace-Beltrami operator over the submanifold corresponding to the $k$-th connected subset $\mathcal{M}_{k}$. In particular, we will use existing results [4, 34] on error estimates by using the spectrum of a random geometric graph to approximate the eigensystem of the Laplace-Beltrami operator in the numerical analysis literature.

Let us first formally define a random geometric graph over $\{X_{i}\}_{i\in G_{k}^{\ast}}$ as i.i.d. samples from the compact connected $q_{k}$-dimensional Riemannian submanifold $\mathcal{M}_{k}$ in $\mathbb{R}^{p}$ with bounded diameter, absolute sectional curvature value, and injectivity radius. Recall that $\mu_{k}$ is a probability measure on $\mathcal{M}_{k}$ that has a Lipschitz density $p_{k}$ with respect to the Riemannian volume measure on $\mathcal{M}_{k}$ satisfying (15). $\{X_{i}\}_{i\in G_{k}^{\ast}}$ can be viewed as a sequence of i.i.d. samples from $\mu_{k}$, and without loss of generality, we may assume $G_{k}^{\ast}=\{1,2,\ldots,n_{k}\}$. Consider the random geometric graph $\mathcal{G}_{k}^{\dagger}=(V_{k},E_{k})$, with $V_{k}=\{X_{i}\}_{i\in G_{k}^{\ast}}$ being its set of vertices and $E_{k}$ set of edges, constructed by putting an edge between $X_{i}$ and $X_{j}$ (write $i\sim j$ and call $X_{i}$ to be a neighbor of $X_{j}$) if and only if $\|X_{i}-X_{j}\|\leqslant h$. We define the natural random walk $\mathcal{W}_{k}^{\dagger}$ as a reversible Markov chain on $V_{k}$ that moves to neighbors of the current state with equal probabilities. In other words, the transition probability matrix $\mathcal{P}_{k}^{\dagger}\in\mathbb{R}^{n_{k}\times n_{k}}$ satisfies

[TABLE]

where $d^{\dagger}_{k,i}:=d^{\dagger}_{k}(X_{i})=\sum_{j=1}^{n_{k}}1(j\sim i)$ denotes the degree of vertex $X_{i}$. It is easy to see that $\pi^{\dagger}_{k}=(d^{\dagger}_{k,1}/d^{\dagger}_{k},d^{\dagger}_{k,2}/d^{\dagger}_{k},\ldots,d^{\dagger}_{k,n_{k}}/d^{\dagger}_{k})^{T}$, where $d^{\dagger}_{k}=\sum_{i=1}^{n_{k}}d^{\dagger}_{k,i}$ denotes the total degree, is the stationary distribution of this random walk. Let $1=\lambda_{0}(\mathcal{P}_{k}^{\dagger})\geqslant\lambda_{1}(\mathcal{P}_{k}^{\dagger})\geqslant\ldots\geqslant\lambda_{n-1}(\mathcal{P}_{k}^{\dagger})\geqslant-1$ denote the eigenvalues of matrix $\mathcal{P}_{k}^{\dagger}$, and $\gamma(\mathcal{P}_{k}^{\dagger})=1-\lambda_{1}(\mathcal{P}_{k}^{\dagger})$ denote its spectral gap.

Let $L_{\mathcal{G}_{k}^{\dagger}}=D_{k}^{\dagger}-A_{k}^{\dagger}\in\mathbb{R}^{n_{k}\times n_{k}}$ denote the graph Laplacian matrix associated with graph $\mathcal{G}_{k}^{\dagger}=(V_{k},E_{k})$, where $D_{k}^{\dagger}\in\mathbb{R}^{n_{k}\times n_{k}}$ is a diagonal matrix with $[D^{\dagger}_{k}]_{ii}=d_{k,i}^{\dagger}$, and $A_{k}^{\dagger}\in\mathbb{R}^{n_{k}\times n_{k}}$ is the adjacency matrix with $[A_{k}^{\dagger}]_{ij}=1(i\sim j)$ for all distinct pair $(i,j)\in[n_{k}]^{2}$. Define the normalized Laplacian of $\mathcal{G}_{k}^{\dagger}$ as $L^{N}_{\mathcal{G}_{k}^{\dagger}}=(D_{k}^{\dagger})^{-1/2}L_{\mathcal{G}_{k}^{\dagger}}(D_{k}^{\dagger})^{-1/2}=I-(D_{k}^{\dagger})^{-1/2}A_{k}^{\dagger}(D_{k}^{\dagger})^{-1/2}$, and denote its ordered eigenvalues by $0\leqslant\lambda_{0}(L^{N}_{\mathcal{G}_{k}^{\dagger}})\leqslant\lambda_{1}(L^{N}_{\mathcal{G}_{k}^{\dagger}})\leqslant\cdots\leqslant\lambda_{n_{k}-1}(L^{N}_{\mathcal{G}_{k}^{\dagger}})$. Since $(D_{k}^{\dagger})^{-1/2}A_{k}^{\dagger}(D_{k}^{\dagger})^{-1/2}=(D_{k}^{\dagger})^{1/2}\mathcal{P}_{k}^{\dagger}(D_{k}^{\dagger})^{-1/2}$ is a similarity transformation of $\mathcal{P}_{k}^{\dagger}$, they share the same eigenvalues. Therefore, we have the relation $\lambda_{j}(L^{N}_{\mathcal{G}_{k}^{\dagger}})=1-\lambda_{j}(\mathcal{P}_{k}^{\dagger})$ for all $j=0,1,\ldots,n_{k}-1$. In particular, by taking $j=1$, we can relate the spectral gap of $\mathcal{P}_{k}^{\dagger}$ with the second smallest eigenvalue of the normalized Laplacian matrix $L^{N}_{\mathcal{G}_{k}^{\dagger}}$ as $\gamma(\mathcal{P}_{k}^{\dagger})=\lambda_{1}(L^{N}_{\mathcal{G}_{k}^{\dagger}})$.

It is known that the eigenvalues of the normalized Laplacian $L^{N}_{\mathcal{G}_{k}^{\dagger}}$ of the geometric random graph $\mathcal{G}_{k}^{\dagger}$ approaches (up to a scaling factor) the eigenvalues of the drift Laplace-Beltrami operator on $\mathcal{M}_{k}$. More concretely, let $L^{2}(\mathcal{M}_{k},\mathrm{d}\mu_{k})$ be the space of all square integrable functions on $\mathcal{M}_{k}$, and $\Delta_{\mathcal{M}_{k}}$ denote the drift Laplace-Beltrami operator on $\mathcal{M}_{k}$ defined in (6) in Section 2.3. Let $0=\lambda_{0}(\mathcal{M}_{k})\leqslant\lambda_{1}(\mathcal{M}_{k})\leqslant\cdots$ denote the sequence of nonnegative eigenvalues of $\Delta_{\mathcal{M}_{k}}$. The connectedness of $\mathcal{M}_{k}$ implies that its second smallest eigenvalue $\lambda_{1}(\mathcal{M}_{k})$ is strictly positive. We will invoke Corollary 2 of [34], which generalizes Theorem 1 of [4] from the uniform density to any Lipschitz continuous density satisfying (15), for relating the spectrum of the drift Laplace-Beltrami operator $\Delta_{\mathcal{M}_{k}}$ on $\mathcal{M}_{k}$ with the spectrum of the discrete normalized graph Laplacian $L^{N}_{\mathcal{G}_{k}^{\dagger}}$, as summarized in the following.

Lemma 6.8* (Convergence of eigenvalues of normalized graph Laplacian).*

Let $\nu_{q_{k}}$ to denote the volume of an unit ball in $\mathbb{R}^{q_{k}}$. For each $j=0,1,\ldots$, suppose the radius $h$ and the $\varepsilon_{n,k}$ to be defined below satisfy $(\sqrt{\lambda_{j}(\mathcal{M}_{k})}+1)\,h+\varepsilon_{n,k}/h\leqslant c_{1}$. Then for any $\xi\geqslant 1$, it holds with probability at least $1-c_{2}\,n_{k}^{-\xi}$

[TABLE]

where constants $c_{1},c_{2},c_{3},C>0$ depend only on the submanifold $\mathcal{M}_{k}$ through $\{q_{k},{\rm Diam_{k}},{\rm Sec_{k}}$, ${\rm Inj_{k}},{\rm Rch_{k}}\}$ and the density $p_{k}$ through $\{\eta_{k},{\rm Lip_{k}}\}$, and

[TABLE]

In particular, this lemma (taking $j=1$) implies a lower bound on the spectral gap of $\mathcal{P}_{k}^{\dagger}$ as

[TABLE]

as long as $h\geqslant c_{1}\varepsilon_{n,k}$.

Next, we will apply the following comparison theorem to relate the spectral gaps of Markov chains $\mathcal{W}_{n,k}$ and $\mathcal{W}_{k}^{\dagger}$.

Lemma 6.9* (Markov chains comparison theorem (Lemma 13.22 in [15])).*

Let $P$ and $P^{\prime}$ be transition matrices of two reversible Markov chains on the same state space $\Omega$, whose stationary distributions are denoted by $\pi$ and $\pi^{\prime}$, respectively. Let $\mathcal{E}(f)$ and $\mathcal{E}^{\prime}(f)$ denote the Dirichlet forms associated to the pairs $(P,\pi)$ and $(P^{\prime},\pi^{\prime})$, where

[TABLE]

and $\mathcal{E}^{\prime}(f)$ can be similarly defined. If there exists some constant $B>0$ such that $\mathcal{E}^{\prime}(f)\leqslant B\,\mathcal{E}(f)$ for all $f$, then

[TABLE]

where $\gamma(P)$ and $\gamma(P^{\prime})$ denote the spectral gaps associated with $P$ and $P^{\prime}$, respectively.

We will apply this comparison theorem with $P_{n,k}\to P,\mathcal{P}_{k}^{\dagger}\to P^{\prime}$, and $\Omega=G_{k}^{*}$. Let us find the constant $B$ such that $\mathcal{E}^{\prime}(f)\leqslant B\,\mathcal{E}(f)$ for all $f$, where in our setting,

[TABLE]

According to the definition of $\tau$ as $\inf_{x,y\in S:\,\|x-y\|\leqslant h}\kappa(x,y)$, we can simply choose $B=(N_{k}/d_{k}^{\dagger})\,\tau^{-1}$. In addition, we have the bound

[TABLE]

Therefore, we can apply Lemma 6.9 to get

[TABLE]

which combined with inequality (44) leads to

[TABLE]

6.8. Proof of Lemma 6.7

Recall that $\kappa(x,y)=\exp\{-\|x-y\|^{2}/(2h^{2})\}$ is the Gaussian kernel with bandwidth parameter $h>0$. For $x\in S$, define $d_{n,k}(x)=\sum_{j\in G_{k}^{*}}\kappa(x,X_{j})$ as the induced degree function of $x$ within cluster $G_{k}^{*}$. Then for each $i\in G_{k}^{\ast}$ we have $d_{n,k}(X_{i})=1+\sum_{j\in G_{k}^{\ast},\,j\neq i}\kappa(X_{i},X_{j})=:1+\tilde{d}_{n,k}(X_{i})$. Denote $\alpha_{k}(x)=\operatorname{\mathds{E}}_{X\sim p_{k}}\kappa(x,\,X)$ and $v_{k}(x)=\text{Var}_{X\sim p_{k}}[\kappa(x,\,X)]$. Applying Lemma 6.10 with $\mathcal{M}=\mathcal{M}_{k}$ and $f(x)=p_{k}(x)$, we have $v_{k}(x)\leqslant C\alpha_{k}(x)\leqslant Ch^{q_{k}}$ for all $x\in\mathcal{M}_{k}$, where constant $C$ only depends on $\{q_{k},{\rm Sec_{k}},{\rm Inj_{k}},\eta_{k},{\rm Lip_{k}}\}$. Then for any fixed $x\in\mathcal{M}_{k}$, using the bound in (49) and the boundedness of $\kappa$, we may apply Bernstein inequality (cf. Lemma 2.2.9 in [37]) to obtain that for all $u>0$,

[TABLE]

Choosing $u=tn_{k}\alpha_{k}(x)$ for $t\in(0,C]$, we have

[TABLE]

Now choosing $t=c_{1}\sqrt{\log(n_{k})/(\alpha_{k}(x)n_{k})}$ for some large enough constant $c_{1}$ so that $Cc_{1}^{2}\geqslant\xi$, we get

[TABLE]

provided that $\log(n_{k})\leqslant C\alpha_{k}(x)n_{k}\leqslant Cn_{k}h^{q_{k}}$ in view of the uniform bounds (48) and (15). But this is ensured by our bandwidth assumption (16). Thus for any fixed $x\in\mathcal{M}_{k}$, we have with probability at least $1-c_{2}n_{k}^{-c_{3}}$,

[TABLE]

Then it follows that with probability at least $1-c_{2}n_{k}^{-\xi}$,

[TABLE]

This implies that the rescaled degree function n_{k}^{-1}\,\big{(}\sqrt{2\pi}h\big{)}^{-q_{k}}\,d_{n,k}(x) provides a good estimate of the density $p_{k}(x)$ at $x$. Since $X_{i}\in\mathcal{M}_{k}$ are i.i.d. for $i\in G_{k}^{*}$, we have with probability at least $1-c_{2}n_{k}^{-\xi}$,

[TABLE]

Then union bound implies that

[TABLE]

with probability at least $1-c_{2}n_{k}^{-\xi}$.

The second part about the concentration of $d_{k}^{\dagger}(X_{i})$ can be analogously proved by applying the Chernoff bound for sum of i.i.d. Bernoulli random variables. Let

[TABLE]

so that $d_{k}^{\dagger}(X_{i})=1+\tilde{d}_{k}^{\dagger}(X_{i})$ where $\tilde{d}_{k}^{\dagger}(X_{i})=\sum_{j\in G_{k}^{*},j\neq i}1(\|X_{i}-X_{j}\|\leqslant h)$. Note that Section 2.2 of [4] provides a uniform estimate of the expectation $\mathbb{E}_{X\sim p_{k}}[1(\|x-X\|\leqslant h)]$ in terms of the density $p(x)$ as

[TABLE]

where recall that $\nu_{q_{k}}$ denotes the volume of unit ball in $\mathbb{R}^{q_{k}}$. The rest of the proof follows a similar line as the proof of the first part, and we omit the details.

Lemma 6.10*.*

Let $\mathcal{M}$ be a $q$-dimensional compact submanifold in $\mathbb{R}^{p}$ with bounded absolute sectional curvature ${\rm Sec}$ and injective radius ${\rm Inj}$, and $\mbox{Vol}_{\mathcal{M}}$ denote its volume form. Let $f$ be a ${\rm Lip}$-Lipschitz continuous probability density function on $\mathcal{M}$ such that $\eta\leqslant f(x)\leqslant\eta^{-1}$ for some constant $\eta>0$. Let $\alpha(x)=\operatorname{\mathds{E}}_{X\sim f}[\kappa(x,\,X)]$ and $v(x)=\text{Var}_{X\sim f}[\kappa(x,\,X)]$. Then we have

[TABLE]

and $v(x)\leqslant C\alpha(x)$ for all $x\in\mathcal{M}$, where the constant $C$ only depends $\mathcal{M}$ through $\{q,{\rm Sec},{\rm Inj}\}$ and $f$ through $\{\eta,{\rm Lip}\}$. Consequently, we have $\sup_{x\in\mathcal{M}}\alpha(x)\leqslant Ch^{q}$ and $\sup_{x\in\mathcal{M}}v(x)\leqslant Ch^{q}$.

Proof of Lemma 6.10.

Note that for each $x\in\mathcal{M}$ and any Lipschitz probability density function $f$ on $\mathcal{M}$, the expectation $\mathbb{E}[\kappa(x,\,X)^{2}]$ takes the form

[TABLE]

Then (48) follows from Lemma 6.11. Similarly, we can bound the variance of $\kappa(x,\,X)$ via $\mbox{Var}_{X\sim f}[\kappa(x,\,X)]\leqslant\mathbb{E}_{X\sim f}[\kappa(x,\,X)^{2}]$, where

[TABLE]

By a second application of Lemma 6.11 with $h/\sqrt{2}\to h$ to obtain

[TABLE]

This together with the uniform boundedness condition on $f$ and inequality (48) imply an upper bound to the variance by the expectation,

[TABLE]

The bounds $\sup_{x\in\mathcal{M}}\alpha(x)\leqslant Ch^{q}$ and $\sup_{x\in\mathcal{M}}v(x)\leqslant Ch^{q}$ follow from the fact that $h\leqslant c$. $\blacksquare$

Lemma 6.11* (Convolution over manifold).*

Let $\mathcal{M}$ be a $q$-dimensional compact submanifold in $\mathbb{R}^{p}$ with bounded absolute sectional curvature ${\rm Sec}$ and injective radius ${\rm Inj}$, and $\mbox{Vol}_{\mathcal{M}}$ denote its volume form., and $\mbox{Vol}_{\mathcal{M}}$ denote its volume form. Then there exists some constant $h_{0}>0$ only depending on $\{{\rm Sec},{\rm Inj}\}$, such that for all $h\in(0,h_{0}]$ and any ${\rm Lip}$-Lipschitz continuous function $f$ on $\mathcal{M}$, we have

[TABLE]

where constant $C$ only depends on $\{{\rm Sec},{\rm Inj},{\rm Lip}\}$.

Proof of Lemma 6.11.

This lemma follows from Lemma 5.2 on pages 895–898 in [41]. $\blacksquare$

6.9. Proof of Lemma 6.4

Recall that $N_{k}=\sum_{i,j\in G^{\ast}_{k}}\kappa(X_{i},X_{j})$ is the total within-weight in $G^{\ast}_{k}$. According to Lemma 6.7 (note that $N_{k}=\sum_{i\in G_{k}^{\ast}}d_{n,k}(X_{i})$ using the notation therein), it holds with probability at least $1-c_{2}n_{k}^{-\xi}$ that

[TABLE]

According to the sandwiched bound on the density function $p_{k}$, $\{p_{k}(X_{i}):\,i\in G_{k}^{\ast}\}$ are independent and bounded random variables. Therefore, we may apply Hoeffding’s inequality to obtain that

[TABLE]

holds with probability at least $1-c_{2}n_{k}^{-\xi}$, where $\beta_{k}=\mathbb{E}[p_{k}(X_{i})]$. Combining the two preceding inequalities, we obtain that

[TABLE]

which completes the proof.

6.10. Proof of Theorem 4.10

Proof of Part (1):

Consider $X_{i}$, where $i\in G_{k}^{\ast}$ for some $k\in[K]$. Fix $h_{U}=c_{U}(\log n/n_{k})^{1/q_{k}}$ and $h_{L}=c_{L}(\log n/n_{k})^{1/q_{k}}$ for two sufficiently large constant $c_{U}$ and $c_{L}$ with $c_{U}=2c_{L}$. We use notation $N(X_{i},h)$ to denote the number of points from $\{X_{i}\}_{i=1}^{n}$ that is within $h$ distance from $X_{i}$.

Recall that $d_{k}^{\dagger}(X_{i}):=d_{k,h}^{\dagger}(X_{i})$ in Lemma 6.7 is the number of points from $\{X_{i}\}_{i\in G_{k}^{\ast}}$ that is within $h$ distance to $X_{i}$ (we will choose $h=h_{L}$ and $h=h_{U}$ later). Here we put an subscript $h$ in $d_{k,h}^{\dagger}(X_{i})$ to indicate the dependence of $d_{k}^{\dagger}$ on $h$. The condition on $\delta_{kk^{\prime}}$ implies that any point outside the $k$th cluster $D_{k}$ has distance at least $C^{\prime}(\log n/n_{k})^{1/q_{k}}\geqslant\max\{h_{L},h_{U}\}$. Therefore all points that are within $h_{L}(h_{U})$ distance to $X_{i}$ must belong to $\mathcal{M}_{k}$, implying $N(X_{i},h)=d_{k,h}^{\dagger}(X_{i})$. Consequently, from the proof of Lemma 6.7 (take $h=h_{U}$ and $h=h_{L}$ respectively), we have (a concentration inequality as (46) plus the expectation bound (47)),

[TABLE]

where constants $(c_{1},c_{1}^{\prime})$ and $(c_{2},c_{2}^{\prime})$ only depend on the constant $\eta_{k}$ in the two-sided bound on density $p_{k}(\cdot)$ on the $k$th region $\mathcal{M}_{k}$. Let $c_{U}$ be sufficiently large such that $c_{2}c_{U}^{q_{k}}=C$, where $C$ is the constant appearing in the neighborhood parameter $k_{0}=\lfloor C\log n\rfloor$. By taking $t$ such that $(1-t)c_{1}c_{L}^{q_{k}}=C/2$ in the first inequality of the preceding display, and $t=1$ in the second inequality, we obtain that

[TABLE]

where $c_{1}^{\prime\prime}$ and $c_{2}^{\prime\prime}$ are two constants independent of $C$. For large enough $C$, we can make these two probabilities smaller than $n^{-2-\xi}$. Since the event $\{N(X_{i},h_{L})\leqslant k_{0}/2\}\cap\{N(X_{i},h_{U})\geqslant 2k_{0}\}$ implies $h_{i}\in[h_{L},h_{U}]$ for each $i\in[n]$, a union bound argument over $i=1,2,\ldots,n$ leads to the claimed two sided bound for $h_{i}$ (with probability at least $1-n^{-\xi}$).

Proof of Part (2):

Using the two sided bound in the proof of Part (1), the proof follows same steps in the proof of Theorem 4.1, with the only exception in proving a counterpart of Lemma 6.6 for bounding the spectral gap $\gamma(P_{n,k})$ of chain $P_{n,k}$ on $\{X_{i}\}_{i\in G_{k}^{\ast}}$ in cluster $G_{k}^{\ast}$ for $k\in[K]$, since $P_{n,k}$ now has different bandwidth parameter $h_{i}$ at each observed point $X_{i}$ (that is, a bandwidth parameter inhomogeneous chain). More precisely, it remains to show that with high probability, it holds that

[TABLE]

Here, we assume without loss of generality that the absolute spectral gap of $P_{n,k}$ is dominated by one minus its second largest eigenvalue. Otherwise, we can always consider the lazy random walk by replacing $P_{n}$ with $P_{n}/2+I_{n}/2$ in the diffusion $K$-mean SDP, whose absolute spectral gap is $\gamma(P_{n,k})/2$.

Our proof strategy is again based on the Markov chains comparison theorem (Lemma 6.9) by comparing this bandwidth parameter inhomogeneous chain with a bandwidth parameter homogeneous chain with $h_{i}\equiv h_{L}$, for each $i\in G_{k}^{\ast}$, where $h_{L}$ is the lower bound of $h_{i}$ in the proof of Part (1). In particular, a lower bound on the spectral gap of the latter is already derived in Lemma 6.6 as of order $\lambda_{1}(\mathcal{M}_{k})h_{L}^{2}$.

Fix the cluster index $k\in[K]$, and without loss of generality assume $G_{k}^{\ast}=\{1,2,\ldots,n_{k}\}$. To avoid confusion of notation, we put an superscript “IH” indicate the bandwidth parameter inhomogeneous chain, and “H” to denote the bandwidth parameter homogeneous chain with $h_{i}\equiv h_{L}$. For example, $P_{n,k}^{IH}$ and $\mathcal{E}^{IH}$ denote the transition probability matrix and the Dirichlet form (45) (defined in Lemma 6.9), respectively, associated with the bandwidth parameter inhomogeneous chain.

Due to the fact that $h_{i}\geqslant h_{L}$ for all $i\in G^{\ast}_{k}$, we immediately have

[TABLE]

where $k^{H}(X_{i},X_{j})=\exp\{-\|X_{i}-X_{j}\|^{2}/(2h_{L}^{2})\}$ and $\kappa^{IH}(X_{i},X_{j})=\exp\{-\|X_{i}-X_{j}\|^{2}/(2h_{i}h_{j})\}$, and $(N_{k}^{H},N_{k}^{IH})$ are the respective total degrees within $G_{k}^{\ast}$. Moreover, recall that the stationary distributions of $P_{n,k}^{IH}$ and $P_{n,k}^{H}$ are $\pi_{n,k}^{IH}=d_{n,k}^{IH}(X_{i})/N_{k}^{IH}$ and $\pi_{n,k}^{H}=d_{n,k}^{H}(X_{i})/N_{k}^{H}$, for $i\in G_{k}^{\ast}$, respectively, where $d_{n,k}^{IH}(X_{i})=\sum_{j\in G_{k}^{\ast}}\kappa^{IH}(X_{i},X_{j})$, $d_{n,k}^{H}(X_{i})=\sum_{j\in G_{k}^{\ast}}\kappa^{H}(X_{i},X_{j})$ are the node degrees, and $N_{k}^{IH}=\sum_{i\in G_{k}^{\ast}}d_{n,k}^{IH}(X_{i})$, $N_{k}^{H}=\sum_{i\in G_{k}^{\ast}}d_{n,k}^{H}(X_{i})$ are the total degrees.

Now we can apply Lemma 6.9 with $B=N_{k}^{IH}/N_{k}^{H}$ and Lemma 6.6 (to the homogeneous chain $P_{n,k}^{H}$) to obtain

[TABLE]

for some constant $C>0$. Similar to the proof of Lemma 6.7, we can apply the concentration inequality to the nodes degree with bandwidth $h=h_{U}$ and $h=h_{L}$ to obtain that with probability at least $1-c_{2}\,n^{-\xi}$,

[TABLE]

for all $i\in G_{k}^{\ast}$. Combining this with the sandwiched bound for $h_{i}$ in Part (1), we obtain

[TABLE]

Putting all pieces together, we obtain that it holds with probability at least $1-c_{2}\,n^{-\xi}$ that

[TABLE]

6.11. Proof of Lemma 3.2

Proof of Part (1): Since $n\rho>\lambda_{\max}(A)$, the matrix $n\rho I_{n}-A$ is positive definite. For any $Z\in\mathscr{C}$, from the constraint $Z\mathbf{1}_{n}=\mathbf{1}_{n}$, we see that $(1,n^{-1/2}\mathbf{1}_{n})$ is an eigen-pair of $Z$. In addition, since $Z\succeq 0$, all eigenvalues $\lambda_{1},\dots,\lambda_{n}$ of $Z$ are non-negative. Let $U_{1}=n^{-1/2}\mathbf{1}_{n},U_{2},\ldots,U_{n}$ denote the corresponding eigenvectors of $Z$. Thus the objective function

[TABLE]

where the equality holds if and only if $\lambda_{2}=\cdots=\lambda_{n}=0$. Note that $Z^{\diamond}\in\mathscr{C}$ is a feasible solution for (13) and $Z^{\diamond}$ has a non-zero eigenvalue equal to $1$ and $(n-1)$ zero eigenvalues. Therefore, $Z^{\diamond}$ is the unique solution of the SDP (13).

Proof of Part (2): For any $Z\in\mathscr{C}$, since $Z$ is a symmetric matrix satisfying $Z\geqslant 0$ and $Z\mathbf{1}_{n}=\mathbf{1}_{n}$, $Z$ is a stochastic matrix and all its eigenvalues have absolute values less than or equal to one. Moreover, from the positive semi-definiteness of $Z$, all eigenvalues of $Z$ lie in the $[0,1]$ interval. Now since $n\rho<\lambda_{\min}(A)$, the matrix $A-n\rho I_{n}$ is positive definite. From matrix Hölder’s inequality, the objective function satisfies

[TABLE]

where the equality holds if and only if all eigenvalues of $Z$ equal to one (since $A-n\rho I_{n}$ is strictly positive definite). Note that $Z^{\dagger}=I_{n}\in\mathscr{C}$ is a feasible solution for (13). Therefore, $Z^{\dagger}$ is the unique solution of the SDP (13).

Proof of Part (3): By the optimality and feasibility of the solutions $\tilde{Z}_{\rho_{1}}$ and $\tilde{Z}_{\rho_{2}}$, we have

[TABLE]

Adding these two inequalities together yields

[TABLE]

which implies $\operatorname{tr}(\tilde{Z}_{\rho_{1}})\geqslant\operatorname{tr}(\tilde{Z}_{\rho_{2}})$ when $\rho_{1}>\rho_{2}$. Moreover, if at least one of the SDPs has a unique solution, then at least one of the two inequalities in (50) is strict, implying

[TABLE]

and $\operatorname{tr}(\tilde{Z}_{\rho_{1}})>\operatorname{tr}(\tilde{Z}_{\rho_{2}})$.

Appendix A Proofs on spectral decompositions

Proof of Lemma 2.1.

Since $\kappa$ is symmetric and positive semidefinite, so is $\mathscr{R}$. Thus the corresponding operator $\mathscr{R}$ is self-adjoint in $L^{2}(\mathrm{d}\mu)$ and is also compact if (4) holds. Then $\mathscr{R}$ has a discrete set of nonnegative eigenvalues $\lambda_{0}\geqslant\lambda_{1}\geqslant\cdots\geqslant 0$, and has the following eigen-decomposition

[TABLE]

where $\{\phi_{j}\}_{j=0}^{\infty}$ is an orthonormal basis of $L^{2}(\mathrm{d}\mu)$. Denote $Z=\int_{S}d(y)\,\mathrm{d}\mu(y)$ so that $\pi(x)=d(x)/Z$. Since

[TABLE]

it follows that $\lambda_{0}=1$ is the largest eigenvalue of $\mathscr{R}$ with the eigenfunction $\phi_{0}(x)=\sqrt{\pi(x)}$ (of unit length) in $L^{2}(\mathrm{d}\mu)$. In addition, observe that

[TABLE]

This implies a decomposition of the transition probability $p(x,y)$ as

[TABLE]

where $\psi_{j}(x)=\phi_{j}(x)/\sqrt{d(x)}$ and $\widetilde{\psi}_{j}(y)=\phi_{j}(y)\sqrt{d(y)}$. In particular, for each $j=0,1,\ldots$,

[TABLE]

implying that $\{\psi_{j}\}_{j=0}^{\infty}$ are the corresponding (right) eigenfunctions of $\mathscr{P}$, with unit $L^{2}(d\mathrm{d}\mu)$ norm, associated with the same eigenvalues $1=\lambda_{0}\geqslant\lambda_{1}\geqslant\cdots\geqslant 0$. Since $\mathscr{P}$ is the transition operator of a Markov chain, $\lambda_{0}=1$ and $\psi_{0}\equiv Z^{-1/2}$. $\blacksquare$

Proof of Lemma 2.2.

For $t\in\mathbb{N}_{+}$ and $x,y\in S$, let $p_{t}(x,y)$ be the $t$-step transition probability from $x$ to $y$. By Lemma 2.1, we have

[TABLE]

and $\{\widetilde{\psi}_{j}\}_{j=0}^{\infty}$ forms an orthonormal basis of $L^{2}(\mathrm{d}\mu/d)$. Consequently, by viewing $\lambda_{j}^{t}\,\psi_{j}(x)$ as the coefficient associated with $\widetilde{\psi}_{j}$ in the orthogonal expansion of function $p_{t}(x,\cdot)$, we have

[TABLE]

$\blacksquare$

Appendix B Empirical diffusion affinity

Similar to the derivations in Appendix A, if we define matrix $R_{n}\in\mathbb{R}^{n\times n}$ with

[TABLE]

then $\{\lambda_{n,j}\}_{j=0}^{n-1}$ are also the eigenvalues of $R_{n}$. Let $\phi_{n,j}\in\mathbb{R}^{n}$ denote the unit Euclidean norm eigenvector associated with $\lambda_{n,j}$. Then the empirical probability transition matrix $P_{n}$ has the decomposition

[TABLE]

where $\psi_{n,j}=D_{n}^{-1/2}\,\phi_{n,j}\in\mathbb{R}^{n}$ and $\widetilde{\psi}_{n,j}=D_{n}^{1/2}\phi_{n,j}\in\mathbb{R}^{n}$, so that $\psi_{n,j}=D_{n}^{-1}\widetilde{\psi}_{n,j}$ for each $j\in\{0,1,\dots,n-1\}$. In particular, $\psi_{n,j}$ has unit $L^{2}(\mbox{diag}(D_{n}))$ norm, and $\widetilde{\psi}_{n,j}$ has unit $L^{2}(\mbox{diag}(D_{n}^{-1}))$ norm, for each $j\in\{0,1,\dots,n-1\}$.

In addition, we have the following relation between the diffusion affinity and $P_{n}^{2t}$,

[TABLE]

Appendix C Technical lemmas

In this appendix, we collect some technical lemmas used in the proofs of our main results.

Lemma C.1*.*

Let $Z^{*}$ be defined in (12). Then for any $Z\in\mathscr{C}$ defined in (13), we have

[TABLE]

In addition, if $Z$ also satisfies $\operatorname{tr}(Z)=\operatorname{tr}(Z^{\ast})$, or $Z\in\mathscr{C}_{K}$, where $\mathscr{C}_{K}$ is defined in (11), then

[TABLE]

Proof of Lemma C.1.

Inequalities (51) and (53) follow from Lemma 1 in [14]. Inequality (52) is due to inequality (57) in [3]. $\blacksquare$

Lemma C.2*.*

Let $Z^{*}$ be defined in (12). Then for any $Z\in\mathscr{C}$ defined in (13), we have

[TABLE]

and (51) holds for $Z\in\mathscr{C}$.

Proof of Lemma C.2.

The first inequality follows from inequality (57) in [3]. The second one follows from the first,  (51), and the following decomposition,

[TABLE]

with inequality $\|U\|_{1}\leqslant n{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|U\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\mbox{\scriptsize op}}\leqslant n{\left|\kern-1.07639pt\left|\kern-1.07639pt\left|U\right|\kern-1.07639pt\right|\kern-1.07639pt\right|}_{\ast}$ for any matrix $U\in\mathbb{R}^{n\times n}$. $\blacksquare$