Compressive Spectral Clustering

Introduction

SC is mainly used in two contexts: $1)$ if the $N$ data points show particular structures (e.g., concentric circles) for which naive $k$-means clustering fails; $2$) if the input data is directly a graph $\mathcal{G}$ modeling a network (White & Smyth, 2005), such as social, neuronal, or transportation networks. SC suffers nevertheless from three main computational bottlenecks for large $N$ and/or $k$: the creation of the similarity matrix $\mathsf{W}$; the partial eigendecomposition of the graph Laplacian matrix $\mathsf{L}$; and $k$-means.

Circumventing these bottlenecks has raised a significant interest in the past decade. Several authors have proposed ideas to tackle the eigendecomposition bottleneck, e.g., via the power method (Boutsidis & Gittens, 2015; Lin & Cohen, 2010), via a careful optimisation of diagonalisation algorithms in the context of SC (Liu et al., 2007), or via matrix column-subsampling such as in the Nyström method (Fowlkes et al., 2004), the nSPEC and cSPEC methods of (Wang et al., 2009), or in (Chen & Cai, 2011; Sakai & Imiya, 2009). All these methods aim to quickly compute feature vectors, but $k$-means is still applied on $N$ feature vectors. Other authors, inspired by research aiming at reducing k-means complexity (Jain, 2010), such as the line of work on coresets (Har-Peled & Mazumdar, 2004), have proposed to circumvent $k$-means in high dimension by subsampling a few data points out of the $N$ available ones, applying SC on its reduced similarity graph, and interpolating the results back on the complete dataset. One can find similar methods in (Yan et al., 2009) and (Wang et al., 2009)’s eSPEC proposition, where two different interpolation methods are used. Both methods are heuristic: there is no proof that these methods approach the results of SC. Also, let us mention (Dhillon et al., 2007) that circumvents both the eigendecomposition and the $k$-means bottlenecks: the authors reduce the graph’s size by successive aggregation of nodes, apply SC on this small graph, and propagate the results on the complete graph using kernel $k$-means to control interpolation errors. The kernel is computed so that kernel $k$-means and SC share the same objective function (Filippone et al., 2008). Finally, we mention works (Boutsidis et al., 2011; Cohen et al., 2015) that concentrate on reducing the feature vectors’ dimension in the $k$-means problem, but do not sidestep the eigendecomposition nor the large $N$ issues.

2 Contribution: compressive clustering

The first ingredient builds upon recent works (Tremblay et al., 2016; Ramasamy & Madhow, 2015) that avoid the costly computation of the eigenvectors of $\mathsf{L}$ by filtering $O(\log(k))$ random signals on $\mathcal{G}$ that will then serve as feature vectors to perform clustering. We show in this paper how to incorporate the effects of non-ideal, but computationally efficient, graph filters on the quality of the feature vectors used for clustering.

The second ingredient uses a recent sampling theory of bandlimited graph-signals (Puy et al., 2015) to reduce the computational complexity of $k$-means. Using the fact that the indicator vectors of each cluster are approximately bandlimited on $\mathcal{G}$, we prove that clustering a random subset of $O(k\log(k))$ nodes of $\mathcal{G}$ using random features vectors of size $O(\log(k))$ is sufficient to infer rapidly and accurately the cluster label of all $N$ nodes of the graph. Note that the complexity of $k$-means is reduced to $O(k^{2}\log^{2}(k))$ instead of $O(Nk^{2})$ for SC. One readily sees that this method scales easily to large datasets, as will be demonstrated on artifical and real-world datasets containing up to $N=10^{6}$ nodes.

The proposed compressive spectral clustering method can be summarised as follows:

generate a feature vector for each node by filtering $O(\log(k))$ random Gaussian signals on $\mathcal{G}$;

sample $O(k\log(k))$ nodes from the full set of nodes;

interpolate the cluster indicator vectors back to the complete graph.

Background

Let $\mathcal{G}=(\mathcal{V},\mathcal{E},\mathbf{W})$ be an undirected weighted graph with $\mathcal{V}$ the set of $N$ nodes, $\mathcal{E}$ the set of edges, and $\mathbf{W}$ the weighted adjacency matrix such that $W_{ij}=W_{ji}\geqslant 0$ is the weight of the edge between nodes $i$ and $j$.

