Achieving Optimal Misclassification Proportion in Stochastic Block Model

Introduction

Network data analysis has become one of the leading topics in statistics. In fields such as physics, computer science, social science and biology, one observes a network among a large number of subjects of interest such as particles, computers, people, etc. The observed network can be modeled as an instance of a random graph and the goal is to infer structures of the underlying generating process. A structure of particular interest is community: there is a partition of the graph nodes in some suitable sense so that each node belongs to a community. Starting with the proposal of a series of methodologies , we have seen a tremendous literature devoted to algorithmic solutions to uncovering community structure and great advances have also been made in recent years on the theoretical understanding of the problem in terms of statistical consistency and thresholds for detection and exact recoveries. See, for instance, , among others. In spite of the great efforts exerted on this “community detection” problem, its state-of-the-art solution has not yet reached the comparable level of maturity as what statisticians have achieved in other high dimensional problems such as nonparametric estimation , high dimensional regression and covariance matrix estimation , etc. In these more well-established problems, not only do we know the fundamental statistical limits, we also have computationally feasible algorithms to achieve them. The major goal of the present paper is to serve as a step towards such maturity in network data analysis by proposing a computationally feasible algorithm for community detection in stochastic block model with provable statistical optimality.

To describe network data with community structure, we focus on the stochastic block model (SBM) proposed by . Let $A\in\left\{0,1\right\}^{n\times n}$ be the symmetric adjacency matrix of an undirected random graph generated according to an SBM with $k$ communities. Then the diagonal entries of $A$ are all zeros and each $A_{uv}=A_{vu}$ for $u>v$ is an independent Bernoulli random variable with mean $P_{uv}=B_{\sigma(u)\sigma(v)}$ for some symmetric connectivity matrix $B\in^{k\times k}$ and some label function $\sigma:[n]\rightarrow[k]$, where for any positive integer $m$, $[m]=\left\{1,\dots,m\right\}$. In other word, if the $u{{}^{\rm th}}$ node and the $v{{}^{\rm th}}$ node belong to the $i{{}^{\rm th}}$ and the $j{{}^{\rm th}}$ community respectively, then $\sigma(u)=i$, $\sigma(v)=j$ and there is an edge connecting $u$ and $v$ with probability $B_{ij}$. Community detection then refers to the problem of estimating the label function $\sigma$ subject to a permutation of the community labels $\left\{1,\dots,k\right\}$. A natural loss function for such an estimation problem is the proportion of wrong labels (subject to a permutation of the label set $[k]$), which we shall refer to as misclassification proportion from here on.

In ground breaking works by Mossel et al. and Massoulié , the authors established sharp threshold for the regimes in which it is possible and impossible to achieve a misclassification proportion strictly less than $\frac{1}{2}$ when $k=2$ and both communities are of the same size (so that it is better than random guess), which solved the conjecture in that was only justified in physics rigor. For some recent progress on the general case of fixed $k$ and possibly unequal sized communities, see . On the other hand, Abbe et al. , Mossel et al. and Hajek et al. established the necessary and sufficient condition for ensuring zero misclassification proportion (usually referred to as “strong consistency” in the literature) with high probability when $k=2$ and community sizes are equal, and was later generalized to a larger set of fixed $k$ by . Arguably, what is of more interest to statisticians is the intermediate regime between the above two cases, namely when the misclassification proportion is vanishing as the number of nodes grows but not exactly zero. This is usually called the regime of “weak consistency” in the network literature.

To achieve weak (and strong) consistency, statisticians have proposed various methods. One popular approach is spectral clustering which is motivated by the observation that the rank of the $n\times n$ matrix $P=(P_{uv})=(B_{\sigma(u)\sigma(v)})$ is at most $k$ and its leading eigenvectors contain information of the community structure. The application of spectral clustering on network data goes back to , and its performance under the stochastic block model has been investigated by , among others. To further improve the performance, various ways for refining spectral clustering have been proposed, such as those in which lead to strong consistency or convergence rates that are exponential in signal-to-noise ratio, while studied the problem of minimizing a non-vanishing misclassification proportion. However, in the regime of weak consistency, these refinement methods are not guaranteed to attain the optimal misclassification proportion to be introduced below. Another important line of research is devoted to the investigation of likelihood-based methods, which was initiated by and later extended to more general settings by . To tackle the intractability of optimizing the likelihood function, an EM algorithm using pseudo-likelihood was proposed by . Another way to overcome the intractability of the maximum likelihood estimator (MLE) is by convex relaxation. Various semi-definite relaxations were studied by , and the aforementioned sharp threshold for strong consistency in was indeed achieved by semi-definite programming. Recently, Zhang and Zhou established the minimax risk for misclassification proportion in SBM under weak conditions, which is of the form

if all $k$ communities are of equal sizes, where $I^{*}$ is the minimum Rényi divergence of order $\frac{1}{2}$ of the within and the between community edge distributions. See Theorem 1 below for a more general and precise statement of the minimax risk. Unfortunately, Zhang and Zhou used MLE for achieving the risk in (1) which was hence computationally intractable. Moreover, none of the spectral clustering based method or tractable variants of the likelihood based method has a known error bound that matches (1) with the sharp constant $1+o(1)$ on the exponent.

The main contribution of the current paper lies in the proposal of a computationally feasible algorithm that provably achieves the optimal misclassification proportion established in adaptively under weak regularity conditions. It covers the cases of both finite and diverging number of communities and both equal and unequal community sizes and achieves both weak and strong consistency in the respective regimes. In addition, the algorithm is guaranteed to compute in polynomial time even when the number of communities diverges with the number of nodes. Since the error bound of the algorithm matches the optimal misclassification proportion in under weak conditions, it achieves various existing detection boundaries in the literature. For instance, for any fixed number of communities, the procedure is weakly consistent under the necessary and sufficient condition of , and strongly consistent under the necessary and sufficient condition of . Moreover, it could match the optimal misclassification proportion in even when $k$ diverges. To the best of our limited knowledge, this is the first polynomial-time algorithm that achieves minimax optimal performance. In other words, the proposed procedure enjoys both statistical and computational efficiency.

The core of the algorithm is a refinement scheme for community detection motivated by penalized maximum likelihood estimation. As long as there exists an initial estimator that satisfies a certain weak consistency criterion, the refinement scheme is able to obtain an improved estimator that achieves the optimal misclassification proportion in (1) with high probability. The key to achieve the bound in (1) is to optimize the local penalized likelihood function for each node separately. This local optimization step is completely data-driven and has a closed form solution, and hence can be computed very efficiently. The additional penalty term is indispensable as it plays a key rule in ensuring the optimal performance when the community sizes are unequal and when the within community and/or between community edge probabilities are unequal.