Denote by $\mathsf{H}_{\lambda_{c}}$ the graph filter operator associated to $h_{\lambda_{c}}$.

2 Spectral clustering

We choose here Ng et al.’s method (2002) based on the normalized Laplacian as our standard SC method. The input is the adjacency matrix $\mathsf{W}$ representing the pairwise similarity of all the $N$ objects to clusterIn network analysis, the raw data is directly $\mathsf{W}$. In the case where one starts with a set of data points $(\bm{x}_{1},\ldots,\bm{x}_{N})$, the first step consists in deriving $\mathsf{W}$ from the pairwise similarities $s(\bm{x}_{i},\bm{x}_{j})$. See (von Luxburg, 2007) for several choices of similarity measure $s$ and several ways to create $\mathsf{W}$ from the $s(\bm{x}_{i},\bm{x}_{j})$.. After computing its Laplacian $\mathsf{L}$, follow Alg. 1 to find $k$ classes.

Principles of CSC

Compressive spectral clustering (CSC) circumvents two of SC’s bottlenecks, the partial diagonalisation of the Laplacian and the high-dimensional $k$-means, thanks to the following ideas.

2) Run $k$-means on $n$ randomly selected feature vectors out of the $N$ available ones - thus clustering the corresponding $n$ nodes into $k$ groups - and interpolate the result back on the full graph. To guarantee robust reconstruction, we take advantage of our recent results on random sampling of $k$-bandlimited graph signals. In Sec. 3.2, we explain why these results are applicable to clustering and show that it is sufficient to sample $n=O(k\log{k})$ features only! Note that to cluster data into $k$ groups, one needs at least $k$ samples. This result is thus optimal up to the extra $\log{k}$ factor.

The following theorem shows that, for large enough $d$,

is a good estimation of $D_{ij}$ with high probability.

Let $\epsilon\in]0,1]$ and $\beta>0$ be given. If $d$ is larger than

then with probability at least $1-N^{-\beta}$, we have

The proof is provided in the supplementary material.

In Sec. 4.2, we generalize this result to the real-world case where the low-pass filter is approximated by a finite order polynomial; we also prove that, as announced in the introduction, one only needs $d=O(\log{k})$ features when using the downsampling scheme that we now detail.

2 Downsampling and interpolation

For a simple regular (with nodes of same degree) graph of $k$ disconnected clusters, it is easy to check that $\{\bm{c}_{1},\ldots,\bm{c}_{k}\}$ form a set of orthogonal eigenvectors of $\mathsf{L}$ with eigenvalue . All indicator vectors $\bm{c}_{j}$ therefore live in ${\rm span}{(\mathsf{U}_{k})}$. For general graphs, we assume that the indicator vectors $\bm{c}_{j}$ live close to ${\rm span}{(\mathsf{U}_{k})}$, i.e., the difference between any $\bm{c}_{j}$ and its orthogonal projection onto ${\rm span}(\mathsf{U}_{k})$ is small. Experiments in Section 5 will confirm that it is a good enough model to recover the cluster indicator vectors.

In graph signal processing words, one can say that $\bm{c}_{j}$ is approximately $k$-bandlimited, i.e., its $k$ first graph Fourier coefficients bear most of its energy. There has been recently a surge of interest around adapting classical sampling theorems to such bandlimited graph signals (Chen et al., 2015; Anis et al., 2015; Tsitsvero et al., 2015; Marques et al., 2015). We rely here on the random sampling strategy proposed in (Puy et al., 2015) to select a subset of $n$ nodes.

2.2 Sampling and interpolation

To recover $\bm{c}_{j}$ from its $n$ observations $\bm{c}_{j}^{r}$, Puy et al. (2015) show that the solution to the optimisation problem

2.3 How many features to sample?

We terminate this section by providing the theoretical number of features $n$ one needs to sample in order to make sure that the indicator vectors can be faithfully recovered. This number is driven by the following quantity.