To obtain a qualified initial estimator, we show that both spectral clustering and its normalized variant could satisfy the desired condition needed for subsequent refinement, though the refinement scheme works for any other method satisfying a certain weak consistency condition. Note that spectral clustering can be considered as a global method, and hence our two-stage algorithm runs in a “from global to local” fashion. In essence, with high probability, the global stage pinpoints a local neighborhood in which we shall search for solution to each local penalized maximum likelihood problem, and the subsequent local stage finds the desired solution. From this viewpoint, one can also regard our approach as an “optimization after localization” procedure. Historically, this idea played a key role in the development of the renowned one-step efficient estimator . It has also led to recent progress in non-convex optimization and localized gradient descent techniques for finding optimal solutions to high dimensional statistical problems. Examples include but are not limited to high-dimensional linear regression , sparse PCA , sparse CCA , phase retrieval and high dimensional EM algorithm . A closely related idea has also found success in the development of confidence intervals for regression coefficients in high dimensional linear regression. See, for instance, and the references therein. Last but not least, even when viewed as a “spectral clustering plus refinement” procedure, our method distinguishes itself from other such methods in the literature by provably achieving the minimax optimal performance over a wide range of parameter configurations.

The rest of the paper is organized as follows. Section 2 formally sets up the community detection problem and presents the two-stage algorithm. The theoretical guarantees for the proposed method are given in Section 3, followed by numerical results demonstrating its competitive performance on both simulated and real datasets in Sections 4 and 5. A discussion on the results in the current paper and possible directions for future investigation is included in Section 6. Section 7 presents the proofs of main results with some technical details deferred to the appendix.

Problem formulation and methodology

In this section, we give a precise formulation of the community detection problem and present a new method for it. The method consists of two stages: initialization and refinement. We shall first introduce the second stage, which is the main algorithm of the paper. It clusters the network data by performing a node-wise penalized neighbor voting based on some initial community assignment. Then, we will discuss several candidates for the initialization step including a new greedy algorithm for clustering the leading eigenvectors of the adjacency matrix or of the graph Laplacian that is tailored specifically for stochastic block model. Theoretical guarantees for the algorithms introduced in the current section will be presented in Section 3.

Recall that a stochastic block model is completely characterized by a symmetric connectivity matrix $B\in^{k\times k}$ and a label vector $\sigma\in[k]^{n}$. One widely studied parameter space of SBM is

where $\beta\geq 1$ is an absolute constant. This parameter space $\Theta_{0}(n,k,a,b,\beta)$ contains all SBMs in which the within community connection probabilities are all equal to $a\over n$ and the between community connection probabilities are all equal to $b\over n$. In the special case of $\beta=1$, all communities are of nearly equal sizes.

Assuming equal within and equal between connection probabilities can be restrictive. Thus, we also introduce the following larger parameter space

Throughout the paper, we treat $\beta\geq 1$ and $\alpha\geq 1$ as absolute constants, while $k,a,b$ and $\lambda$ should be viewed as functions of the number of nodes $n$ which can vary as $n$ grows. Moreover, we assume $0<\frac{b}{n}<\frac{a}{n}\leq 1-\epsilon$ throughout the paper for some numeric constant $\epsilon\in(0,1)$. Thus, the parameter space $\Theta(n,k,a,b,\lambda,\beta;\alpha)$ requires that the within community connection probabilities are bounded from below by $a\over n$ and the connection probabilities between any two communities are bounded from above by $b\over n$. In addition, it requires that the sizes of different communities are comparable. In order to guarantee that $\Theta(n,k,a,b,\lambda,\beta;\alpha)$ is a larger parameter space than $\Theta_{0}(n,k,a,b,\beta)$, we always require $\lambda$ to be positive and sufficiently small such that

According to Proposition 1 in the appendix, a sufficient condition for (4) is $\lambda\leq\frac{a-b}{2\beta k}$. We assume (4) throughout the rest of the paper.

The labels on the $n$ nodes induce a community structure $[n]=\cup_{i=1}^{k}\mathcal{C}_{i}$, where $\mathcal{C}_{i}=\left\{u\in[n]:\sigma(u)=i\right\}$ is the $i{{}^{\rm th}}$ community with size $n_{i}=|\mathcal{C}_{i}|$. Our goal is to reconstruct this partition, or equivalently, to estimate the label of each node modulo any permutation of label symbols. Therefore, a natural error measure is the misclassification proportion defined as

where $S_{k}$ stands for the symmetric group on $[k]$ consisting of all permutations of $[k]$.

2 Main algorithm

We are now ready to present the main method of the paper – a refinement algorithm for community detection in stochastic block model motivated by penalized local maximum likelihood estimation.

Indeed, for any SBM in the parameter space $\Theta_{0}(n,k,a,b,1)$ with equal community size, the MLE for $\sigma$ is

which is a combinatorial optimization problem and hence is computationally intractable. However, node-wise optimization of (6) has a simple closed form solution. Suppose the values of $\{\sigma(u)\}_{u=2}^{n}$ are known and we want to estimate $\sigma(1)$. Then, (6) reduces to

For each $i\in[k]$, the quantity $\sum_{\{v\neq 1:\sigma(v)=i\}}A_{1v}$ is exactly the number of neighbors that the first node has in the $i{{}^{\rm th}}$ community. Therefore, the most likely label for the first node is the one it has the most connections with when all communities are of equal sizes. In practice, we do not know any label in advance. However, we may estimate the labels of all but the first node by first applying a community detection algorithm $\sigma^{0}$ on the subnetwork excluding the first node and its associated edges, the adjacency matrix of which is denoted by $A_{-1}$ since it is the $(n-1)\times(n-1)$ submatrix of $A$ with its first row and first column removed. Once we estimate the remaining labels, we can apply (7) to estimate $\sigma(1)$ but with $\{\sigma(v)\}_{v=2}^{n}$ replaced with the estimated labels.

For any $u\in[n]$, let $A_{-u}$ denote the $(n-1)\times(n-1)$ submatrix of $A$ with its $u{{}^{\rm th}}$ row and $u{{}^{\rm th}}$ column removed. Given any community detection algorithm $\sigma^{0}$ which is able to cluster any graph on $n-1$ nodes into $k$ categories, we present the precise description of our refinement scheme in Algorithm 1.

The algorithm works in two consecutive steps. The first step carries out the foregoing heuristics on a node by node basis. For each fixed node $u$, we first leave the node out and apply the available community detection algorithm $\sigma^{0}$ on the remaining $n-1$ nodes and the edges among them (as summarized in the matrix $A_{-u}\in\left\{0,1\right\}^{(n-1)\times(n-1)}$) to obtain an initial community assignment vector $\sigma^{0}_{u}$. For convenience, we make $\sigma^{0}_{u}$ an $n$-vector by fixing $\sigma^{0}_{u}(u)=0$, though applying $\sigma^{0}$ on $A_{-u}$ does not give any community assignment for $u$. We then assign the label of the $u{{}^{\rm th}}$ node according to (10), which is essentially (7) with $\sigma$ replaced with $\sigma^{0}_{u}$ except for the additional penalty term. The additional penalty term is added to ensure the optimal performance even when both the diagonal and the off-diagonal entries of the connectivity matrix $B$ are allowed to take different values and the community sizes are not necessarily equal. To determine the penalty parameter $\rho_{u}$ in an adaptive way as spelled out in (11) – (12), we first estimate the connectivity matrix $B$ based on $A_{-u}$ in (8) – (9). After we obtain the community assignment for $u$, we organize the assignment for all $n$ vertices into an $n$-vector $\widehat{\sigma}_{u}$. We call this step “penalized neighbor voting” since the first term on the RHS of (10) counts the number of neighbors of $u$ in each (estimated) community while the second term is a penalty term proportional to the size of each (estimated) community.