The global cumulative coherence of order $k$ of the graph $\mathcal{G}$ is $\nu_{k}~{}:=~{}\sqrt{N}~{}\cdot~{}\max_{1\leqslant i\leqslant N}\left\{v_{k}(i)\right\}.$

It is shown in (Puy et al., 2015) that $\nu_{k}\in[k^{1/2},N^{1/2}]$.

Let $\mathsf{M}$ be a random sampling matrix constructed as in (6). For any $\delta,\epsilon\in\;]0,1[$,

for all $\bm{x}\in{\rm span}(\mathsf{U}_{k})$ with probability at least $1-\epsilon$ provided that

The above theorem presents a sufficient condition on $n$ ensuring that $\mathsf{M}$ satisfies the restricted isometry property (8). This condition is required to ensure that the solution of (7) is an accurate estimation of $\bm{c}_{j}$. The above theorem thus indicates that sampling $O(\nu_{k}^{2}\log{k})$ features is sufficient to recover the cluster indicator vectors.

For a simple regular graph $\mathcal{G}$ made of $k$ disconnected clusters, we have seen that $\mathsf{U}_{k}=(\bm{c}_{1},\ldots,\bm{c}_{k})$ up to a normalisation of the vectors. Therefore, $\nu_{k}=N^{1/2}/\min_{i}\{N_{i}^{1/2}\}$, where $N_{i}$ is the size of the $i^{\text{th}}$ cluster. If the clusters have the same size $N_{i}=N/k$ then $\nu_{k}=k^{1/2}$, the lower bound on $\nu_{k}$. In this simple optimal scenario, sampling $O(\nu_{k}^{2}\log{k})=O(k\log{k})$ features is thus sufficient to recover the cluster indicator vectors.

The attentive reader will have noticed that for graphs where $\nu_{k}^{2}\approx N$, no downsampling is possible. Yet, a simple solution exists in this situation: variable density sampling. Indeed, it is proved in (Puy et al., 2015) that, whatever the graph $\mathcal{G}$, there always exists an optimal sampling distribution such that $n=O(k\log{k})$ samples are sufficient to satisfy Eq. (8). This distribution depends on the profile of the local cumulative coherence and can be estimated rapidly (see (Puy et al., 2015) for more details). In this paper, we only consider uniform sampling to simplify the explanations, but keep in mind that in practice results will always be improved if one uses variable density sampling. Note also that one cannot expect to sample less than $k$ nodes to find $k$ clusters. Up to the extra $\log(k)$, our result is optimal.

CSC in practice

We have detailed the two fundamental theoretical notions supporting our algorithm, presented in Alg. 2. However, some steps in Alg. 2 still need to be clarified. In particular, Sec. 4.2 provides an extension of Theorem 3.2 that takes into account the use of a non-ideal low-pass filter (to handle the practical case where the order of the polynomial approximation is finite). This theorem in fine explains and justifies step 4 of Alg. 2. Then, in Sec. 4.3, important details are discussed such as the estimation of $\lambda_{k}$ (step 1) and the choice of the polynomial approximation (step 2). We finish this section with complexity considerations.

2 Non-ideal filtering of random signals

where the $\{\bm{\delta}_{i}^{r}\}$ are here Diracs in $n$ dimensions.

The normalisation of Step 4 in Alg. 2 approximates the action of $\mathsf{V}_{k}$ in the above equation. More details and justifications are provided in the “Important remark” at the end of this section. The distance between any two features reads

Approximation error. Denote $e(\lambda)$ the approximation error of the ideal low-pass filter:

In the form of graph filter operators, one has

We model the error $e$ using two parameters: $e_{1}$ (resp. $e_{2}$) the maximal error for $\lambda\leqslant\lambda_{k}$ (resp. $\lambda>\lambda_{k}$). We have

The resolution parameter. In some cases, the ideal reduced spectral distance $D_{ij}^{r}$ may be null. In such cases, approximating $D_{ij}^{r}=0$ using a non-ideal filter is not possible. In fact, non-ideal filtering introduces an irreducible error on the estimation of the feature vectors that is not possible to compensate in general. We thus introduce a resolution parameter $D_{min}^{r}$ below which the distances $D_{ij}^{r}$ do not need to be approximated exactly, but should remain below $D_{min}^{r}$ (up to a tolerated error).

Let $D_{min}^{r}\in\left]0,\sqrt{2}\right]$ be a chosen resolution parameter. For any $\delta\in\,]0,1]$, $\beta>0$, if $d$ is larger than