Once we complete the above procedure for each of the $n$ nodes, we obtain $n$ vectors $\widehat{\sigma}_{u}\in[k]^{n}$, $u=1,\dots,n$, and turn to the second step of the algorithm. The basic idea behind the second step is to obtain a unified community assignment by assembling $\{\widehat{\sigma}_{u}(u):u\in[n]\}$ and the immediate hurdle is that each $\widehat{\sigma}_{u}$ is only determined up to a permutation of the community labels. Thus, the second step aims to find the right permutations by (13) before we assemble the $\widehat{\sigma}_{u}(u)$’s. We call this step “consensus” since we are essentially looking for a consensus on the community labels for $n$ possibly different community assignments, under the assumption that all of them are close to the ground truth up to some permutation.

3 Initialization via spectral methods

Another important issue in spectral clustering lies in the subsequent clustering method used to cluster the eigenvectors. A popular choice is $k$-means clustering. However, finding the global solution to the $k$-means problem is NP-hard . Kumar et al. proposed a polynomial time algorithm for achieving $(1+\epsilon)$ approximation to the $k$-means problem for any fixed $k$, which was utilized in to establish consistency for spectral clustering under stochastic block model with fixed number of communities. However, a closer look at the complexity bound suggests that the smallest possible $\epsilon$ is proportional to $k$. Thus, applying the algorithm and the associated bound in directly in our settings can lead to inferior error bounds when $k\to\infty$ as $n\to\infty$. To address this issue under stochastic block model, we propose a greedy clustering algorithm in Algorithm 2 inspired by the fact that the clustering centers under stochastic block model are well separated from each other on the population level. It is straightforward to check that the complexity of Algorithm 2 is polynomial in $n$.

Last but not least, we would like to emphasize that one needs not limit the initialization algorithm to the spectral methods introduced in this section. As Theorem 2 below shows, Algorithm 1 works for any initialization method that satisfies a weak consistency condition.

Theoretical properties

Before stating the theoretical properties of the proposed method, we first review the minimax rate in , which will be used as the optimality benchmark. The minimax risk is governed by the following critical quantity,

which is the Rényi divergence of order $\frac{1}{2}$ between $\mbox{Bern}\left(\frac{a}{n}\right)$ and $\mbox{Bern}\left(\frac{b}{n}\right)$, i.e., Bernoulli distributions with success probabilities $\frac{a}{n}$ and $\frac{b}{n}$ respectively. Recall that $0<\frac{b}{n}<\frac{a}{n}\leq 1-\epsilon$ is assumed throughout the paper. It can be shown that $I^{*}\asymp\frac{(a-b)^{2}}{na}$. Moreover, when $\frac{a}{n}=o(1)$,

where $H^{2}(P,Q)=\frac{1}{2}\int(\sqrt{{\rm d}P}-\sqrt{{\rm d}Q})^{2}$ is the squared Hellinger distance between two distributions $P$ and $Q$. The minimax rate for the parameter spaces (2) and (3) under the loss function (5) is given in the following theorem.

When $\frac{(a-b)^{2}}{ak\log k}\rightarrow\infty$, we have

for both $\Theta=\Theta_{0}(n,k,a,b,\beta)$ and $\Theta=\Theta(n,k,a,b,\lambda,\beta;\alpha)$ with any $\lambda\leq\frac{a-b}{2\beta k}$ and any $\beta\in[1,\sqrt{5/3})$, where $\eta=\eta_{n}\rightarrow 0$ is some sequence tending to as $n\rightarrow\infty$.

The assumption $\beta\in[1,\sqrt{5/3})$ is needed in for some technical reason. Here, the parameter $\beta$ enters the minimax rates when $k\geq 3$ since the worst case is essentially when one has two communities of size $\frac{n}{\beta k}$, while for $k=2$, the worst case is essentially two communities of size $\frac{n}{2}$. For all other results in this paper, we allow $\beta$ to be an arbitrary constant no less than $1$.

To this end, let us show that the two-stage algorithm proposed in Section 2 achieves the optimal misclassification proportion. The essence of the two-stage algorithm lies in the refinement scheme described in Algorithm 1. As long as any initialization step satisfies a certain weak consistency criterion, the refinement step directly leads to a solution with optimal misclassification proportion. To be specific, the initialization step needs to satisfy the following condition.

There exist constants $C_{0},\delta>0$ and a positive sequence $\gamma=\gamma_{n}$ such that

Under Condition 1, we have the following upper bounds regarding the performance of the proposed refinement scheme.

Suppose as $n\to\infty$, $\frac{(a-b)^{2}}{ak\log k}\to\infty$, $a\asymp b$ and Condition 1 is satisfied for

and $\Theta=\Theta_{0}(n,k,a,b,\beta)$. Then there is a sequence $\eta\to 0$ such that

If in addition Condition 1 is satisfied for $\gamma$ satisfying both (16) and

and $\Theta=\Theta(n,k,a,b,\lambda,\beta;\alpha)$, then the conclusion in (17) continues to hold for $\Theta=\Theta(n,k,a,b,\lambda,\beta;\alpha)$.

Theorem 2 assumes $a\asymp b$. The case when $a\asymp b$ may not hold is considered in Section 6. Compared with Theorem 1, the upper bounds (17) achieved by Algorithm 1 is minimax optimal. The condition (16) for the parameter space $\Theta_{0}(n,k,a,b,\beta)$ is very mild. When $k=O(1)$, it reduces to $\gamma=o(1)$ and simply means that the initialization should be weakly consistent at any rate. For $k\rightarrow\infty$, it implies that the misclassification proportion within each community converges to zero. Note that if the initialization step gives wrong labels to all nodes in one particular community, then the misclassification proportion is at least $1/k$. The condition (16) rules out this situation. For the parameter space $\Theta(n,k,a,b,\lambda,\beta;\alpha)$, an extra condition (18) is required. This is because estimating the connectivity matrix $B$ in $\Theta(n,k,a,b,\lambda,\beta;\alpha)$ is harder than in $\Theta_{0}(n,k,a,b,\beta)$. In other words, if we do not pursue adaptive estimation, (18) is not needed.

Theorem 2 is an adaptive result without assuming the knowledge of $a$ and $b$. When these two parameters are known, we can directly use $a$ and $b$ in (11) of Algorithm 1. By scrutinizing the proof of Theorem 2, the conditions (16) and (18) can be weakened as $\gamma=o(k^{-1})$ in this case.

Given the results of Theorem 2, it remains to check the initialization step via spectral clustering satisfies Condition 1. For matrix $P=(P_{uv})=(B_{\sigma(u)\sigma(v)})$ with $(B,\sigma)$ belonging to either $\Theta_{0}(n,k,a,b,\beta)$ or $\Theta(n,k,a,b,\lambda,\beta;\alpha)$, we use $\lambda_{k}$ to denote $\lambda_{k}(P)$. Define the average degree by

Assume $e\leq a\leq C_{1}b$ for some constant $C_{1}>0$ and

for some sufficiently small $c\in(0,1)$. Consider USC$(\tau)$ with a sufficiently small constant $\mu>0$ in Algorithm 2 and $\tau=C_{2}\bar{d}$ for some sufficiently large constant $C_{2}>0$. For any constant $C^{\prime}>0$, there exists some $C>0$ only depending on $C^{\prime},C_{1},C_{2}$ and $\mu$ such that

with probability at least $1-n^{-C^{\prime}}$. If $k$ is fixed, the same conclusion holds without assuming $a\leq C_{1}b$.

Theorem 3 improves the error bound for spectral clustering in . While requires the assumption $a>C\log n$, our result also holds for $a=o(\log n)$. A result close to ours is that by , but their clustering step is different from Algorithm 2. Moreover, the conclusion of Theorem 3 holds with probability $1-n^{-C^{\prime}}$ for an arbitrary large $C^{\prime}$, which is critical because the initialization step needs to satisfy Condition 1 for the subsequent refinement step to work. On the other hand, the bound in is stated with probability $1-o(1)$.

When $k=O(1)$, Theorem 2 and Theorem 3 jointly imply the following result.

Consider Algorithm 1 initialized by $\sigma^{0}$ with USC$(\tau)$ for $\tau=C\bar{d}$, where $C$ is a sufficiently large constant. Suppose as $n\rightarrow\infty$, $k=O(1)$, $\frac{(a-b)^{2}}{a}\rightarrow\infty$ and $a\asymp b$. Then, there exists a sequence $\eta\rightarrow 0$ such that

where the parameter space is $\Theta=\Theta_{0}(n,k,a,b,\beta)$.

Compared with Theorem 1, the proposed procedure achieves the minimax rate under the condition $\frac{(a-b)^{2}}{a}\rightarrow\infty$ and $a\asymp b$. When $k=O(1)$, the condition $\frac{(a-b)^{2}}{a}\rightarrow\infty$ is necessary and sufficient for weak consistency in view of Theorem 1. More general results including the case of $k\rightarrow\infty$ are stated and discussed in Section 6.

The following theorem characterizes the misclassification rate of normalized spectral clustering.

Assume $e\leq a\leq C_{1}b$ for some constant $C_{1}>0$ and

for some sufficiently small $c\in(0,1)$. Consider NSC$(\tau)$ with a sufficiently small constant $\mu>0$ in Algorithm 2 and $\tau=C_{2}\bar{d}$ for some sufficiently large constant $C_{2}>0$. Then, for any constant $C^{\prime}>0$, there exists some $C>0$ only depending on $C^{\prime},C_{1},C_{2}$ and $\mu$ such that

with probability at least $1-n^{-C^{\prime}}$. If $k$ is fixed, the same conclusion holds without assuming $a\leq C_{1}b$.

A slightly different regularization of normalized spectral clustering is studied by only for the dense regime, while Theorem 4 holds under both dense and sparse regimes. Moreover, our result also improves that of due to our tighter bound on $\|L(A_{\tau})-L(P_{\tau})\|_{\rm op}$ in Lemma 7 below. We conjecture that the $\log a$ factor in both the assumption and the bound of Theorem 4 can be removed.

Note that Theorem 3 and Theorem 4 are stated in terms of the quantity $\lambda_{k}$. We may specialize the results into the parameter spaces defined in (2) and (3). By Proposition 1, $\lambda_{k}\geq\frac{a-b}{2\beta k}$ for $\Theta_{0}(n,k,a,b,\beta)$ and $\lambda_{k}\geq\lambda$ for $\Theta(n,k,a,b,\lambda,\beta;\alpha)$. The implications of Theorem 3 and Theorem 4 and their uses as the initialization step for Algorithm 1 are discussed in full details in Section 6.

Numerical results

In this section we present the performance of the proposed algorithm on simulated datasets. The experiments cover three different scenarios: (1) dense network with communities of equal sizes; (2) dense network with communities of unequal sizes; and (3) sparse network. Recall the definition of $\bar{d}$ in (19). For each setting, we report results of Algorithm 1 initialized with four different approaches: USC$(\infty)$, USC$(2\bar{d})$, NSC$(0)$ and NSC$(\bar{d})$, the description of which can all be found in Section 2.3. For all these spectral clustering methods, Algorithm 2 was used to cluster the leading eigenvectors. The constant $\mu$ in the critical radius definition was set to be $0.5$ in all the results reported here. For each setting, the results are based on 100 independent draws from the underlying stochastic block model.

To achieve faster running time, we also ran a simplified version of Algorithm 1. Instead of obtaining $n$ different initializers $\{\sigma_{u}\}_{u\in[n]}$ to refine each node separately, the simplified algorithm refines all the nodes with a single initialization on the whole network. Thus, the running time can be reduced roughly by a factor of $n$. Simulation results below suggest that the simplified version achieves similar performances to that of Algorithm 1 in all the settings we have considered. For the precise description of the simplified algorithm, we refer readers to Algorithm 3 in the appendix.

In this setting, we generate networks with 2500 nodes and 10 communities, each of which consists of 250 nodes, and we set $B_{ii}=0.48$ for all $i$ and $B_{ij}=0.32$ for all $i\neq j$. Figure 1 shows the boxplots of the number of misclassified nodes. The first four boxplots correspond to the four different spectral clustering methods, in the order of USC$(\infty)$, USC$(2\bar{d})$, NSC$(0)$ and NSC$(\bar{d})$. The middle four correspond to the results achieved by applying the simplified refinement scheme with these four initialization methods, and the last four show the results of Algorithm 1 with these four initialization methods. Regardless of the initialization method, Algorithm 1 or its simplified version reduces the number of misclassified nodes from around 30 to around 5.

Imbalanced case

In this setting, we generate networks with 2000 nodes and 4 communities, the sizes of which are $200,400,600$ and $800$, respectively. The connectivity matrix is

Hence, the within-community edge probability is no smaller than 0.45 while the between-community edge probability is no greater than 0.35, and the underlying SBM is inhomogeneous. Figure 2 shows the boxplots of the number of misclassified nodes obtained by different initialization methods and their refinements, and the boxplots are presented in the same order as those in Figure 1. Similarly, we can see refinement significantly reduces the error.

Sparse case