then, for all $(i,j)\in\{1,\ldots,n\}^{2}$,

with probability at least $1-2n^{-\beta}$ provided that

The proof is provided in the supplementary material.

Consequence of Theorem 4.1. All distances smaller (resp. larger) than the chosen resolution parameter $D_{min}^{r}$ are estimated smaller than $(1+\delta)D_{min}^{r}$ (resp. correctly estimated up to a relative error $\delta$). Moreover, for a fixed distance estimation error $\delta$, the lower we decide to fix $D_{min}^{r}$, the lower should also be the errors $e_{1}$ and/or $e_{2}$ to ensure that Eq. (9) still holds, which implies an increase of the order $p$ of the polynomial approximation of the ideal filter $h_{\lambda_{k}}$, and ultimately, that means a higher computation time for the filtering operation of the random signals.

The polynomial approximation. Theorem 4.1 uses a separate control on $e(\lambda)$ below $\lambda_{k}$ (with $e_{1}$) and above $\lambda_{k}$ (with $e_{2}$). To have such a control in practice, one would need to use rational filters (ratio of two polynomials) to approximate $h_{\lambda_{k}}$. Such filters have been introduced in the graph context (Shi et al., 2015), but they involve another optimisation step that would burden our main message. We prefer to simplify our analysis by using polynomials for which only the maximal error can be controlled. We write

In this easier case, one can show that Theorem 4.1 is still valid if Eq. (9) is replaced by

In our experiments, we could follow (Shuman et al., 2011) and use truncated Chebychev polynomials to approximate the ideal filter, as these polynomials are known to require a small degree to ensure a given tolerated maximal error $e_{m}$. We prefer to follow (Napoli et al., 2013) who suggest to use Jackson-Chebychev polynomials: Chebychev polynomials to which are added damping multipliers to alleviate the unwanted Gibbs oscillations around the cut-off frequency $\lambda_{k}$.

The polynomial’s order $p$. For a fixed $\delta$, $D_{min}^{r}$, and $\min_{i}\{v_{k}(i)\}$, one should use the Jackson-Chebychev polynomial of smallest order $p^{*}$ ensuring that $e_{m}$ satisfies Eq. (11), in order to optimize the computation time while making sure that Theorem 4.1 applies. Studying $p^{*}$ theoretically without computing the Laplacian’s complete spectrum (see Eq. (10)) is beyond the scope of this paper. Experimentally, $p=50$ yields good results (see Fig. 2c).

Estimation of $\lambda_{k}$. The fast filtering step is based on the polynomial approximation of $h_{\lambda_{k}}$, which is itself parametrized by $\lambda_{k}$. Unless we compute the first $k$ eigenvectors of $\mathsf{L}$, thereby partly loosing our efficiency edge on other methods, we cannot know the value of $\lambda_{k}$ with infinite precision. To estimate it efficiently, we use eigencount techniques (Napoli et al., 2013): based on low-pass filtering with a cut-off frequency at $\lambda$ of random signals, one obtains an estimation of the number of enclosed eigenvalues in the interval $[0,\lambda]$. Starting with $\lambda=2$ and proceeding by dichotomy on $\lambda$, one stops the algorithm as soon as the number of enclosed eigenvalues equals $k$. For each value of $\lambda$, in order to have a proper estimation of the number of enclosed eigenvalues, we choose to filter $2\log{N}$ random signals with Jackson-Chebychev polynomial approximation of the ideal low-pass filters.

4 Complexity considerations