In this setting we consider a much sparser stochastic block model than the previous two cases. In particular, each simulated network has 4000 nodes, divided into 10 communities all of size 400. We set all $B_{ii}=0.032$ and all $B_{ij}=0.005$ when $i\neq j$. The average degree of each node in the network is around 30. Figure 3 shows the boxplots of the number of misclassified nodes obtained by different initialization methods and their refinements, and the boxplots are presented in the same order as those in Figure 1. Compared with either USC or NSC initialization, refinement reduces the number of misclassified nodes by 50%.

Summary

In all three simulation settings, for all four initialization approaches considered, the refinement scheme in Algorithm 1 (and its simplified version) was able to significantly reduce the number of misclassified nodes, which is in agreement with the theoretical properties presented in Section 3.

Real data example

We now compare the results of our algorithm and some existing methods on a political blog dataset . Each node in this network represents a blog about US politics and a pair of nodes is connected if one blog contains a link to the other. There were 1490 nodes to start with, each labeled liberal or conservative. In what follows, we consider only the 1222 nodes located in the largest connected component of the network. This pre-processing step is the same as what was done in . After pre-processing, the network has 586 liberal blogs and 636 conservative ones which naturally form two communities. As shown in the right panel of Figure 4, nodes are more likely to be connected if they have the same political ideology.

Table 1 summarizes the results of Algorithm 1 and its simplified version on this network with four different initialization methods, as well as the performances of directly applying the four methods on the dataset. The average degree of the network $\bar{d}$ is 27, which is used as the tuning parameter for regularized NSC. For regularized USC, we set $\tau$ equals to twice the average degree, leading to the removal of 196 most connected nodes. The result of directly applying any of the four spectral clustering based initializations was unsatisfactory with at least 30% nodes misclassified. Despite the unsatisfactory performance of the initializers, Algorithm 1 and its simplified version are able to significantly reduce the number of misclassified nodes except for the case of NSC$(0)$, and the performance of the two are close to each other regardless of the initialization method.

An interesting observation is that if we apply the refinement scheme multiple times, the number of misclassified nodes keeps decreasing until convergence and the further reduction of misclassification proportion compared to a single refinement can be sizable. Figure 5 plots the numbers of misclassified nodes for multiple iterations of refinement via the simplified version of Algorithm 1. We are able to achieve 61, 58 or 63 misclassified nodes out of 1222 depending on which initialization method is used. For the three initialization methods included in the figure, the number of misclassified nodes converges within several iterations. NSC with $\tau=0$ is not included in Figure 5 due to the relatively inferior initialization, but its error also converges to around 60/1222 after 20 iterations. For comparison, state-of-the-art method such as SCORE was reported to achieve a comparable error of 58/1222. It is worth noting that SCORE was designed under the setting of degree-corrected stochastic block model, which fits the current dataset better than SBM due to the presence of hubs and low-degree nodes. The regularized spectral clustering implemented by , which was also designed under the degree-corrected stochastic block model, was reported to have an error of $(80\pm 2)/1222$. The semi-definite programming method by achieved 63/1222.

To summarize, our algorithm leads to significant performance improvement over several popular spectral clustering based methods on the political blog dataset. With repeated refinements, it demonstrates competitive performance even when compared with methods designed for models that better fit the current dataset.

Discussion

In this section, we discuss a few important issues related to the methodology and theory we have presented in the previous sections.

In Section 3, we established upper bounds on misclassification proportion under the assumption of $a\asymp b$. The following theorem shows that slightly weaker upper bounds can be obtained even when $a\asymp b$ does not hold. To state the result, recall that we assume throughout the paper $\frac{a}{n}\leq 1-\epsilon$ for some numeric constant $\epsilon\in(0,1)$.

Suppose as $n\to\infty$, $\frac{(a-b)^{2}}{ak\log k}\to\infty$ and Condition 1 is satisfied for $\gamma$ satisfying (16) and $\Theta=\Theta_{0}(n,k,a,b,\beta)$. Then for some positive constants $c_{\epsilon}$ and $C_{\epsilon}$ that depend only on $\epsilon$, for any sufficiently small constant $\epsilon_{0}\in(0,c_{\epsilon})$, if we replace the definition of $t_{u}$’s in (11) with

where $I^{*}$ is defined as in (14). In particular, we can set $C_{\epsilon}=\frac{10}{3}\frac{2-\epsilon}{\frac{\epsilon}{2}\log\frac{2}{\epsilon}}$ and $c_{\epsilon}=\min(\frac{1}{10C_{\epsilon}},\frac{\epsilon}{2-\epsilon})$.

If in addition Condition 1 is satisfied for $\gamma$ satisfying both (16) and (18) and $\Theta=\Theta(n,k,a,b,\lambda,\beta;\alpha)$, then the same conclusion holds for $\Theta=\Theta(n,k,a,b,\lambda,\beta;\alpha)$.

Compared with the conclusion (17) in Theorem 2, the vanish sequence $\eta$ in the exponent of the upper bound is replaced by $C_{\epsilon}\epsilon_{0}$, which is guaranteed to be smaller than $\min(0.1,\frac{2}{\log(2/\epsilon)})$ and can be driven to be arbitrarily small by decreasing $\epsilon_{0}$. To achieve this, the $t_{u}$’s used in defining the penalty parameters in the penalized neighbor voting step need to be truncated at the value $\log\frac{1}{\epsilon_{0}/2}$.

2 Implications of the results

We now discuss some implications of the results in Theorems 2 – 5.

When using USC as initialization for Algorithm 1, we obtain the following results by combining Theorem 2, Theorem 3 and Theorem 5. Recall that $\bar{d}$ is the average degree of nodes in $A$ defined in (19).

Consider Algorithm 1 initialized by $\sigma^{0}$ with USC$(\tau)$ with $\tau=C\bar{d}$ for some sufficiently large constant $C>0$. If as $n\rightarrow\infty$, $a\asymp b$ and

then there is a sequence $\eta\rightarrow 0$ such that (17) holds with $\Theta=\Theta_{0}(n,k,a,b,\beta)$. If as $n\rightarrow\infty$, $a\asymp b$ and

then (17) holds for $\Theta=\Theta(n,k,a,b,\lambda,\beta;\alpha)$. If for either parameter space, $a\asymp b$ may not hold but $k$ is fixed and (24) or (25) holds respectively, then (23) holds as long as $t_{u}$ is replaced by (22) in Algorithm 1.

regardless of the behavior of $\frac{a}{b}$. On the other hand, Theorem 1 shows that it is impossible to achieve strong consistency if