The complexity of steps 2, 3 and 5 of Alg. 2 are not detailed as they are insignificant compared to the others. First, note that fast filtering a graph signal costs $O(p\,\#\mathcal{E})$.Recall that $p$ is the order of the polynomial filter. Therefore, Step 1 costs $O(p\,\#\mathcal{E}\log{N})$ per iteration of the dichotomy, and Step 4 costs $O(p\,\#\mathcal{E}\log{n})$ (as $d=O(\log{n})$). Step 7 requires to solve Eq. (7) with the polynomial approximation of $g_{\lambda_{k}}(\lambda)=1-h_{\lambda_{k}}(\lambda)$. When solved, e.g., by conjugate gradient or gradient descent, this step costs a fast filtering operation per iteration of the solver and for each of the $k$ classes. Step 7 thus costs $O(p\,\#\mathcal{E}k)$. Also, the complexity of $k$-means to cluster $Q$ feature vectors of dimension $r$ into $k$ classes is $O(kQr)$ per iteration. Therefore, Step 6 with $Q=n$ and $r=d=O(\log(n))$ costs $O(kn\log{n})$. CSC’s complexity is thus $O\left(kn\log{n}+p\,\#\mathcal{E}\left(\log{N}+\log{n}+k\right)\right).$ In practice, we are interested in sparse graphs: $\,\#\mathcal{E}=O(N)$. Using the fact that $n=O(k\log{k})$, CSC’s complexity simplifies to

SC’s $k$-means step has a complexity of $O(Nk^{2})$ per iteration. In many casesRoughly, all cases for which $k^{2}>p(\log{N}+k)$. this sole task is more expensive than the CSC algorithm. On top of this, SC has the additional complexity of computing the first $k$ eigenvectors of $\mathsf{L}$, for which the cost of ARPACK - a popular eigenvalue solver - is $O(k^{3}+Nk^{2})$ (see, e.g., Sec. 3.2 of (Chen et al., 2011)).

This study suggests that CSC is faster than SC for large $N$ and/or $k$. The above algorithms’ number of iterations are not taken into account as they are difficult to predict theoretically. Yet, the following experiments confirm the superiority of CSC over SC in terms of computational time.

Experiments

We first perform well-controlled experiments on the Stochastic Block Model (SBM), a model of random graphs with community structure, that was showed suitable as a benchmark for SC in (Lei & Rinaldo, 2015). We also show performance results on a large real-world network. Implementation was done in Matlab R2015a, using the built-in function kmeans with 20 replicates, and the function eigs for SC. Experiments were done on a laptop with a 2.60 GHz Intel i7 dual-core processor running OS Fedora release 22 with 16 GB of RAM. The fast filtering part of CSC uses the gsp$\_$cheby$\_$op function of the GSP toolbox (Perraudin et al., 2014). Equation (7) is solved using Matlab’s gmres function. All our results are reproducible with the CSCbox downloadable at http://cscbox.gforge.inria.fr/.

What distinguishes the SBM from Erdos-Renyi graphs is that the probability of connection between two nodes $i$ and $j$ is not uniform, but depends on the community label of $i$ and $j$. More precisely, the probability of connection between nodes $i$ and $j$ equals $q_{1}$ if they are in the same community, and $q_{2}$ if not. In a first approach, we look at graphs with $k$ communities, all of same size $N/k$. Furthermore, instead of considering the probabilities, one may fully characterize a SBM by providing their ratio $\epsilon=\frac{q_{2}}{q_{1}}$, as well as the average degree $s$ of the graph. The larger $\epsilon$, the more difficult the community structure’s detection. In fact, Decelle et al. (2011) show that a critical value $\epsilon_{c}$ exists above which community detection is impossible at the large $N$ limit: $\epsilon_{c}=(s-\sqrt{s})/(s+\sqrt{s}(k-1))$.

2 Performance results

In Figs. 2 a-d), we compare the recovery performance of CSC versus SC for different parameters. The performance is measured by the Adjusted Rand similarity index (Hubert & Arabie, 1985) between the SBM’s ground truth and the obtained partitions. It varies between $-1$ and $1$. The higher it is, the better is the reconstruction. These figures show that the performance of CSC saturates at the default values of $n,d,p$ and $\gamma$ (see top of Alg. 2). Experiments on the SBM with heterogeneous community sizes are provided in the supplementary material and show similar results.