When $\frac{a}{n}=o(1)$, $nI^{*}=(1+o(1))(\sqrt{a}-\sqrt{b})^{2}$ and so one can replace ${nI^{*}}$ in (26) – (27) with ${(\sqrt{a}-\sqrt{b})^{2}}$. In the literature, Abbe et al. , Mossel et al. and Hajek et al. obtained comparable strong consistency results via efficient algorithms for the special case of two communities of equal sizes, i.e., $k=2$ and $\beta=1$. Under the additional assumption of $a\asymp b\asymp{\log n}$, Hajek et al. later achieved the result via efficient algorithm for the case of fixed $k$ and $\beta=1$, and Abbe and Sandon investigated the case of fixed $k$ and $\beta\geq 1$. In comparison, our result holds for any fixed $k$ and any $\beta\geq 1$ without assuming $a\asymp b\asymp\log n$. In the weak consistency regime, in terms of misclassification proportion, for the special case of $k=2$ and $\beta=1$, Yun and Proutiere achieved the optimal rate for $\Theta_{0}(n,2,a,b,1)$ when $a\asymp b\asymp a-b$, while the error bounds in other papers are typically off by a constant multiplier on the exponent. In comparison, Theorem 6 provides optimal results (17) and near optimal results (23) for a much broader class of models under much weaker conditions. Last but not least, our algorithm can provably achieve strong consistency and minimax optimal performance even for growing $k$, which to our limited knowledge, is the first in the literature.

The performance of Algorithm 1 initialized by NSC can be summarized as the following theorem by combining Theorem 2, Theorem 4 and Theorem 5. In this case, the sufficient condition for achieving minimax optimal performance is slightly stronger than when USC is used for initialization.

Consider Algorithm 1 initialized by $\sigma^{0}$ with NSC$(\tau)$ with $\tau=C\bar{d}$ for some sufficiently large constant $C>0$. If as $n\rightarrow\infty$, $a\asymp b$ and

then there is a sequence $\eta\rightarrow 0$ such that (17) holds with $\Theta=\Theta_{0}(n,k,a,b,\beta)$. If as $n\rightarrow\infty$, $a\asymp b$ and

then (17) holds for $\Theta=\Theta(n,k,a,b,\lambda,\beta;\alpha)$. If for either parameter space, $a\asymp b$ may not hold but $k$ is fixed and (28) or (29) holds respectively, then (23) holds as long as $t_{u}$ is replaced by (22) in Algorithm 1.

Last but not least, we would like to point out that when the key parameters $a$ and $b$ are known, we can obtain the desired performance guarantee under weaker conditions as summarized in the following theorem.

Suppose $a,b$ are known. Consider Algorithm 1 initialized by $\sigma^{0}$ with USC($\tau$) with $\tau=Ca$ for some sufficiently large constant $C>0$ and $\widehat{a}_{u}=a$, $\widehat{b}_{u}=b$ in (9) for all $u\in[n]$. If as $n\to\infty$, $a\asymp b$ and

then there is a sequence $\eta\to 0$ such that (17) holds with $\Theta=\Theta_{0}(n,k,a,b,\beta)$. If as $n\to\infty$, $a\asymp b$ and

then (17) holds with $\Theta=\Theta(n,k,a,b,\lambda,\beta;\alpha)$. If for either parameter space without assuming $a\asymp b$, (30) or (31) holds respectively, then (23) holds if in addition $t_{u}$ is replaced by (22).

If instead NSC($\tau$) is used for initialization with $\tau=Ca$ for some sufficiently large constant $C>0$, then the above conclusions hold if we replace (30) with $\frac{(a-b)^{2}}{ak^{3}\log a}\to\infty$ and (31) with $\frac{\lambda^{2}}{ak\log a}\to\infty$, respectively.

3 Potential future research problems

In simulation studies, we experimented a simplified version of Algorithm 1 (with precise description as Algorithm 3 in appendix) and showed that it provided similar performance to Algorithm 1 on simulated datasets. Moreover, for the political blog data, we showed that iterative application of this simplified refinement scheme kept driving down the number of misclassified nodes till convergence. It is of great interest to see if comparable theoretical results to Theorem 2 could be established for the simplified and/or the iterative version, and if the iterative version converges to a local optimum of certain objective function for community detection. Though answering these intriguing questions is beyond the scope of the current paper, we think it can serve as an interesting future research problem.

Data-driven choice of k𝑘k

The knowledge of $k$ is assumed and is used in both methodology and theory of the present paper. Date-driven choice of $k$ is of both practical importance and contemporary research interest, and researchers have proposed various ways to achieve this goal for stochastic block model, including cross-validation , Tracy–Widom test , information criterion , likelihood ratio test , etc. Whether these methods are optimal and whether it is possible to select $k$ in a statistically optimal way remains an important open problem.

More general models

The results in this paper cover a large range of parameter spaces for stochastic block models and we show the competitive performance of the proposed algorithm both in theory and on numerical examples. Despite its popularity, stochastic block model has its own limits for modeling network data. Therefore, an important future research direction is to design computationally feasible algorithms that can achieve statistically optimal performance for more general network models, such as degree-corrected stochastic block models.

Proofs of main results

The main result of the paper, Theorem 2, is proved in Section 7.1. Theorem 3 and Theorem 4 are proved in Section 7.2 and Section 7.3 respectively. The proofs of the remaining results, together with some auxiliary lemmas, are given in the appendix.

We first state a lemma that guarantees the accuracy of parameter estimation in Algorithm 1.

Let $\Theta=\Theta(n,k,a,b,\lambda,\beta;\alpha)$. Suppose as $n\to\infty$, $\frac{(a-b)^{2}}{ak}\to\infty$ and Condition 1 holds with $\gamma$ satisfying (16) and (18). Then there is a sequence $\eta\to 0$ as $n\to\infty$ and a constant $C>0$ such that

For $\Theta=\Theta_{0}(n,k,a,b,\beta)$, the conclusion (32) continues to hold even when the assumption (18) is dropped.

$1^{\circ}$ Let $\Theta=\Theta(n,k,a,b,\lambda,\beta;\alpha)$. For any community assignments $\sigma_{1}$ and $\sigma_{2}$, define

Fix any $(B,\sigma)\in\Theta$ and $u\in[n]$. Define event

Let ${\mathcal{C}}_{i}^{\prime}$ be any deterministic subset of $[n]$ such that (35) holds with $\widetilde{\mathcal{C}}_{i}^{u}$ replaced by ${\mathcal{C}}_{i}^{\prime}$. By definition, there are at most

different subsets with this property where $C_{1}>0$ is an absolute constant. Let ${\mathcal{E}}_{i}^{\prime}$ be the edges within ${\mathcal{C}}_{i}^{\prime}$. Then $|{\mathcal{E}}_{i}^{\prime}|$ consists of independent Bernoulli random variables, where at least $(1-\beta\gamma k)^{2}$ proportion of them follow the Bern$(B_{ii})$ distribution, at most $(\beta\gamma k)^{2}$ proportion that are stochastically smaller than Bern$(\frac{\alpha a}{n})$ and stochastically larger than Bern$(\frac{a}{n})$, and at most $2\beta\gamma k$ proportion are stochastically smaller than Bern$(\frac{b}{n})$. Therefore, we obtain that

Note that the LHS is $(1-(2+o(1))\beta\gamma k)B_{ii}$. On the other hand, under condition (18), the RHS is attained at $t=0$ and equals $B_{ii}$ exactly. Thus, we conclude that

for some $\eta^{\prime}\to 0$ that depends only on $a,k,\alpha,\beta$ and $\gamma$, where the last inequality is due to (18).

On the other hand, by Bernstein’s inequality, for any $t>0$,

where we the second inequality holds since $\frac{\log x}{x}$ is monotone decreasing as $x$ increases and so $\gamma\log\gamma^{-1}\geq\frac{1}{n}\log n$ for any $\gamma\geq\frac{1}{n}$, which is the case of most interest since $\gamma<\frac{1}{n}$ leads to $\gamma=0$ and so the initialization is already perfect. Even when $\gamma=0$, we can still continue to the following arguments by replacing every $\gamma$ with $\frac{1}{n}$ and all the steps continue to hold. Thus, we obtain that for positive constant $C_{\alpha,\beta,\delta}$ that depends only on $\alpha,\beta$ and $\delta$,

Thus, with probability at least $1-\exp\left\{-C_{1}\gamma n\log\gamma^{-1}\right\}n^{-(3+\delta)}$,

where $\eta^{\prime}\to 0$ depends only on $a,k,\alpha,\beta,\gamma$ and $\delta$. Here, the last inequality holds since

where $\sqrt{ak}\ll a-b$ since $\frac{(a-b)^{2}}{ak}\to\infty$ and $k\gamma\log\gamma^{-1}=O(1)$, and

We combine (37) and (39) and apply the union bound to obtain that for a sequence $\eta\to 0$ that depends only on $a,k,\alpha,\beta,\gamma$ and $\delta$, with probability at least $1-n^{-(3+\delta)}$

The proof for $B_{ij}$ estimation is analogous and hence is omitted. A final union bound on $i,j\in[k]$ leads to the desired claim since all the constants and vanishing sequences in the above analysis depend only on $a,b,k,\alpha,\beta,\gamma$ and $\delta$, but not on $u$, $B$ or $\sigma$.

$2^{\circ}$ If $\Theta=\Theta_{0}(n,k,a,b,\beta)$, then condition (18) on $\gamma$ is no longer needed. This is because (36) can be replaced by

where the LHS equals $\frac{a}{n}-(1-\beta\gamma k(1+o(1)))\frac{a-b}{n}=\frac{a}{n}+o(\frac{a-b}{n})$ and the RHS equals $\frac{a}{n}$. Thus, no additional condition is needed to guarantee (37) in the foregoing arguments. This completes the proof. ∎

The next two lemmas establish the desired error bound for the node-wise refinement.

Let $\Theta_{0}$ be defined as in (2) and $k\geq 2$. Suppose as $n\to\infty$, $\frac{(a-b)^{2}}{ak}\to\infty$ and $a\asymp b$. If there exists two sequences $\gamma=o(1/k)$ and $\eta=o(1)$, constants $C,\delta>0$ and permutations $\left\{\pi_{u}\right\}_{u=1}^{n}\subset S_{k}$ such that

Then for $\widehat{\sigma}_{u}(u)$ defined as in (10) with $\rho=\rho_{u}$ in (12), there is a sequence $\eta^{\prime}=o(1)$ such that for $k=2$,

In what follows, let $E_{u}$ denote the event in (41). For the sake of brevity, we let $p=a/n$, $q=b/n$, $\widehat{p}_{u}=\widehat{a}_{u}/n$ and $\widehat{q}_{u}=\widehat{b}_{u}/n$. Moreover, let $\sigma_{u}=\pi_{u}(\sigma)$, $n_{i}=|\{v:\sigma_{u}(v)=i\}|$, $m_{i}=|\{v:\sigma^{0}_{u}(v)=i\}|$ and $m_{i}^{\prime}=|\{v:\sigma^{0}_{u}(v)=\sigma_{u}(v)=i\}|$. Without loss of generality, let $\sigma_{u}(u)=1$.

Now we bound each $p_{l}$. By the independence structure and Chernoff bound, we have

We are going to give bounds for the terms in (45) and (45) respectively. Before doing that, we need some preparatory inequalities. Define $t^{*}$ through the equation

for some constant $C_{2}>0$. Therefore, for the term in (45), on the event $E_{u}$,

By (46), the term in (49) is upper bounded by

By (47), the term in (49) is upper bounded by

Now we provide an upper bound for (45). By (47), on $E_{u}$,

When $k=2$, $\min_{l\neq 1}\left(\frac{n_{1}+n_{l}}{2}\right)=\frac{n}{2}$, and when $k\geq 3$, $\min_{l\neq 1}\left(\frac{n_{1}+n_{l}}{2}\right)\geq\frac{n}{\beta k}$. Thus, the proof is complete. ∎

Let $\Theta$ be defined as in (3) and $k\geq 2$. Suppose as $n\to\infty$, $\frac{(a-b)^{2}}{ak}\to\infty$ and $a\asymp b$. If there exists two sequences $\gamma=o\left(\frac{a-b}{ak}\right)$ and $\eta=o(1)$, constants $C,\delta>0$ and permutations $\left\{\pi_{u}\right\}_{u=1}^{n}\subset S_{k}$ such that (41) holds. Then for $\widehat{\sigma}_{u}(u)$ defined as in (10) with $\rho=\rho_{u}$ in (12), the conclusions of Lemma 2 continue to hold.

The proof is similar to that of Lemma 2 and we use the same notation as there. First, we give a bound for $p_{l}$ defined in (42). Let $X_{j}\sim\mbox{Bern}(q)$, $Y_{j}\sim\mbox{Bern}(p)$ and $Z_{j}\sim\mbox{Bern}(\alpha p)$, $j\geq 1$, be mutually independent. Then, a stochastic order argument gives

Note that the term in (54) is the same as that in (45), and thus it can be upper bounded by (50) as before. To bound for (54), observe that by (47),

Thus, under the assumption $\gamma=o\left(\frac{p-q}{kp}\right)$, the term (54) is bounded by $\exp\left(o(1)\frac{n_{1}+n_{l}}{2}I^{*}\right)$. The remaining proof is the same as that of Lemma 2. ∎

Finally, we need a lemma to justify the consensus step in Algorithm 1.

For any community assignments $\sigma$ and $\sigma^{\prime}$: $[n]\to[k]$, such that for some constant $C\geq 1$

Thus, what remains to be shown is that $\xi\in S_{k}$, i.e., $\xi(l_{1})\neq\xi(l_{2})$ for any $l_{1}\neq l_{2}$. To this end, note that if for some $l_{1}\neq l_{2}$, $\xi(l_{1})=\xi(l_{2})$, then there would exist some $l_{0}\in[k]$ such that for any $l\in[k]$, $\xi(l)\neq l_{0}$, and so

This is in contradiction to the second last display, and hence $\xi\in S_{k}$. This completes the proof. ∎

Let $\Theta=\Theta(n,k,a,b,\lambda,\beta;\alpha)$, and fix any $(B,\sigma)\in\Theta$. For any $u\in[n]$, by Condition 1 and the fact that $\sigma^{0}_{u}$ and $\widehat{\sigma}_{u}$ differ only at the community assignment of $u$, for $\gamma^{\prime}=\gamma+1/n$, there exists some $\pi_{u}\in S_{k}$ such that