Fig. 2 e) shows the estimation results of $\lambda_{k}$ for different values of $\epsilon$ : it is overestimated in the SBM context. As long as the estimated value stays under $\lambda_{k+1}$, this overestimation does not have a strong impact on the method. On the other hand, as $\epsilon$ becomes larger than $\sim 0.06$, our estimation of $\lambda_{k}$ is larger than $\lambda_{k+1}$, which means that our feature vectors start to integrate some unwanted information from eigenvectors outside of ${\rm span}{(\mathsf{U}_{k})}$. Even though the impact of this additional information is application-dependent and in some cases insignificant, further efforts to improve the estimation of $\lambda_{k}$ would be beneficial to our method.

In Figs. 2 f-g) we fix $\epsilon$ to $\epsilon_{c}/4$, $n,d,p$ and $\gamma$ to the values given in Alg. 2, and vary $N$ and $k$. We compare the recovery performance and the time of computation of CSC, SC and Boutsidis’ power method (Boutsidis & Gittens, 2015). The power method (PM), in a nutshell, 1) applies the Laplacian matrix to the power $r$ to $k$ random signals, 2) computes the left singular vectors of the $N\times k$ obtained matrix, to extract feature vectors, 3) applies $k$-means in high-dimension (like SC) with these feature vectors. In our experiments, we use $r=10$. The recovery performances are nearly identical in all situations, even though CSC is only a few percents under SC and PM (Fig. f is zoomed around the high values of the recovery score). For the time of computation, the experiments confirm that all three methods are roughly linear in $N$ and polynomial in $k$ (Fig. g is plotted in log-log), with a lower exponent for CSC than for SC and PM; such that SC and PM are faster for $k=20$ but CSC becomes up to an order of magnitude faster as $k$ increases to 200. Note that the SBM is favorable to SC as Matlab’s function eigs converges very fast in this case, e.g., for $N=10^{5}$, it finds the first $k=200$ eigenvectors in less than 2 minutes! PM sidesteps successfully the cost of eigs, but the cost of $k$-means in high-dimension is still a strong bottleneck.

We finally compare CSC and SC on a real-world dataset: the Amazon co-purchasing network (Yang & Leskovec, 2015). It is an undirected connected graph comprising $N=334\,863$ nodes and $\#\mathcal{E}=925\,872$ edges. The results are presented in Fig.2 h) for three values of $k$. As there is no clear ground truth in this case, we use the modularity (Newman & Girvan, 2004) to measure the algorithm’s clustering performance, a well-known cost function that measures how well a given partition separates a network in different communities. Note that the 20 replicates of $k$-means would not converge for SC with the default maximum number of iterations set to $100$. For a fair comparison with CSC, we used only 2 replicates with a maximum number of iterations set to $1000$ for SC’s $k$-means step. We see that for the same clustering performance, CSC is much faster than SC, especially as $k$ increases. The PM algorithm on this dataset does not perform well: even though the features are estimated quickly, they apparently do not form clear classes such that its $k$-means step takes even longer than SC’s. For the three values of $k$, we stopped the PM algorithm after a time of computation exceeding SC’s.

Conclusion

By graph filtering $O(\log k)$ random signals, we construct feature vectors whose interdistances approach the standard SC feature distances. Then, building upon compressive sensing results, we show that one can sample $O(k\log k)$ nodes from the set of $N$ nodes, cluster this reduced set of nodes and interpolate the result back to the whole graph. If the low-dimensional $k$-means result is correct, i.e., if Eq. (3) is verified, we guarantee that the interpolation is a good approximation of the SC result. To improve the clustering result of the reduced set of nodes, one could consider the concept of community cores (Seifi et al., 2013). In fact, as the filtering and the low-dimensional clustering steps are fairly cheap to compute, one could repeat these steps for different random signals, keep the sets of nodes that are always classified together and use only these stable “cores” for interpolation. Our experiments show that even without such potential improvements, CSC proves efficient and accurate in synthetic and real-world datasets; and could be preferred to SC for large $N$ and/or $k$.

Acknowledgments

This article was submitted when G. Puy was with INRIA Rennes - Bretagne Atlantique, France. This work was partly funded by the European Research Council, PLEASE project (ERC-StG-2011-277906), and by the Swiss National Science Foundation, grant 200021-154350/1 - Towards Signal Processing on Graphs.