In addition, (56) implies with probability at least $1-Cn^{-(1+\delta)}$, we have

When $k\geq 3$, Lemma 1, (16) and (18) imply that the condition of Lemma 3 is satisfied, which in turn implies that for a sequence $\eta^{\prime}=o(1)$,

where the last inequality holds since $\frac{\,nI^{*}}{k}\asymp\frac{(a-b)^{2}}{ak}\to\infty$. Thus, Markov’s inequality leads to

If $(k-1)\exp\left\{-(1-\eta)\frac{nI^{*}}{\beta k}\right\}\geq n^{-(1+\delta/2)}$, then

If $(k-1)\exp\left\{-(1-\eta)\frac{nI^{*}}{\beta k}\right\}<n^{-(1+\delta/2)}$, then

Here, the second last inequality holds since $\eta>\eta^{\prime}$ and so $(k-1)\exp\left\{-(1-\eta^{\prime})nI^{*}/(\beta k)\right\}<(k-1)\exp\left\{-(1-\eta)nI^{*}/(\beta k)\right\}<n^{-(1+\delta/2)}$. We complete the proof for the case of $\Theta(n,k,a,b,\lambda,\beta;\alpha)$ and $k\geq 3$ by noting that $(k-1)\exp\left\{-(1-\eta)\frac{nI^{*}}{\beta k}\right\}=\exp\left\{-(1-\eta^{\prime\prime})\frac{nI^{*}}{\beta k}\right\}$ for another sequence $\eta^{\prime\prime}=o(1)$ under the assumption $\frac{(a-b)^{2}}{ak\log k}\rightarrow\infty$ and no constant or sequence in the foregoing arguments involves $B,\sigma$ or $u$. When $\Theta=\Theta(n,k,a,b,\lambda,\beta;\alpha)$ and $k=2$, the foregoing arguments continue to hold with $\beta$ and $k$ replaced with $1$ and $2$ respectively.

When $\Theta=\Theta_{0}(n,k,a,b,\beta)$, we can run the foregoing arguments with Lemma 3 replaced by Lemma 2 to reach the conclusion in (17), which does not require condition (18). This completes the proof. ∎

2 Proof of Theorem 3

Consider a symmetric adjacency matrix $A\in\{0,1\}^{n\times n}$ and a symmetric matrix $P\in^{n\times n}$ satisfying $A_{uu}=0$ for all $u\in[n]$ and $A_{uv}\sim\text{Bernoulli}(P_{uv})$ independently for all $u>v$. For any $C^{\prime}>0$, there exists some $C>0$ such that

with probability at least $1-n^{-C^{\prime}}$ uniformly over $\tau\in[C_{1}(np_{\max}+1),C_{2}(np_{\max}+1)]$ for some sufficiently large constants $C_{1},C_{2}$, where $p_{\max}=\max_{u\geq v}P_{uv}$.

For $P=(P_{uv})=(B_{\sigma(u)\sigma(v)})$, we have SVD $P=U\Lambda U^{T}$, where

with $\Delta=\text{diag}(\sqrt{n_{1}},...,\sqrt{n_{k}})$, $Z\in\{0,1\}^{n\times k}$ is a matrix with exactly one nonzero entry in each row at $(i,\sigma(i))$ taking value $1$ and $W\in O(k,k)$.

and observe that $Z\Delta^{-1}\in O(n,k)$. Apply SVD to the matrix $\Delta B\Delta^{T}=W\Lambda W^{T}$ for some $W\in O(k,k)$, and then we have $P=U\Lambda U^{T}$ with $U=Z\Delta^{-1}W\in O(k,k)$. ∎

where $V=Z\Delta^{-1}W_{2}\in O(n,k)$ for some $W_{2}\in O(k,k)$. Combining (60), Lemma 5 and the conclusion $\tau\in[C_{1}a,C_{2}a]$, we have

with probability at least $1-n^{-C^{\prime}}$. The definition of $V$ implies that

By definition, $T_{i}\cap T_{j}=\varnothing$ when $i\neq j$, and we also have

where the last inequality is by (61). After rearrangement, we have

In other words, most nodes are close to the centers and are in the set (63). Note that the sets $\{T_{i}\}_{i\in[k]}$ are disjoint. Suppose there is some $i\in[k]$ such that $|T_{i}|<|\sigma^{-1}(i)|-\left|\left(\cup_{i\in[k]}T_{i}\right)^{c}\right|$, we have $\left|\cup_{i\in[k]}T_{i}\right|=\sum_{i\in[k]}|T_{i}|<n-\left|\left(\cup_{i\in[k]}T_{i}\right)^{c}\right|=\left|\cup_{i\in[k]}T_{i}\right|$, which is impossible. Thus, the cardinality of $T_{i}$ for each $i\in[k]$ is lower bounded as

where the last inequality above is by the assumption (20). Intuitively speaking, except for a negligible proportion, most data points in $\{\widehat{U}_{u*}\}_{u\in[n]}$ are very close to the population centers $\{Q_{i*}\}_{i\in[k]}$. Since the centers are at least $\sqrt{\frac{2k}{\beta n}}$ away from each other and $\{T_{i}\}_{i\in[k]}$ and $\{\widehat{C}_{i}\}_{i\in[k]}$ are both defined through the critical radius $r=\mu\sqrt{\frac{k}{n}}$ for a small $\mu$, each $\widehat{C}_{i}$ should intersect with only one $T_{i}$ (see Figure 6). We claim that there exists some permutation $\pi$ of the set $[k]$, such that for $\widehat{C}_{i}$ defined in Algorithm 2,

In what follows, we first establish the result of Theorem 3 by assuming (66). The proof of (66) will be given in the end. Note that for any $i\neq j$, $T_{\pi(i)}\cap\widehat{\mathcal{C}}_{j}=\varnothing$, which is deduced from the fact that $\widehat{\mathcal{C}}_{j}\cap T_{\pi(j)}\neq\varnothing$ and the definition of $\widehat{\mathcal{C}}_{j}$. Therefore, $T_{\pi(i)}\subset\widehat{\mathcal{C}}_{j}^{c}$ for all $j\neq i$. Combining with the fact that $T_{\pi(i)}\cap\widehat{\mathcal{C}}_{i}^{c}\subset\widehat{\mathcal{C}}_{i}^{c}$, we get $T_{\pi(i)}\cap\widehat{\mathcal{C}}_{i}^{c}\subset(\cup_{i\in[k]}\widehat{\mathcal{C}}_{i})^{c}$. Therefore,

Since $T_{i}\cap T_{j}=\varnothing$ for $i\neq j$, we deduce from (67) that

By definition, $\widehat{\mathcal{C}}_{i}\cap\widehat{\mathcal{C}}_{j}=\varnothing$ for $i\neq j$, we deduce from (66) that