Appendix A Proof of Theorem 3.2

where the $\bm{f}_{i}$ are the standard SC feature vectors. Applying Theorem 1.1 of (Achlioptas, 2003) (an instance of the Johnson-Lindenstrauss lemma) to $\left\|\mathsf{R}^{\intercal}\mathsf{U}_{k}(\bm{f}_{i}-\bm{f}_{j})\right\|$, the following holds. If $d$ is larger than:

then with probability at least $1-N^{-\beta}$, we have, $\forall(i,j)\in\{1,\ldots,N\}^{2}$:

As the columns of $\mathsf{U}_{k}$ are orthonormal, we end the proof:

Appendix B Proof of Theorem 4.1

We continue the proof by bounding $\left\|\mathsf{R}^{\intercal}\mathsf{H}_{\lambda_{k}}^{\intercal}\mathsf{V}_{k}^{\intercal}\mathsf{M}^{\intercal}\,\bm{\delta}_{ij}^{r}\right\|$ and $\left\|\mathsf{R}^{\intercal}\mathsf{E}^{\intercal}\mathsf{V}_{k}^{\intercal}\mathsf{M}^{\intercal}\,\bm{\delta}_{ij}^{r}\right\|$ separately.

Let $\delta\in]0,1]$. To bound $\left\|\mathsf{R}^{\intercal}\mathsf{H}_{\lambda_{k}}^{\intercal}\mathsf{V}_{k}^{\intercal}\mathsf{M}^{\intercal}\bm{\delta}_{ij}^{r}\right\|$, we set $\epsilon=\delta/2$ in Theorem 3.2. This proves that if $d$ is larger than

then with probability at least $1-n^{-\beta}$,

for all $(i,j)\in\{1,\ldots,n\}^{2}$. To bound $\left\|\mathsf{R}^{\intercal}\mathsf{E}^{\intercal}\mathsf{V}_{k}^{\intercal}\mathsf{M}^{\intercal}\bm{\delta}_{ij}^{r}\right\|$, we use Theorem $1.1$ in (Achlioptas, 2003). This theorem proves that if $d>d_{0}$, then with probability at least $1-n^{-\beta}$,

for all $(i,j)\in\{1,\ldots,n\}^{2}$. Using the union bound and (B), we deduce that, with probability at least $1-2n^{-\beta}$,

for all $(i,j)\in\{1,\ldots,n\}^{2}$ provided that $d>d_{0}$.

Then, as $e$ is bounded by $e_{1}$ on the first $k$ eigenvalues of the spectrum and by $e_{2}$ on the remaining ones, we have

Define, for all $(i,j)\in\{1,\ldots,n\}^{2}$:

Thus, the above inequality may be rewritten as:

for all $(i,j)\in\{1,\ldots,n\}^{2}$, which combined with (B) yields

for all $(i,j)\in\{1,\ldots,n\}^{2}$, with probability at least $1-2n^{-\beta}$ provided that $d>d_{0}$.

Let us now separate two cases. In the case where $D_{ij}^{r}\geqslant D_{min}^{r}>0$, we have

provided that Eq. (7) of the main paper holds. Combining the last inequality with (B) proves the first part of the theorem.

In the case where $D_{ij}^{r}<D_{min}^{r}$, we have

provided that Eq. (7) of the main paper holds. Combining the last inequality with (B) terminates the proof. ∎

Appendix C Experiments on the SBM with heterogeneous community sizes

We perform experiments on a SBM with $N=10^{3},k=20,s=16$ and hetereogeneous community sizes. More specifically, the list of community sizes is chosen to be: $5$, $10$, $15$, $20$, $25$, $30$, $35$, $40$, $45$, $50$, $50$, $55$, $60$, $65$, $70$, $75$, $80$, $85$, $90$ and $95$ nodes. In this scenario, there is no theoretical value of $\epsilon$ over which it is proven that recovery is impossible in the large $N$ limit. Instead, we vary $\epsilon$ between and $0.2$ and show the recovery performance results with respect to $n$, $d$, $p$ and $\gamma$ in Fig. 2. Results are similar to the homogeneous case presented in Fig. 1(a-d) of the main paper.

References