Lorenzo Dall'amico
email: lorenzo.dall-amico@gipsa-lab.fr

Romain Couillet
email: romain.couillet@univ-grenoble-alpes.fr

Nicolas Tremblay
email: nicolas.tremblay@gipsa-lab.fr

Tina Eliassi-Rad

A Unified Framework for Spectral Clustering in Sparse Graphs

Keywords: community detection, sparsity, heterogeneous degree distribution, spectral clustering, unsupervised learning ©2021 Lorenzo Dall'Amico, Romain Couillet

This article considers spectral community detection in the regime of sparse networks with heterogeneous degree distributions, for which we devise an algorithm to efficiently retrieve communities. Specifically, we demonstrate that a well parametrized form of regularized Laplacian matrices can be used to perform spectral clustering in sparse networks without suffering from its degree heterogeneity. Besides, we exhibit important connections between this proposed matrix and the now popular non-backtracking matrix, the Bethe-Hessian matrix, as well as the standard Laplacian matrix. Interestingly, as opposed to competitive methods, our proposed improved parametrization inherently accounts for the hardness of the classification problem. These findings are summarized under the form of an algorithm capable of both estimating the number of communities and achieving high-quality community reconstruction.

Introduction

Graph theory has found many applications in a variety of domains that span from modern biology, to technology and social sciences [START_REF] Barabási | Network science[END_REF]. Although the underlying represented systems may be fundamentally different, some common features emerge in complex networks. Some of these are [START_REF] Mark Ej Newman | The structure and function of complex networks[END_REF]:

1. Their heterogeneous degree distribution: the number of connections of each node (their degree) is far from following a homogeneous distribution and, instead, is typically a broad distribution [START_REF] Barabási | Emergence of scaling in random networks[END_REF] in which few nodes (called hubs) have a very large degree while the majority has only few connections. 2. Sparsity: the average degree is typically much smaller than the size of the graph and a small fraction of all the possible connections are present [START_REF] Barabási | Network science[END_REF]; consequently, the average degree can be considered as roughly independent of the size of the graph.

3. the clustering effect: as a consequence of the nodes affinity encoded by edges [START_REF] Stephen | On network theory[END_REF], in real graphs typically emerge small groups of densely connected nodes, called clusters or communities.

The remainder of the article focuses on community detection, that is in determining a good community assignment of the nodes of a given graph. While point 3. of the list above is of central concern, the (nuisance) properties 1. and 2. need to be accounted for to guarantee high performance community assignments on real-world graphs.

Community Detection

One of the most natural tasks in graph theory is community detection, i.e., the identification of similarity groups on a given graph. Practically, for an unweighted and undirected graph G(V, E) with |V| = n nodes and |E| edges, community detection consists in finding a nonoverlapping partition of the nodes that identifies underlying communities in a completely unsupervised manner. There is no unique definition of a community, but a general criterion is to impose that nodes in the same community have more inter-connections than nodes in different communities, as a consequence of the stronger affinity among members of the same community [START_REF] Fortunato | Community detection in graphs[END_REF].

A possible way of formalizing this intuition consists in defining the class labels as the solution of an optimization problem such as MinCut, RatioCut, and NormalizedCut [START_REF] Von | A tutorial on spectral clustering[END_REF]. These definitions do not make any assumption on the underlying graph but only on what a satisfactory community assignment should be like. The aforementioned optimization problems are, however, NP-hard and only approximate solutions can be obtained. Another possible approach (adopted in the following) is instead to consider community detection as a statistical inference problem. The graph G(V, E) is seen as a realization of a random process in which the class-label assignment is encoded by some hidden parameters of the generative model that have to be inferred. To simultaneously account for sparsity and heterogeneity that typically characterize real graphs, in this article we consider graphs generated according to a k-class sparse degree-corrected stochastic block model (DC-SBM) [START_REF] Karrer | Stochastic blockmodels and community structure in networks[END_REF], that is

P(A ij = 1| i , j , θ i , θ j ) = θ i θ j C i , j n 1 ≤ i < j ≤ n (1)
in which = {1, • • • , k} n denotes the node-wise labelling vector ( i = p if node i is in class p), θ is the vector of intrinsic connection "probabilities" which are used to produce an arbitrary degree distribution and are independent of the label vector . The entries of θ are positive independent random variables θ i ∈ [θ min , θ max ], with θ min > 0, having first moment E[θ i ] = 1 and second moment E[θ 2 i ] = Φ. We denote with C ∈ R k×k the symmetric class-affinity matrix, with entries independent of n. The term 1/n in (1) bounds the average degree d = 1 n 1 T n A1 n = O n (1) to an n-independent value, setting the problem in the sparse regime. A major advantage of a well defined generative model results in statistically tractable conclusions, like the existence in asymptotically large graphs (n → ∞) generated from the (degree-corrected) stochastic block model of a limiting detectability threshold. Specifically, for k = 2 classes, one can in general identify a parameter α, function of the community statistics and such that, beyond a threshold (α > α c ), partial label reconstruction can be theoretically achieved, whereas below the threshold (α < α c ) no algorithm can perform better than random guess. For k > 2 classes, the situation is more involved: it is conjectured [START_REF] Decelle | Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications[END_REF] that there exists an easy detection zone in which a non-trivial clustering can be found in polynomial time, a hard detection zone in which non-trivial clustering can be found but only in exponential time and an impossible detection zone in which no algorithm can find a non-trivial partition.

In the very specific setting where k = 2 and both classes are of equal size (n/2), the detectability threshold assumes a simple expression. Letting C i , j = c in if i = j and c out otherwise, in the stochastic block model (for which θ = 1 n ) it was initially conjectured [START_REF] Decelle | Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications[END_REF] and later proved [START_REF] Mossel | Stochastic block models and reconstruction[END_REF][START_REF] Massoulié | Community detection thresholds and the weak ramanujan property[END_REF] that the condition to non-trivial clustering is given by

α ≡ c in -c out √ c > 2 ≡ α c , (2)
where c = (c in + c out )/2. Equation ( 2) was later generalized in [START_REF] Gulikers | An impossibility result for reconstruction in the degree-corrected stochastic block model[END_REF](Gulikers et al., , 2017a) ) to the DC-SBM but now with α c = 2/ √ Φ (we recall that Φ = E[θ 2 i ]). If on one side random generative models like the DC-SBM give a solid theoretical understanding, the applicability to real graphs of a DC-SBM inspired clustering algorithm is a serious concern to be taken into account and which will be addressed here.

Spectral Clustering: Related Work

Regardless of the definition of a community considered (as a solution of an optimization problem or of an inference problem), different classes of algorithms can be considered to find exact or approximate solutions for community detection. Spectral clustering is a class of algorithms according to which each node of the graph is mapped to a low dimensional space. Spectral clustering hence provides a small dimensional graph embedding from which nodes are partitioned using common clustering algorithms, like k-means. The name spectral comes from the fact that the embedding is obtained using some eigenvectors (later referred to as informative) of a suitable graph representation matrix. Spectral clustering allows one to obtain an approximate fast solution of the optimal class assignment both in the case in which clusters are defined as the solution of an NP-hard optimization problem [START_REF] Von | A tutorial on spectral clustering[END_REF][START_REF] Zhang | Understanding regularized spectral clustering via graph conductance[END_REF] or of an inference problem [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF].

Algorithm 0 A typical spectral clustering algorithm on a graph G with k classes 1: Input : G, k 2: Choose the matrix M ∈ R n×n , a suited representation of the graph 3: Find the k largest (or smallest) eigenvalues of M and stack the corresponding eigenvectors in the columns of a matrix X ∈ R n×k 4: (Optional ) Normalize the rows of X 5: Estimate community labels as the output of some small dimensional clustering method (such as k-means or expectation maximization) performed on the points in R k defined by the rows of X 6: return Estimated label community vector. Algorithm 0 provides the general structure of a k-class spectral clustering algorithm based on an "affinity" matrix M . Algorithm 0 relies on the fact that some of the eigenvectors of these appropriate affinity matrices contain the class structure of the graph. Step 3 of Algorithm 0 underlines that the informative eigenvalues are generally found to be the largest or smallest eigenvalues, so that there is an algorithmically clear definition of which eigenvectors should be used for clustering.

In order to review the most relevant spectral clustering algorithms for community detection, let us define the adjacency matrix A ∈ {0, 1} n×n of the graph G by A ij = 1 (ij)∈E where 1 x equals 1 if the condition x is verified and zero otherwise, and the diagonal degree matrix D = diag(A1 n ), where 1 n is the vector with all entries equal to one. In [START_REF] Fiedler | Algebraic connectivity of graphs[END_REF], and in the case of k = 2 communities, it was proposed to reconstruct communities using the eigenvector corresponding to the second smallest eigenvalue of the combinatorial graph Laplacian matrix L = D -A. It was then shown (see e.g. [START_REF] Von | A tutorial on spectral clustering[END_REF] that this eigenvector provides a relaxed solution of the RatioCut problem. Based on this result, one can build a spectral algorithm that, referring to Algorithm 0, sets M = L; in this case, the informative eigenvectors correspond to the smallest eigenvalues and step 4 is not performed. Similarly, the normalized graph Laplacian matrices (L rw = D -1 A and L sym = D -1/2 AD -1/2 ) can be used to solve a relaxed form of the NCut problem [START_REF] Von | A tutorial on spectral clustering[END_REF]. In [START_REF] Shi | Normalized cuts and image segmentation[END_REF], M = L rw , the informative eigenvectors correspond to the k largest eigenvalues and step 4 is again not performed, while in [START_REF] Andrew | On spectral clustering: Analysis and an algorithm[END_REF], M = L sym , the informative eigenvectors correspond to the k largest eigenvalues and step 4 is performed. Another classical choice of matrix M is in [START_REF] Newman | Modularity and community structure in networks[END_REF], where the authors propose an alternative method to define communities through a NP-hard optimization problem by maximizing the so-called modularity of the graph. For k = 2 classes, a relaxed solution can be obtained by considering the largest eigenvector of the modularity matrix A -dd T 2|E| , where d = A1 n is the degree vector. For more than two communities, the algorithm is slightly more involved and deviates from the general structure of Algorithm 0.

Many results justifying the earlier methods exist in the asymptotic limit where n 1, in the so-called dense regime, in which the average degree (number of connections) grows with the size (n) of the network. However, as we anticipated in the introduction, real networks tend instead to be sparse as well as heterogeneous in their degree distribution, following broad distributions such as power laws. These two characteristics (sparsity and heterogeneity) make it difficult to theoretically predict the behavior and performances of Algorithm 0 and to choose a proper matrix M to work with. From a linear algebra viewpoint, sparsity in the node degrees has been shown to spread the eigenvalues of the Laplacian matrix, with the deleterious effect to "swallow" the isolated (smallest or largest) informative eigenvalues within the so-called "bulk" of uninformative eigenvalues: as a result, the corresponding informative eigenvectors are no longer found to be associated with dominant (largest or smallest) eigenvalues, as shown in Figure 1B, and as opposed to Figure 1A. Besides, by losing the isolation of informative eigenvalues, the associated eigenvectors tend to merge with the eigenvectors associated to close-by (non-informative) eigenvalues. Heterogeneity in the degrees also induces a severe negative effect: it modifies the amplitude of the i-th entry (for all 1 ≤ i ≤ n) of the informative eigenvectors by a non-trivial function of the degree of node i. This further compromises the efficiency of the last classification step (usually performed through k-means), as shown in Figure 1D, to be opposed to Figure 1C. Several contributions have thus tackled the challenging problem of devising efficient spectral clustering in the sparse regime and for an arbitrary degree distribution.

One prominent line of research was devoted to defining spectral algorithms capable of retrieving communities as soon as theoretically possible (α > α c , Equation (2)) on sparse graphs generated from the SBM, for which we recall that there is a sharp transition between the "detectable" and "undetectable" situations. The authors of [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF] proposed an algorithm based on the (non-symmetric) non-backtracking matrix B (ij)(lm) = δ jl (1 -δ im ), with B ∈ {0, 1} 2|E|×2|E| . In [START_REF] Massoulié | Community detection thresholds and the weak ramanujan property[END_REF] it was indeed shown that the eigenvectors associated to the largest (in modulus) eigenvalues of B have a non-trivial correlation with the actual underlying communities, as soon as α > α c . A closely related algorithm is the one proposed in [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF] which instead uses the eigenvectors attached to the smallest eigenvalues of the Bethe-Hessian matrix H r = (r 2 -1)I n + D -rA ∈ R n×n , for r = ρ(B) (B being the non-backtracking matrix just defined and ρ(•) the spectral radius). These algorithms have the benefit of achieving non-trivial classification down to the detectability threshold, but they have only been investigated under a sparse SBM setting; hence they do not cope with degree heterogeneity. The extension to the sparse DC-SBM case was treated in (Gulikers et al., 2017a), in which the authors show that spectral clustering on B also works down to the (generalized) threshold.

All these results are powerful as they propose algorithms capable of reaching the information-theoretic threshold, but they also have inherent weaknesses as they only guarantee a positive correlation of their output classification with the underlying true structure. Specifically, these positive correlations do not imply that the classification performance is maximal. In particular, even in (Gulikers et al., 2017a) where spectral clustering on the DC-SBM is studied, the problem of eigenvector pollution due to degree heterogeneity (which is one of our main focuses in this article) is not discussed and a fortiori not corrected.

Another line of research which studies the reconstruction of sparse DC-SBM networks suggests to exploit regularization (i.e., using D τ = D + τ I n and A τ = A + τ 1 n 1 T n in place of D and A) as a solution to maintain the location of the informative eigenvalues in their dominant (smallest or largest) positions [START_REF] Arash | Pseudo-likelihood methods for community detection in large sparse networks[END_REF][START_REF] Joseph | Impact of regularization on spectral clustering[END_REF][START_REF] Lei | Consistency of spectral clustering in stochastic block models[END_REF][START_REF] Can M Le | Concentration of random graphs and application to community detection[END_REF]. A particularly interesting method in this direction is proposed in [START_REF] Qin | Regularized spectral clustering under the degree-corrected stochastic blockmodel[END_REF] which uses the regularized symmetric reduced Laplacian

matrix L sym τ = D -1/2 τ AD -1/2 τ
, for τ = d, the average degree. In terms of Algorithm 0, the algorithm proposed in [START_REF] Qin | Regularized spectral clustering under the degree-corrected stochastic blockmodel[END_REF] sets M = L sym τ , searches for the k eigenvectors corresponding to the k largest eigenvalues of M , and then performs the normalization step 4 on the rows of the resulting matrix X. While [START_REF] Qin | Regularized spectral clustering under the degree-corrected stochastic blockmodel[END_REF]) deals with graphs with arbitrary degree distribution (generated from the DC-SBM), the authors do not discuss whether their algorithm can achieve non-trivial clustering down to the threshold, i.e. as soon as α > α c . We will bring new conclusions in this direction in Section 3.

Contributions

In the state-of-the-art methods we presented, while sparsity is not properly accounted for in [START_REF] Shi | Normalized cuts and image segmentation[END_REF][START_REF] Andrew | On spectral clustering: Analysis and an algorithm[END_REF], the line of research based on B (and consequently H r ) [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF][START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF][START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF]Gulikers et al., 2017a) does account for sparsity and provides methods to reach the detectability threshold; yet, all they guarantee is the possibility to obtain a positive, possibly suboptimal, correlation between the algorithm output and the underlying community structure. As for the works on regularization [START_REF] Qin | Regularized spectral clustering under the degree-corrected stochastic blockmodel[END_REF][START_REF] Can | Estimating the number of communities in networks by spectral methods[END_REF], they only establish theoretical results of perfect community recovery far from the (most interesting) detection threshold.

The present work solves both the issues of sparsity and heterogeneity at once by proposing a simple spectral clustering scheme which, under a sparse DC-SBM setting, provably performs non-trivial clustering as soon as α > α c and is robust to degree heterogeneity, as it retrieves eigenvectors not infected by the node degrees. This claim is supported by three parallel arguments having the side advantage to unify in a joint framework the ideas of [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF][START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF][START_REF] Shi | Normalized cuts and image segmentation[END_REF][START_REF] Fiedler | Algebraic connectivity of graphs[END_REF][START_REF] Qin | Regularized spectral clustering under the degree-corrected stochastic blockmodel[END_REF].

We further discuss how our algorithm, analyzed under the DC-SBM assumption, can be extended to real graphs. In particular, the algorithm provides an accurate estimate of the number of communities when unknown, and is observed to systematically achieve good scores (both in terms of modularity or log-likelihood of the DC-SBM posterior probability) on all real-world graphs which we experimented on.

In summary, our main contributions are as follows:

1. We devise a practical spectral algorithm exploiting the eigenvectors corresponding to the smallest eigenvalues of the Bethe-Hessian matrix H r for a set of properly chosen values of r. We further provide strong arguments and extensive numerical simulations to claim that the proposed method achieves non-trivial reconstruction down to the DC-SBM detectability threshold. 2. The proposed algorithm is extensively tested on synthetic DC-SBM graphs (for which the performance is measured in terms of overlap between the estimated and groundtruth labels) as well as on real networks (for which the performance is measured in terms of modularity and DC-SBM log-likelihood). In both cases our algorithm outperforms or is on par with the standard competing spectral techniques for sparse graphs. 3. In passing, we provide a new vision to spectral clustering, and in particular a compelling new connection between our proposed approach and all aforementioned standard spectral methods, so far treated independently. In particular, we show that the regularized reduced random walk Laplacian matrix L rw τ = D -1 τ A (that shares the same eigenvalues as L sym τ ) can be used to perform non-trivial clustering down to the detectability threshold for a choice of τ strongly related to our proposed algorithm.

A Python implementation of the codes needed to reproduce the results of this article can be found at lorenzodallamico.github.io/codes/: all data, experiments and running times refer to this implementation. On top of the Python codes, we developed a more efficient Julia implementation of the proposed Algorithm 2, available in the CoDeBetHe.jl package which gains a factor x5 ∼ x10 in terms of computational time over the Python implementation.

The remainder of the article is organized as follows: Section 2 formally introduces the problem at hand and provides three complementary supporting arguments to the claim that the eigenvectors corresponding to the smallest eigenvalues of a set of Bethe-Hessian matrices H r , with different values of r, can be used to reconstruct classes. Section 3 draws the connection between the above H r and the regularized Laplacian and proves that, under a correct choice of τ related to r above, the matrix L rw τ is also a suitable candidate to reconstruct communities. Section 4 introduces the final form of the proposed algorithm and provides simulation outputs on real networks. Section 5 finally closes the article.

Notation

• Matrices are indicated with standard font capital letters (M ). The only exception is X, the matrix concatenating the eigenvectors of the spectral clustering algorithm. • M •,p indicates the p-th column of M and M p,• its p-th row.

• Vectors are denoted in bold v p , while scalar and vector elements in standard font v p,i .

• We denote by s ↑ i (M ) the i-th smallest eigenvalue of a Hermitian matrix M , i.e.,

s ↑ 1 (M ) ≤ s ↑ 2 (M ) ≤ • • • ≤ s ↑ n (M )
, and s ↓ i (M ) the i-th largest. For a non-Hermitian matrix M , we denote by s

↑/↓ R,i (M ) the i-th smallest/largest real eigenvalue, while s ↑/↓ |•| i (M ) indicates the i-th smallest/largest in modulus. When using the notation s i (M ) we consider a generic eigenvalue of M . • S(M ) = {s i (M ), 1 ≤ i ≤ n} denotes the set of eigenvalues of M ∈ R n×n .
• The notation 1 n indicates the vector of size n containing all ones:

1 T n = (1, 1, • • • , 1).
• We denote the neighborhood of the node i of a graph with adjacency matrix A by ∂i = {j ∈ V : A ij = 1}.

Informative Eigenvectors of B and H r in the Sparse Regime

Consider an undirected and unweighted graph G(V, E) with |V| = n nodes and |E| edges. Our objective is to devise a community detection algorithm on G which is resilient to the typically heterogeneous and sparse nature of real graphs, and also accounts for the fact that the actual grouping of nodes in communities is not clear-cut in a real graph (i.e., there usually exists no "optimal" number of communities and no "ground-truth" for the allocation of nodes into communities). While our proposed algorithm (Algorithm 2) is applicable to any arbitrary graph (as discussed in Section 4), for analytic purposes, we start the article by assuming that G is constituted of exactly k communities and generated as a sparse DC-SBM graph; that is, the edges of G are drawn independently according to Equation (1). Before delving into the technical details though, we first provide a short digest of our main result.

Informal Statement of the Main Result

Let us present our central result (Claim 1) and the resulting new algorithm for community detection in simple terms. The result revolves around the shape of the eigenvectors attached to the smallest (algebraic) eigenvalues of the Bethe-Hessian matrix H r ∈ R n×n , defined as

H r = (r 2 -1)I n + D -rA, (3)
for carefully chosen values of r ≥ 1; we recall that A and D are the adjacency and degree matrices of G(V, E), respectively, while I n is the identity matrix of size n.

To each k-class graph, community detection is performed by first associating a set of k Bethe-Hessian matrices {H ζp } p=1,...,k , with ζ p defined so that the p-th smallest eigenvalue of H ζp is equal to zero. We then extract from each matrix the eigenvector x p attached to the zero eigenvalue (so that H ζp x p = 0). The eigenvectors {x p } p=1,...,k are stacked in the columns of the matrix X defined in Algorithm 0 and used to produce the small dimensional node embedding on which the k-means algorithm is applied.

Claim 1 justifies the relevance of the above procedure as a highly performing community detection method. Specifically, for a k-class DC-SBM graph, Claim 1 states that:

• The largest value of ζ p that can be defined is ζ k . We show that, as a consequence, the (k + 1)-th smallest eigenvalue of H r is always positive for all r and the maximal number of negative eigenvalues of H r as function of r is precisely equal to k. This allows one to build an estimator for the number of classes. • The eigenvectors {x p } p=1,...,k and the k dominant eigenvectors of the non-backtracking matrix B are correlated to the class structure under the same hypothesis. This, for k = 2 classes of equal size, means that {x 1 , x 2 } can be used to reconstruct communities as soon as theoretically possible, i.e. when α > α c , (Equation ( 2)). • The entries of {x p } p=1,...,k are not polluted by the degree heterogeneity; hence do not suffer the problem shown in Figure 1D that hinders the performance of k-means.

The properties of {x p } allow us to define an efficient algorithm for community detection in sparse graphs with a heterogeneous degree distribution, solving at once the two main challenges discussed in Section 1.

Main Result

In this section, we proceed with the formal statement of Claim 1, together with a clear definition of the hypotheses under which it is formulated.

Model and Setting

Let G be a k class DC-SBM generated graph. We consider graphs G whose classes are not necessarily of equal size and denote by 0 < π p < 1, for every class 1 ≤ p ≤ k, the fraction of nodes belonging to class p. We correspondingly define the diagonal matrix Π = diag(π 1 , . . . , π k ) (note that in particular Tr(Π) = 1). We next denote by (ν p , v p ) the eigenvalue-eigenvector pairs of CΠ, i.e., CΠv p = s ↓ p (CΠ)v p ≡ ν p v p , where we recall s ↓ p (CΠ) denotes the p-th largest eigenvalue of CΠ, hence the eigenvalues are sorted as ν 1 ≥ . . . ≥ ν k . Note that the eigenvalues of CΠ are all real as CΠ has the same spectrum as Π 1/2 CΠ 1/2 , which is a symmetric and real matrix (thus diagonalizable in R). In order to set ourselves under an asymptotically non-trivial community detection setting, the following assumptions on matrix CΠ need to be posed:

Assumption 1 The symmetric matrix C and the diagonal matrix Π with Tr(Π) = 1 have entries independent of n and satisfy the following hypotheses:

1. (CΠ) pq > 0, ∀ 1 ≤ p, q ≤ k 2. CΠ1 k = c1 k , for some positive c 3. For 1 ≤ p ≤ k, ν p > c Φ .
Assumption 1 is fundamental to our analysis. Its three key requests deserve a detailed motivation and explanation:

1. (CΠ) pq > 0: this implies that there is a non-null probability of connection between any two classes. As a consequence, one may apply the Perron-Frobenius theorem on CΠ: the eigenvalue of CΠ of largest modulus, denoted with c > 0, is positive, simple and its corresponding eigenvector is the only one that can be chosen with all positive entries. By definition of the eigenvalues {ν p } k p=1 , we have in particular c = ν 1 .

2. CΠ1 k = c1 k : this assumption imposes that the Perron-Frobenius eigenvector of CΠ is the vector of all ones, 1 k . From Equation (1) and the independence of the entries of A, one can show that this assumption implies that E[ d], the expected value of the average degree, is equal to c. In the two-class symmetric case, c = (c in + c out )/2 as introduced in Equation (2). Assumption 1 also importantly ensures that the average degree does not depend on the class; in fact, the expected degree of each node i belonging to class p equals θ i (CΠ1 k ) p = θ i c, which is independent of p. this is a standard assumption in the literature and sets the problem under a scenario for which the degree of each node, being independent of its label, does not bring any information on the class structure [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF][START_REF] Decelle | Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications[END_REF][START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF].

3. ν p > c/Φ : this assumption is made for consistency between the actual number of classes k and the number of classes k that can be detected with the spectral methods described in the following and consists in imposing k = k. Recalling that ν 1 = c, this assumption notably implies that cΦ > 1, which is a necessary and sufficient condition to have a giant component in G, as stated in Theorem 1 in Appendix A.

In the case of two classes of equal sizes, the condition ν p > c/Φ is equivalent to setting the problem above the detectability threshold α > α c . As a consequence of this assumption, for all large n with high probability, the spectrum of the affinity matrices under study (B, H r , L rw τ , etc.) can be decomposed as the union of a set of contiguous (non-informative) eigenvalues, collectively referred to as the bulk, and of a set of k (informative) isolated eigenvalues s 1 , . . . , s k found at non-vanishing distance from the bulk (that is, up to multiplicity, s i is isolated if for some > 0, we have for all n and j = i, |s i -s j | > ). As opposed to bulk eigenvalues, the isolated eigenvalues are referred to as informative as their corresponding eigenvectors are non-trivially correlated to the community structure.

It must be noted that, from Assumption 1, we necessarily have ν p > 0 for all 1 ≤ p ≤ k. However, the results of this paper remain valid if one replaces the third point of Assumption 1 by |ν p | > c Φ (that is, ν p for p ≥ 2 may be of arbitrary sign and only its modulus may be lower-bounded; ν 1 = c, being the Perron-Frobenius eigenvalue, is necessarily positive in any case). As the notations to describe this more general case are more cumbersome, we prefer here to focus on the simpler case where ν p > c Φ in the core of the article and defer the discussion of the general case to Section 3.2.2.

Characterization of the Informative Eigenvectors of H r and B

Before stating the main claim of the article, we first need to introduce the values r = ζ p which must be carefully selected to best operate on the Bethe Hessian H r .

Definition 1 (ζ p ) Consider an arbitrary graph G composed of n CC connected components. Let G (j) be the subgraph reduced to the j-th connected component, and H (j) r its associated Bethe-Hessian matrix, for r ∈ R (Equation (3)). Let p be an integer between 1 and n. We define, if it exists, ζ (j) p as:

ζ (j) p =min r≥1 {r : s ↑ p (H (j) r ) = 0}, (4)
where s ↑ p (•) is the p-th smallest eigenvalue. In words, ζ

p is, if it exists, the smallest value of the parameter r ≥ 1 such that the p-th smallest eigenvalue of H (j) r is null.

Here we list a few important properties of these values: j) is the combinatorial Laplacian of the subgraph associated to the j-th connected component. As this subgraph is by definition connected, it is well known [START_REF] Chung | Spectral graph theory[END_REF] q for q > p do not exist.

1. At r = 1, H (j) 1 = D (j) -A (
3. Empirically, on as many connected graphs as we could think of, we have observed that the function s ↑ p (H r ) for r ≥ 1 either never crosses zero (in which case ζ p does not exist), crosses zero exactly twice and is convex between these two crossings (in which case ζ p is the lowest of the two values), or, in very symmetric cases, touches zero exactly once without crossing it (in which case ζ p is that value).

We are now in position to state our main result.

Claim 1 Let G(V, E) be generated according to Equation (1), that is, a DC-SBM with k communities. Let D and A be its degree and adjacency matrices, and H r = (r2 -1)I n +D-rA, for r ∈ R its Bethe-Hessian matrix. Provided that Assumption 1 is satisfied, we have with high probability for large n that:

• There is only one connected component j for which ζ

ζ p =min r>1 {r : s ↑ p (H r ) = 0}. (5)
• The largest p for which ζ p exists is equal to k, the number of underlying communities of the DC-SBM.

2 One has 1 = ζ 1 < ζ 2 ≤ . . . ≤ ζ k ≤ √ cΦ. More precisely: ∀ p s.t. 1 ≤ p ≤ k, ζ p = c ν p + o n (1). (6)
• For 2 ≤ p ≤ k and ζ p ≤ r ≤ √ cΦ, the p smallest eigenvalues of H r are isolated. In particular, zero is an isolated eigenvalue of H ζp . Its corresponding eigenvector x p is correlated to the community structure and the entries of x p are in expectation (over realizations of A) independent of the degrees of the graph.

In simple terms, Claim 1 predicts that, in a graph of k communities, spectral clustering can be successfully performed by successively retrieving the eigenvector associated to the null eigenvalue of each of the matrices H ζ 2 , . . . , H ζ k . Note that these eigenvectors have null entries for nodes that are not connected to the giant component. Thus, all the nodes not connected to the giant component will be arbitrarily classified: this is a structural consequence of sparsity, for which complete recovery is not feasible [START_REF] Mossel | Belief propagation, robust reconstruction and optimal recovery of block models[END_REF].

In the specific case of k = 2 classes of equal size with C i , j = c in if i = j and c out if not, Claim 1 states that the ζ's exist only up to p = k = 2, and that

3 ζ 2 = (c in + c out )/(c in -c out ) + o n (1)
. Also, the relevant community information is found in the 1. This statement allows to define ζp with respect to the Bethe-Hessian matrix of the whole graph, instead of the Bethe-Hessian matrix of its giant component. The fact that Eq. ( 5) is verified with high probability is not evident. In fact, considering all the disconnected components at once might change the ordering of the eigenvalues. More details are to be found in Appendix B. 2. Precisely, we will see that s ↑ k+1 (Hr) > 0 for all r > 1 with high probability, so that ζ k+1 is not defined; in fact it was interestingly shown in [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF] and we will verify that s ↑ k+1 (H √ cΦ ) ↓ 0 + as n → ∞, but the limit is never reached. 3. One can indeed easily show that ν2 = c in -c out eigenvector x 2 associated with the zero eigenvalue (the second smallest) of H r for r = ζ 2 . As we will see next, this value of r differs from the choice r = ρ(B) advocated by [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF]

, unless α = c in -cout √ c = α c = 2 √ Φ , i.e.
, exactly at the phase transition point. Figure 3 provides a visual representation of the typical spectrum of H r for r not too far from a ζ p : the eigenvalue of H r closest to zero is clearly isolated and the eigenvector associated to this eigenvalue is strongly aligned to the community structure. In Figure 4, for a typical realization of a DC-SBM, s ↑ p (H r ) is represented as a function of r for p = 1, 2, 3 and 4 in solid lines and p = 5 in dotted line. In this instance, ζ 5 does not exist (and all subsequent ζ p 's do not exist either), andζ 2 , ζ 3 andζ 4 exist andlie between 1 and2. The non-obvious parts of the claim are (i) of course that such ζ p 's do exist up to p = k, (ii) that they are concentrated around a deterministic value depending on the statistics of the model, (iii) importantly, that zero is indeed an isolated eigenvalue in the spectrum of H ζp , (iv) that the associated eigenvector is informative for the underlying community structure and that is not infected by the degrees of the graph.

ζ 1 = 1,
The statement of Claim 1 is formulated as a claim in the sense that, while efforts have been made to rigorously prove parts of the result (see, e.g., [START_REF] Coste | Eigenvalues of the non-backtracking operator detached from the bulk[END_REF], the mathematical tools required to fully prove Claim 1, to the best of our knowledge, do not exist yet. Instead, the remainder of this section will propose three complementary supporting elements, arising from non-rigorous but convincing approximations, in particular borrowing arguments from the field of statistical physics.

Supporting Arguments to Claim 1

This section details our heuristic arguments in support of Claim 1. Specifically, we will successively show

• under Section 2.3.1 that the vectors x p , solution of H ζp x p = 0, are informative in the sense that they are correlated with the community labels.

• under Section 2.3.2, that the informative null eigenvalue of H ζp (associated to the informative eigenvector x p ) is located in p-th smallest position and is isolated. This result is algorithmically crucial to determine ζ p itself.

• under Section 2.3.3, that the entry x p,i is essentially independent of d i , the degree of node i; precisely, it is strictly independent of d i on average (over random allocation of the labels) and only depends on d i through a noise term of order 1/ √ d i otherwise.

Linearization of the Belief Propagation Equations

We first proceed to our agenda by arguing that the eigenvectors x p of H ζp are correlated with the community labels.

A first approach to the question of community detection in sparse graphs consists in estimating the community allocation probabilities P( i |A) using belief propagation (BP) approximation, based on the local-tree approximation of sparse graphs. The BP node marginal is given by (Decelle et al., 2011, Equation 28), which we recall here:

P( i |A) ≈ π i Z i e -h i j∈∂i j C i , j η ij ( j ), (7)
in which Z i is the probability normalization constant, and the η ji ( i ) (the "messages") and h( i ) are the solutions to (see Decelle et al., 2011, Equations 26,27):

η ji ( i ) = 1 Z ji π i e -h i m∈∂i\j m C m, i η im ( m ) (8)
h i = 1 n k∈V k C i , k P( k |A). (9)
Equation ( 8) admits the trivial η ji ( i ) = π i which corresponds to the so-called paramagnetic fixed point [START_REF] Zhang | Nonbacktracking operator for the ising model and its applications in systems with multiple states[END_REF]. Let us perform a linear expansion of Equation ( 8) around this trivial fixed point by denoting

η ji ( i ) = π i + b ji ( i ). Letting b T = [b T 1 , . . . , b T k ] ∈ R 2|E|k , with b p ∈ R 2|E| containing the entries of b ij (p)
, we obtain as in [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF]:

(T ⊗ B)b = b + o n (1) ( 10
)
where ⊗ is the Kronecker product,

T = ΠC c -Π1 k 1 T k and B ∈ {0, 1} 2|E|×2|E|
is the nonbacktracking matrix, which we recall is defined on the set of directed edges E d of G as:

B (ij)(kl) = δ jk (1 -δ il ), ∀ (ij), (kl) ∈ E d , (11)
in which δ is the Kronecker symbol. Before proceeding, we recall an important property of the spectrum of B rigorously established in [START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF]Gulikers et al., 2017a) and which we will extensively use in the following. Under Assumption 1, with high probability, the matrix B has k isolated real eigenvalues equal to4 ν p Φ + o n (1), 1 ≤ p ≤ k, while all other (complex) eigenvalues are contained within a disk of the complex plane of radius

ρ(B) = √ cΦ + o n (1).
From Equation (10), our interest is in the eigenvalues of (T ⊗ B) equal (or close) to one. From the properties of the Kronecker product, the set of eigenvalues of (T ⊗ B) is

S(T ⊗ B) = {s i (B)s j (T ) : 1 ≤ i ≤ 2|E|, 1 ≤ j ≤ k}.
This induces a one-to-one relation between the k -1 non-zero eigenvalues5 of T and the k -1 eigenvalues of B, that need to satisfy the relation s q (T )s t (B) = 1 for some t, q. From the expression of T , it comes that the p-th largest eigenvalue of T satisfies s ↓ p (T ) =

ν p+1 c for 1 ≤ p ≤ k -1 (see Appendix C for details) so that, for 2 ≤ p ≤ k, there must exist ζ p ∈ S(B) such that: ζ p = c ν p + o n (1). ( 12
)
To these eigenvalues correspond k -1 exact eigenvectors g p ∈ R 2|E| of B satisfying Bg p = ζ p g p (with g p = b p + o n (1)) that are naturally "informative" as they are associated to structural eigenvalues of T . Also, from a belief propagation standpoint, these eigenvectors are small deviations from the uninformative fixed point, so must point towards informative directions. Note importantly here that Equation (12) defines ζ p as a real eigenvalue of the matrix B. So far, this definition needs not correspond to the claimed definition of ζ p as introduced in Claim 1: Section 2.3.2 will show that these two definitions are indeed equivalent. In addition, by Assumption 1 and the Perron-Frobenius theorem, ν 1 = c has unit multiplicity, so that ζ p > 1 + o n (1) for all p > 1. Along with Equation ( 12), we thus have

1 + o n (1) < ζ p ≤ √ cΦ + o n (1), and thus ζ p is asymptotically confined inside the bulk of radius ρ(B) = √ cΦ + o n (1) of B. Figure 2 confirms,
in agreement with the formal result of [START_REF] Coste | Eigenvalues of the non-backtracking operator detached from the bulk[END_REF] obtained in the slightly non sparse regime log(n)/c = o n (1), that inside the disk of radius ρ(B), these are the only informative eigenvalues of B. Specifically, the real eigenvalues of B inside the bulk are divided between (i) the (non-informative) eigenvalues -1, 0, 1, and (ii) the k

eigenvalues γ = ζ p ≈ c νp ≤ √ cΦ just described.
Having identified the k -1 informative eigenvectors g 2≤p≤k ∈ R 2|E| of B, we now consider how to process g p in order to obtain a vector of size n whose entries are in one to one mapping with the nodes of the graph. From the linearization of the marginal probability distribution of the labels (Equation ( 7)), the term j∈∂i b ij ( i ) ≈ j∈∂i g i ,ij ≡ x i ,i brings information on the label of node i, hence on the class structure. Now, let g be the solution of Bg = γg for some γ and g in i = j∈∂i g ij . Then, thanks to the Ihara-Bass formula [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF][START_REF] Terras | Zeta functions of graphs: a stroll through the garden[END_REF], we find that

[(γ 2 -1)I n + D -γA]g in = 0. (13)
The matrix in brackets appears to be the Bethe-Hessian matrix H r = (r 2 -1)I n + D -rA, here for r = γ. This provides an explicit link between B and H r . According to our earlier discussion, the informative eigenvectors of the matrix B have corresponding eigenvalues γ = ζ p and there are consequently k -1 informative eigenvectors x p for 2 ≤ p ≤ k such that

H ζp x p = 0. ( 14
)
Note that, for ζ 1 = 1 we have H 1 = D -A and x 1 is the vector 1 n , which is irrelevant to reconstruct communities.

Let us specifically focus on the particular case of two classes of equal size and consider the question of the detectability threshold. In this case the matrix CΠ has a unique non-trivial eigenvalue equal to ν 2 = (c in -c out )/2. In this setting, B has two informative eigenvalues ≈ ν 2 Φ and ≈ c ν 2 on either side of the disk (bulk) of radius √ cΦ. As the detection problem becomes harder, the two eigenvalues get closer together until they (asymptotically) meet right at the detectability transition where ν 2 Φ = √ cΦ (and α = α c ). Further increasing the detection difficulty (that is, for α < α c ), the two eigenvalues now become complex, each being the complex conjugate of the other. This behavior is shown in Figure 2 (right panel). Summarizing, we have

ν 2 Φ = c ν 2 = √ cΦ, when α = α c ν 2 Φ > √ cΦ > c ν 2 > 1, when α > α c .
The fact that the second eigenvalue of largest amplitude ≈ ν 2 Φ of B remains isolated down to the detectability threshold is the strongest argument in favor of the algorithm proposed by [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF]. The authors in [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF] however ignored the effect on the corresponding eigenvalue ≈ c ν 2 which, from our present discussion, is similar. This behavior can be extended to more than two classes, obtaining ζ p ≈ c νp ≤ ν p Φ with equality at the transition point where ν p Φ = √ cΦ. This can be summarized as follows:

Element 1 Under Assumption 1 for all large n with high probability, the non-backtracking matrix B of a graph G(V, E) generated from a k-class DC-SBM has k -1 isolated real eigenvalues inside the disk of radius √ cΦ (bulk). These eigenvalues are found at positions

1 < ζ p = c νp + o n (1) for 2 ≤ p ≤ k.
Besides, the eigenvectors g p of B associated with these eigenvalues ζ p bring information about the community structure, which can be extracted through the vectors x p ∈ R n , defined as:

x p,i = j∈∂i g p,ij , Bg p = ζ p g p ;
or equivalently satisfying

H ζp x p = 0.
Element 1 provides a first statement of Claim 1, according to which the class information should be retrieved from the vector x p solution to H ζp x p = 0, and that ζ p = c/ν p + o n (1). This argument however does not specify the location (in the ordered list of the eigenvalues) of the null eigenvalue of H ζp to which the eigenvector x p corresponds nor the structure of the eigenvector x p , and in particular its dependence on the degrees of the graph.

The subsequent sections will cover these aspects.

The Excited States of the Ising Hamiltonian on G

In this section, through a statistical physics mapping between the graph G and a system of interacting spins, we aim to informally justify why the smallest eigenvalues of the matrix H r correspond to "informative states" of the system and why the zero eigenvalue of H ζp (associated by definition with the informative eigenvectors x p ) should be isolated and in the p-th smallest position of the spectrum of H ζp .

The Smallest Eigenvalues of H r are Informative

Consider the Ising model of interacting spins σ i ∈ {±1} on the graph G(V, E), in which the intensity of spin coupling is controlled by a parameter r. Its Hamiltonian H(σ) expresses the energy of a given configuration σ as (see e.g. [START_REF] Mezard | Information, physics, and computation[END_REF]:

H(σ) = - (ij)∈E ath 1 r σ i σ j . (15)
The parameter ath(1/r) is related to the temperature of the system and it appears through an inverse hyperbolic tangent for computational ease. The vector σ is a random variable, distributed according to the Boltzmann distribution µ(σ):

µ(σ) = 1 Z e -H(σ) . (16)
Equation ( 16) implies that the configurations with a low energetic cost are realized with high probability. Since σ is a random variable, to properly characterize its stable average configuration, we will focus our study on its statistical mean. The mean m * = σ , as well as the covariance χ * = { σ i σ j } i,j and all subsequent moments of µ (here • denotes an expectation taken over the distribution ( 16)), can be obtained from the explicit expression of Z which, seen as a function of the temperature, is a moment generating function.

(m * , χ * ) ≈ arg min m,χ -log Z Bethe (m, χ) ≡ arg min m,χ F Bethe (m, χ). ( 17
)
For r sufficiently large (that is, at high temperature) the free energy favors disordered configurations with high entropy. In this case, the function F Bethe (m, χ) has a unique Figure 3: Spectrum of the matrix H r for r = 3/2 √ cΦ. On the x axis the eigenvalues, on the y axis the respective histogram. In purple the bulk of uninformative eigenvalues and in orange the two isolated eigenvalues-zoomed to simplify the readability-. The two drawings give a sketchy representation of the ferromagnetic ground state (all black spins) and of the informative state (half black and half white). For this network,

n = 5000, k = 2, π = 1 k /2, c in = 8, c out = 3, r = 6.4, θ i ∼ [U(3, 5)] 4 .
minimum at m * = 0 and the system is said to be in the paramagnetic phase. For smaller (more interesting) values of r, the free energy tends to favor configurations with a small energetic cost (minimizing Equation ( 15)). In this case the minima satisfy m * = 0 and m = 0 becomes a saddle point [START_REF] Leone | Ferromagnetic ordering in graphs with arbitrary degree distribution[END_REF]: the stable average configuration of the spins has a non trivial ordering. In Appendix D we show that, on a graph G with k communities, there exist exactly k directions from the point m = 0 along which the free energy can exhibit a local minimum. The limitation to k directions is a direct consequence of the fact that rank(E[A]) = k. These directions are those along which a non-trivial organization of the spins becomes stable at sufficiently low temperature: they are naturally correlated to the underlying community structure of G.

In order to formally carry out the stability analysis (i.e., to identify the above non-trivial directions of free energy descent), one needs to study the eigenvalues and eigenvectors of the Hessian matrix ∇ 2 m F Bethe of the function F Bethe at the paramagnetic point m = 0. In [START_REF] Watanabe | Graph zeta function in the bethe free energy and loopy belief propagation[END_REF][START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF], the authors show that this Hessian matrix at m = 0 is strictly proportional to the Bethe-Hessian, H r . This induces a natural link to H r , which in the previous section was defined from the non-backtracking matrix B (arising from a linearization of the belief propagation equations). The eigenvectors associated to the negative eigenvalues of the Hessian H r precisely correspond to the sought-for directions towards the local (or global) minima of the Bethe free energy. According to our earlier discussion, only k such directions may exist, so that only up to k eigenvalues (the k smallest) of H r can be negative (i.e., in the limit of large n, s ↑ k+1 (H r ) ≥ 0, ∀ r > 1), and their corresponding eigenvectors are all correlated to the class labels.

We are thus left to showing why specifically the p-th smallest eigenvalue of H r for the particular choice of r = ζ p is of utmost importance, and why it corresponds to the null isolated eigenvalue of H ζp .

Zero is an Isolated Eigenvalue of H ζp

This result follows from the observation, reported in Figure 4, that the k-th smallest eigenvalues s ↑ 1 (H r ), . . . , s ↑ k (H r ) of H r successively equal zero in this order as r increases. Starting from r = 1, H r = H 1 = D -A for which we know that s ↑ j (H 1 ) = 0 for 1 ≤ j ≤ n CC with n CC the number of connected components and s ↑ i (H 1 ) ≥ 0, i > n CC . Increasing r beyond r = 1, the successive smallest eigenvalues of H r first all increase and remain equal to the left edge of the bulk before escaping, each in turn (and in the order s ↑ 2 (H r ), . . . , s ↑ k (H r )), the bulk of non-informative eigenvalues. At their point of escape, they successively shift until they cross zero: this is where the successive values ζ 2 ≤ ... ≤ ζ k are defined. This in particular implies that the p-th smallest eigenvalue of H ζp coincides with (i) the null eigenvalue of H ζp , as well as with (ii) its largest isolated (so the last informative) eigenvalue. This allows us to redefine ζ p as in Equation ( 5) of Claim 1.

Note in passing that, letting r further increase beyond ζ k , the left edge of the noninformative bulk of H r progressively shifts back (from positive values) towards zero until it reaches asymptotically zero in the limit where r = ρ(B) [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF], before increasing again away for r > ρ(B). These last findings may be summarized as follows.

Element 2 Given a graph G(V, E) with k classes, generated from the DC-SBM, the p smallest eigenvalues of the matrix H r are isolated for ζ p ≤ r ≤ √ cΦ, for all 2 ≤ p ≤ k. The entries of the p smallest eigenvectors are correlated with the class labels. Besides, in the specific case where r = ζ p , the p-th smallest eigenvalue of H ζp is equal to zero.

Element 2 answers the problem of locating the informative eigenvalue-eigenvector pair (0, x p ) of H ζp introduced in Element 1, and of justifying why it is isolated. The last step, developed in the next section, consists in better understanding the structural content of the informative eigenvector x p .

Parametrization to Provide Resilience to Degree Heterogeneity

From a purely algebraic standpoint, the Bethe-Hessian matrix H r may be seen as a regularized combinatorial Laplacian. In [START_REF] Dall | Revisiting the bethe-hessian: improved community detection in sparse heterogeneous graphs[END_REF] we studied the problem of a good parametrization r of H r for which the informative eigenvectors are resilient to degree heterogeneity. We briefly report here the main conclusions, based on a local approximation of the neighborhood of each node. See [START_REF] Dall | Revisiting the bethe-hessian: improved community detection in sparse heterogeneous graphs[END_REF], Section 2) for further details.

The argument goes as follows: exploiting sparsity, the graph G can be shown to locally converge to a Galton-Watson tree in which the offsprings are statistically independent [START_REF] Dembo | Gibbs measures and phase transitions on sparse random graphs[END_REF][START_REF] Mossel | Reconstruction and estimation in the planted partition model[END_REF][START_REF] Salez | Some implications of local weak convergence for sparse random graphs[END_REF][START_REF] Mossel | Belief propagation, robust reconstruction and optimal recovery of block models[END_REF][START_REF] Gulikers | An impossibility result for reconstruction in the degree-corrected stochastic block model[END_REF].

Fixing A (and thus the degrees d i ), we now perform a Bayesian analysis on a random allocation of the class labels. Considering node j, labeled as j , as the root of the tree, the probability for offspring i to have label i reads

P( i | j , A ij = 1) = P( i , j |A ij = 1) P( j |A ij = 1) = dθ i dθ j P( i , j , θ i θ j |A ij = 1) P( j ) = dθ i dθ j P(A ij = 1|θ i , θ j , i , j )P( i )P( j )P(θ i )P(θ j ) Zπ j = π i C i , j c = (ΠC) i , j c = (CΠ) j , i c . (18)
Recalling the definition of the eigenvalue-eigenvector pairs (ν p , v p ) of CΠ (i.e., CΠv p = ν p v p ), let us define u p ∈ R n with u p,i = v p, i the n-dimensional class-wise expansion of v p : this vector u p is inherently random as the class allocations 1 , . . . , n are here considered random. Then, for n large, under the limiting tree approximation with conditionally independent offsprings, we take an expectation over the random allocation of the labels different from i ( \i ), with A and i known:

E \i [(Au p ) i |A, i ] = j∈∂i E \i [u p,j |A, i ] = j∈∂i E \i [v p, j |A, i ] ≈ j∈∂i j P( j | i , A ij = 1)v p, j = 1 c j∈∂i j (CΠ) i , j v p, j = 1 c j∈∂i (CΠv p ) i = d i c ν p v p, i = d i ν p c u p,i (19)
where the approximation follows from the fact that conditional independence of the neighbors of a given node only holds asymptotically; in particular, the approximation becomes an equality on a tree rooted at node i. As a consequence of Equation ( 19), using H r = (r 2 -1) + D -rA, we find that

E \i [H r u p |A, i ] ≈ (r 2 -1)I n + D 1 -r ν p c u p . (20)
In order for this equation to be an approximate eigenvector equation for arbitrary degrees in D, the right hand-side term proportional to D must vanish. That is, one must select

r = c ν p ≈ ζ p
(this last approximation having been introduced and discussed in Section 2.3.1). This result implies that the eigenvectors x p defined as H ζp x p = 0, for 1 ≤ p ≤ k, correspond to a noisy version of u p which is not affected on average over the class allocation by the degree distribution. In [START_REF] Dall | Revisiting the bethe-hessian: improved community detection in sparse heterogeneous graphs[END_REF], we confirmed -so long that the average node degree is not too small-that the approximation holds, beyond the average, for every typical realization of the class allocation. In a nutshell, this behavior unfolds from the following remark: denoting ξ p the "noise" vector satisfying x p = u p + ξ p , and thus

H ζp (u p + ξ p ) = 0, Var \i [(H ζp u p ) i |A, i ] = Var \i [(H ζp ξ p ) i |A, i ] = Var \i [(ζ 2 p -1 + d i )ξ p,i -ζ p (Aξ p,i )|A, i ] ≈ O d i (d 2 i )Var \i [ξ p,i |A, i ] + O d i (d i )Var \i [ξ p,i |A, i ] ≈ O d i (d 2 i )Var \i [ξ p,i |A, i ] (21)
where we recall that the underlying random variable is the random class allocation. There we used the fact that the random variables ξ p,i and [Aξ p ] i are essentially independent and that the variance of the sum of the d i asymptotically independent variables ξ p,1 , .

) i |A, i ] = O d i (d i ).
Combining with Equation ( 21), we get that, for

d i sufficiently large, Var \i [ξ p,i |A, i ] = O d i (d -1 i ).
The vector x p can therefore be written as the sum of the deterministic information and a noise with amplitude inversely proportional to the square root of the degree, consistently predicting that nodes with higher degrees are easier to classify.

These results are summarized under the form of our last argument.

Element 3 The eigenvector x p (1 ≤ p ≤ k), solution to H ζp x p = 0, is a noisy version of the vector u p , defined as u p,i = v p, i , where CΠv p = ν p v p , where the noise for entry i scales as 1/ √ d i for d i sufficiently large and is zero on average. Consequently, the entries of x p do not, to first order, depend on the degree distribution but only on the labels.

Element 3 combined with Elements 1, 2 completes our justification of Claim 1.

Comments and Performance Comparison

Let us now discuss how Claim 1 can be exploited in practice to obtain an efficient spectral clustering algorithm for sparse graphs with a heterogeneous degree distribution.

From Claim 1 and our detailed analysis in Sections 2.3.1-2.3.3, for all large n with high probability, under Assumption 1, the k -1 Bethe-Hessian matrices H ζ 2 , . . . , H ζ k have an eigenvalue equal to zero, which is isolated. The corresponding eigenvectors x 2 , . . . , x k (which are not necessarily orthogonal as they correspond to distinct matrices) are all informative in the sense that they are noisy realizations of piece-wise constant versions of the eigenvectors v 2 , . . . , v k of the matrix CΠ (each "piece" identifying each class). Besides, x 2 , . . . , x k are, to first order, resilient to the degree distribution in the network.

This so far assumes that the number of classes k is known, and that all classes satisfy the separability condition of Assumption 1. Yet, in practice, k is generally unknown. The following remark provides a consistent estimator for k, as proposed in [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF][START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF].

ζ p ← r such that s ↑ p (H r ) = 0 on 1 < r ≤ √ cΦ 6: X •,p ← x p such that H ζp x p = 0 7: Estimate community labels ˆ from the node embedding X = [X •,1 , . . . , X •, k].
8: return Estimated number k of communities and label vector ˆ .

Remark 1 (Estimation of the number of classes) Under Assumption 1, for all large n with high probability, kB = kH = k where

kB = i, R[s i (B)] > √ cΦ kH = i, s i (H r )| r= √ cΦ < 0 ,
where R indicates the real part.

The remark is justified by the fact that r = √ cΦ is the (limiting) right edge of the bulk of the non-backtracking matrix B and that the eigenvalues of B whose real part exceeds √ cΦ (in the limit n → ∞) are the isolated largest real eigenvalues of B. These are mapped, according to Figure 4 and the discussion of Section 2.3.2, to the number of communities. Likewise, from the connection between H r and B, for r = √ cΦ, the negative eigenvalues of H r are one-to-one mapped to the largest real eigenvalues of B.

The arguments above naturally lead to Algorithm 1 for community detection using the Bethe-Hessian matrix. Algorithm 1 is a meta-algorithm, written for the user readability. Section 4 will provide an efficient implementation (Algorithm 2) of Algorithm 1 accounting for deeper algorithmic considerations.

Applying Algorithm 1, Figure 5 compares the overlap performance (defined in Equation ( 22)), achieved by Algorithm 1 versus competing methods for synthetic DC-SBM graphs as a function of the hardness of the problem α.

Overlap = max ˆ ∈P 1 1 -1 k n i=1 δ( ˆ i , i ) - 1 k , ( 22
)
where ˆ is the estimate of the label eigenvector and P the set of the permutations of ˆ .

On the left hand-side of Figure 5 is depicted a symmetric two-class scenario, for which α is defined as per Equation (2) and community reconstruction is asymptotically feasible if and only if α > α c . On the right hand-side a more involved multiple class scenario is devised: the matrix C is obtained here by drawing random Gaussian numbers with mean c out and variance f = c out /k, before being symmetrized and tuning the diagonal elements to guarantee CΠ1 k = c1 k . In this case, we still define α = 2(c -c out )/ √ c but α c now only represents the transition conjectured in [START_REF] Decelle | Asymptotic analysis of the stochastic block model for modular networks and its algorithmic applications[END_REF] of the planted partition model Blue sharp diamonds are the result given by Algorithm 1, red diamonds the algorithm using the Bethe-Hessian of [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF], green squares are for the non-backtracking of [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF], yellow pentagons the adjacency matrix, purple stars the random walk Laplacian, black dots are the algorithm of [START_REF] Qin | Regularized spectral clustering under the degree-corrected stochastic blockmodel[END_REF]. For both graphs n = 50.000,

θ i ∼ [U(3, 10)] 4 , π ∝ 1 k , c = 5. (Left): k = 2, c out : 0.05 → 4. (Right): k = 5, π p ∼ N (k -1 , (4k) -1 ), C a>b ∼ N (c out , c out k -1 ), c out : 0.05 → 4. α := 2(c -c out )/ √ c and α c = 2/ √ Φ.
Averages over 10 samples.

for which f = 0. In practice, in our simulation with f = 0, there are multiple transitions that are close (but in general different) to α c .

In both cases, Algorithm 1 outperforms the competing algorithms. For k = 2, as α → α c , Algorithm 1, the algorithm of [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF] based on H √ cΦ and the one of [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF] based on B give essentially the same result, which confirms that they are indeed equivalent at the phase transition. For easier detection problems though, except for the algorithm of [START_REF] Qin | Regularized spectral clustering under the degree-corrected stochastic blockmodel[END_REF] which is only slightly less accurate, the performance of all methods is largely improved by Algorithm 1. Interestingly, it can be shown here that

| √ cΦ -ζ 2 | < |ν 2 Φ -ζ 2 |
, giving an intuition on why the Bethe-Hessian method of [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF] performs better than the non-backtracking approach of [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF]. As for the standard spectral clustering algorithm which exploits the dominant eigenvectors of A, it can only perform non trivial community reconstruction for α far beyond the threshold α = α c . The algorithm of [START_REF] Shi | Normalized cuts and image segmentation[END_REF], based on the popular random walk Laplacian L rw , is here incapable of making any non-trivial reconstruction for the considered set of parameters, suggesting that its dominant eigenvectors are not informative. This last point is quite interesting and raises the possibility that the informative eigenvectors of L rw are possibly lost in the bulk. The next section will confirm this affirmation and make the statement more precise.

Relation Between H r and the Regularized Random Walk Laplacian L rw τ

In this section we show that there is a strong connection between the Bethe-Hessian matrix and the regularized random walk Laplacian matrix L rw τ = D -1 τ A, where D τ = D + τ I n . We will specifically prove that, when τ is carefully chosen, the matrix L rw τ has isolated informative eigenvalues in its largest positions. Remarking that L rw τ has the same eigenvalues as the matrix

L sym τ = D -1/2 τ AD -1/2 τ
, we also find a strong connection between L rw τ and the community detection method of [START_REF] Qin | Regularized spectral clustering under the degree-corrected stochastic blockmodel[END_REF].

The mapping between H r and L rw τ unfolds from the following basic remark: for ζ p > 1 (as defined in Equation ( 5)),

[(ζ 2 p -1)I n + D -ζ p A]x p = 0 (23) ⇔ [D + (ζ 2 p -1)I n ] -1 Ax p = 1 ζ p x p (24)
inducing a natural mapping between the Bethe-Hessian matrix at r = ζ p to the regularized random walk Laplacian at τ = ζ 2 p -1. In particular, the vector x p solution of Equation ( 23) coincides with the eigenvector of L rw ζ 2 p -1 associated with the eigenvalue 1/ζ p . In the following this relation is further analyzed and will be used to provide an alternative method to detect communities down to the detectability threshold.

Main Result

The spectral methods of [START_REF] Shi | Normalized cuts and image segmentation[END_REF][START_REF] Andrew | On spectral clustering: Analysis and an algorithm[END_REF] do not lead to good partitions in the sparse regime, as evidenced in Figure 5. The reason, which we briefly commented in Section 1, is that the bulk of uninformative eigenvalues of both L rw and L sym undergo a certain spreading in the sparse regime (see Figure 1B) which "swallows" the informative eigenvalues within the bulk. This is not to say that informative eigenvectors vanish: rather, these eigenvectors are associated to non-isolated eigenvalues; this was already observed in [START_REF] Joseph | Impact of regularization on spectral clustering[END_REF] for the benchmark network Political blogs [START_REF] Lada | The political blogosphere and the 2004 us election: divided they blog[END_REF]. This has two major negative consequences: (i) it may be practically infeasible (especially as the clustering task is more difficult) to identify the correct informative eigenvalue, and (ii) since the eigenvectors associated to close-by eigenvalues (in the bulk, the typical distance between consecutive eigenvalues is O n (1/n)) tend to "spread" across the eigenvectors of these neighboring eigenvalues, not a single eigenvector but a collection of neighboring eigenvectors need to be considered. This is clearly deleterious to spectral clustering.

This problem is partially solved through regularization through a parameter τ (added to D as D τ = D + τ I n ): it was indeed shown in a series of related works that regularization helps clustering in sparse networks [START_REF] Lei | Consistency of spectral clustering in stochastic block models[END_REF][START_REF] Joseph | Impact of regularization on spectral clustering[END_REF][START_REF] Qin | Regularized spectral clustering under the degree-corrected stochastic blockmodel[END_REF]. In all these works, the results obtained (or inferred) are not straightforwardly applicable to the sparse regime. Furthermore, according to the analysis performed in these works (however only far from the transition point), large values of τ would seem preferable; however, in practice, it is rather observed that small values entail better performances.

We answer here to the question of why small parametrizations τ should be preferred to large τ , as well as which is the smallest value of τ for which we can guarantee the existence of isolated informative eigenvectors of L rw τ , down to the detectability threshold. Besides, in accordance to the previous sections, we determine a value for τ which is optimal in its (i) ensuring the existence of these informative eigenvectors all the way to the detectability threshold, and most importantly (ii) which is resilient to degree heterogeneity in the graph. Formally, the result is formulated as follows.

Claim 2 Consider the graph G(V, E), built on a sparse DC-SBM as per Equation (1) with k communities. Let ζ p for 2 ≤ p ≤ k be defined as per Equation (5) (imposing ζ 1 = 1) and τ ∈ R be such that ζ 2 p -1 ≤ τ ≤ cΦ -1. Then, under Assumption 1, for all large n with high probability, the p largest eigenvalues of the matrix L rw τ are isolated. In particular,

s ↓ p (L rw ζ 2 p -1 ) = 1 ζ p .
Note that the proposition is not an obvious consequence of the equivalence between Equation (23) and Equation (24) as it is not clear that the eigenvalue 1/ζ p of L rw ζ 2 p -1 corresponds to the p-th largest and that it remains isolated (as is zero in the spectrum of H ζp ).

Since Claim 2 asserts that the eigenvalue 1/ζ p of L rw ζ 2 p -1 is the p-th largest and is isolated, it can be clearly identified for all finite but large n: its corresponding informative eigenvector x p (from Equation 24) can thus not be confused with other eigenvectors. This eigenvector x p is also the solution to H ζp x p = 0 (from Equation ( 23)); its properties have been extensively discussed in Section 2: it is in particular asymptotically insensitive to the degrees of the graph. Intuitively, by picking τ away from ζ 2 p -1, the entries of the p-th eigenvector are likely to be more polluted by the degrees of the network. This suggests that large values of τ (as studied by the authors of [START_REF] Joseph | Impact of regularization on spectral clustering[END_REF] are likely to lead to sub-optimal partitions.

Claim 2 further asserts that, for any τ in the interval [ζ 2 p -1, cΦ -1], the p dominant eigenvalues of L rw τ are isolated. Since in addition ζ p ≤ ζ k ≤ √ cΦ, regardless of the hardness of the detection problem, the k largest eigenvalues of L rw cΦ-1 (i.e., L rw r 2 -1 for r = √ cΦ) must be isolated. Note in passing that, in [START_REF] Qin | Regularized spectral clustering under the degree-corrected stochastic blockmodel[END_REF], the authors propose the regularization τ = d = 1 n i d i from a heuristic intuition. While suboptimal according to the claim, this choice is somewhat meaningful as τ must indeed essentially grow with the average degree d ≈ c.

Besides, as a corollary, this provides a convenient alternative method to estimate the number of communities, based on the regularized Laplacian matrix.

Remark 2 (Second estimation of the number of classes) As a direct consequence of Claim 2, we have that, for all large n with high probability k = kB = kH = kL , where

kL = i : s i (L rw τ )| τ =cΦ-1 ≥ 1 √ cΦ .
Note that, since S(L rw τ ) = S(L sym τ ), this calculation can be performed on a symmetric matrix, bringing gain in computational efficiency.

The technical supporting elements of Claim 2 can be found in Appendix E. In a nutshell, letting τ (r) = -s ↑ p (D -rA), we show that there is a one-to-one map between the isolated eigenvalues of L rw τ (r) and of H r , only reversed in order (the smallest of H r are mapped to the largest of L rw τ (r) ). This is particularly valid for r = ζ p for which -s ↑ p (D -ζ p A) = ζ 2 p -1. We further show that the function τ (r) is bijective, thereby extending the results of Claim 1 for r

∈ [ζ p , √ cΦ] into τ ∈ [ζ 2 p -1, cΦ -1].

Side Comments

This section provides further interpretations and insights of the results introduced in the previous chapters. Firstly, we relate Algorithm 1 to other commonly adopted methods for spectral clustering and secondly we discuss how our result can be extended in the presence of disassortative communities.

3.2.1 Connection Among the Spectral Algorithms Based on B, H r , L rw τ , L rw , L We showed in the previous sections the deep connection between the belief propagation equations and the regularized Laplacian matrix L rw τ , by successively passing through the non-backtracking B and the Bethe-Hessian matrices H r , so far treated in parallel (and with different tools) in the literature. From a practical perspective, we notably pointed that adding the regularization τ I n to the degree matrix D (i) favors efficient clustering in sparse networks but (ii) is optimally tuned for τ taken as a function of the hardness of the detection task (explicitly, for k = 2, a function of α = (c in -c out )/ √ c). In particular, for more challenging clustering tasks, larger values of r in H r and τ in L rw τ should be employed. On the opposite, in the limit of trivially simple clustering, r → 1 (so that H r → D -A, the standard Laplacian) and τ → 0 (so L rw τ → L rw , the standard random walk Laplacian). This suggests that the classical Laplacian matrices are still "optimal" in the sparse regime, yet only for easy clustering tasks, i.e., possibly far beyond the detectability threshold.

Most importantly, while the aforementioned "difficulty" of a clustering task is of course not accessible to the practitioner, we showed that the optimal tuning of the hyperparameters r and τ can be practically obtained by retrieving isolated eigenvalues in the spectra of H r , B, or L rw τ (through the fundamental quantities ζ 2 , . . . , ζ k ). In passing, the relative distance of the ζ p 's to the bulk of non-informative eigenvalues is a further clue for the practitioner of the level of confidence of the ultimate clustering result.

Extension to Disassortative Communities

Another remark concerns the possible existence of disassortative communities in the graph, i.e., groups of nodes which are identified as a class because they repel (rather than attract) each other. As a concrete example one may think of the vertices of a graph as the words contained in a text with edges if two words appear next to each other: on this graph, adjectives and nouns represent two disassortative communities [START_REF] Mark | Finding and evaluating community structure in networks[END_REF]. Under a two-class DC-SBM model, a disassortative network can be easily generated by imposing c out > c in . In this case, the second largest eigenvalue of CΠ will therefore be negative and equal to ν 2 = (c in -c out )/2 < 0.

The problem of detection of two disassortative communities can be addressed both with the non-backtracking and the Bethe Hessian matrices (as shown in [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF][START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF]. In particular, in the presence of two disassortative communities, s matrix H r can detect two disassortative communities for r < 0, becoming a deformed version (equal in the limit of easy clustering) of the signless Laplacian D + A, frequently used to study bipartite graphs [START_REF] Cvetković | Signless laplacians of finite graphs[END_REF].

↓ |•| 2 (B) < 0,
While defining communities when both assortative and disassortative elements are present is debatable, at least the result of Claim 2 can be extended to the case where the matrix CΠ has possibly both positive and negative eigenvalues, as follows.

Remark 3 (Eigenvalues of CΠ of arbitrary sign) Assume, unlike stipulated in Assumption 1, that the sign of the eigenvalues ν p of CΠ is arbitrary (except for ν 1 = c), but that |ν p | > c Φ in absolute value. We sort the eigenvalues ν p according to their modulus as

ν 1 > |ν 2 | ≥ . . . ≥ |ν k |.
In this setting, Claim 2 generalizes by stating that the p eigenvalues of L rw τ with largest modulus are isolated and that s

↓ |•| p (L rw ζ 2 p -1 ) = 1 ζ p = ν p c + o n (1). ( 25
)
As suggested above, while the mathematical result is well defined for ν p < 0, the interpretation of the negative eigenvalues of CΠ is not straightforward for k > 2 and the very definition of a community might not be evident. As an example, the right frame of Figure 6 displays a network designed from the SBM with k = 4, Π ∝ I k and C defined as:

C = c in c out c out c in ⊗ c out c in c in c out (26)
with '⊗' the Kronecker product. This graph has a hierarchical structure: it can be divided into two assortative communities, each of them composed of two disassortative communities.

In Figure 6 (right) the position of each node is assigned according to its true label, while the color is assigned according to the output of k-means on the vectors x p of L rw

ζ 2
p -1 with largest modulus. We can observe that our algorithm performs well also in this particular case. For the sake of clarity, we close here the discussion on the generalization to cases where ν p may be negative and, unless otherwise stated, we are thus once again under Assumption 1.

Implementation of the Algorithm

Based on the various discussions and results from the previous section, we now formally introduce our proposed algorithm along with pragmatic discussions on implementation cost, optimization and robustness to real-world network configurations.

Applicability to Real Networks

As opposed to the sparse stochastic block model, the DC-SBM model of Equation ( 1) accounts both for the sparsity and heterogeneity of real networks. Yet, it does not capture other fundamental aspects of real graphs. For instance, the presence of many triangle in social networks [START_REF] Paul | Transitivity in structural models of small groups[END_REF]) and more generally of short loops may invalidate the local tree-like approximation that underlies our theoretical findings. Nevertheless, its usability on real data sets suggest that some properties of B are valid for various graph topologies. In particular, the following properties of B seem to hold in general:

• All complex eigenvalues come in pairs of complex conjugates and most of them are bounded by a circle on the complex plane of radius r radius ≈ ρ(B). This statement originates both from empirical observations and from heuristic arguments on the asymptotic density of the eigenvalues of B dicussed in [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF]. • The number of real eigenvalues of B, different from {-1, 0, 1, ρ(B)}, is even. Half of the eigenvalues are larger in modulus than r radius , while the other half lies between 1 and r radius . All of them are isolated.

It is easy to verify that the two earlier properties hold for any graph topology if D = cI n and can be extended to the case in which D ≈ cI n as suggested in [START_REF] Coste | Eigenvalues of the non-backtracking operator detached from the bulk[END_REF]. Given the connection between B and the Bethe Hessian matrix H r , from these two points, we can claim that the steps of Algorithm 1 are all well defined on arbitrary graphs.

Some properties of B are however likely model-dependent and thus not resilient to arbitrary graphs, such as • ρ(B) may be far from i d 2 i i d i , which is a natural estimator for cΦ in the DC-SBM model.

• For an arbitrary network, ζ p (as defined in Claim 1) may be far from ρ(B)

s ↓ R,p (B)
in general.

These observations impose that, when devising a spectral clustering algorithm adapted to real networks, these purely DC-SBM considerations should not be exploited.

A further specificity of real graphs is their possibility to be inherently made of disjoint components. We have instead so far worked with the assumption that G has a giant component and that the communities are contained within the giant component. In practice, communities may live in disconnected subgraphs (in a two-class DC-SBM, this would correspond to setting c out = 0). In this case, one can perform a two-stage clustering. First the connected components are detected, looking into the eigenvectors with zero eigenvalue of H ζ 1 = D -A. Afterwards, each connected component is treated independently. We will thus consider, without loss of generality, that the real graphs which we will consider are connected.

Estimation of k

The problem of estimating the number of communities in an unsupervised manner is in general non-trivial. Methods based on the non-backtracking rather than Bethe-Hessian matrix have been studied and exploited to efficiently recover communities regardless of the generative model [START_REF] Can | Estimating the number of communities in networks by spectral methods[END_REF]. In the course of our argumentation, we suggested different ways to estimate the number of communities, which are asymptotically equivalent (recall Remark 1 and Remark 2):

1. kB = |{s R,i (B) > ρ(B)}| 2. kH = |{s i H √ ρ(B) < 0}| 3. kL = s i L sym ρ(B)-1 > 1 √ ρ(B)
.

From a computational standpoint, note that the eigenvalues of the matrix B different from ±1 are the same as the eigenvalues of the matrix B ∈ R 2n×2n [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF][START_REF] Coste | Eigenvalues of the non-backtracking operator detached from the bulk[END_REF], defined as

B = A I n -D I n 0 ( 27
)
so all computation involving the eigenvalues of B can be performed in an efficient manner on the matrix B . All three estimators above require the value of ρ(B) that can be obtained efficiently by direct computation. But not all estimates perform the same in practice.

Estimator 1 is in particular not efficient as not only real but also complex eigenvalues of B need be evaluated, which sensibly slows down the algorithm. Estimators 2 and 3 do not suffer this limitation. Estimator 3 is nonetheless preferred as estimating eigenvalues with largest, rather than smallest, algebraic value can in general be performed more efficiently. Furthermore, when the matrix CΠ has at the same time both positive and negative eigenvalues, all communities can be detected using the same matrix L sym ρ(B)-1 , as opposed to the Bethe-Hessian for which one needs to consider both H √ ρ(B) and H - √ ρ(B) . Consequently, the estimator kL is the one that will be adopted in our final Algorithm 2. The complexity of estimating k scales as O(n k3 ), as computing the p largest eigenvalues of a sparse matrix of size n × n (that is, containing O(n) non-null elements) costs O(np 2 ) with state-of-the-art methods such as restarted Arnoldi methods (see for instance [START_REF] Saad | Numerical Methods for Large Eigenvalue Problems[END_REF].

Subroutine 1 details how to estimate 6 k and it will be referenced in Algorithm 2.

6. Subroutine 1 is actually a naïve implementation to estimate the number of communities, as it requires to compute several times the same eigenvectors. In the CoDeBetHe.jl package a more efficient implementation is proposed, based random projections, that allows to efficiently estimate k.

Subroutine 1 estimate_number_of_classes 1: Input : adjacency matrix of a connected, undirected graph G 2: Compute ρ(B) ← s

↓ |•| 1 (B ), with B defined in Equation (27) 3: j ← 1 4: while v > 0 do 5: v ← s ↓ j+1 (L sym ρ(B)-1 ) -1 √ ρ(B)
6:

j ← j + 1 7: k ← j 8: return Estimated number k of communities.
Of course, as per Assumption 1 in the DC-SBM setting, since for all 1 ≤ p ≤ k, ν p > c Φ , all communities were claimed "visible". In practice though, this becomes a stringent condition. If, in particular, only κ < k eigenvalues of CΠ exceed the threshold, fewer eigenvectors will be exploited and fewer classes will be looked for by the algorithm (as is also the case of the algorithms in [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF][START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF]. In particular, k -κ communities might be merged to close-by communities or spread across the κ detected communities.

Estimation of {ζ p }

From Section 2.3.1, the values of the ζ p for 1 ≤ p ≤ k can be estimated as

ζ p = c ν p + o n (1) = ρ(B) s ↓ p (B) + o n (1). (28)
This estimation of ζ p (via B ) is computationally efficient but may be quite inaccurate for graphs not generated from the DC-SBM model. Conversely, a naive line search for ζ p ∈ (1, ρ(B)) satisfying s ↑ p (H ζp ) = 0 is computationally inefficient. We propose here a faster method, motivated by a Courant-Fischer theorem argument. The details and proof of convergence of the proposed method are given in Appendix F.

In a few words, starting from an initial guess r 0 , we devise an iterative sequence r 0 , r 1 , . . . , such that r t → ζ p . In Appendix F we show that the convergence is guaranteed when setting r 0 = ζ p+1 , under the convention ζ k+1 = ρ(B). The values of ζ p are then estimated from the largest, ζ k , to the smallest, ζ 2 .

The algorithm builds on two parts: (i) the computation of a p × p matrix obtained from the eigenvectors of H rt with computational cost of O(np 2 ), (ii) a subsequent line-search using this matrix, with computational cost O(p 3 ). The advantage of this method is that the line search (which requires many iterations) is computationally cheap, while the most expensive part of the algorithm needs to be performed much fewer times with respect to the greedy line-search to obtain the same accuracy. The total complexity of the algorithm needed to compute the vector ζ = (ζ 1 , . . . , ζ k ) T , scales as O(nk 3 ).

The proposed algorithm is described in the following subroutine. We indicate with X r ∈ R n×p the matrix containing in its columns the p smallest eigenvectors of H r and S r ∈ R p×p the diagonal matrix with the corresponding eigenvalues. For further details, the reader is referred to Appendix F.

Subroutine 2 compute_ζ

1: Input : adjacency matrix of a connected, undirected graph G, number of classes k

2: r ← s ↓ |•| 1 (B )
, where B is defined in Equation ( 27)

3: p ← k 4: while p > 1 do 5: repeat 6: Compute S r = diag s ↑ 1 (H r ), . . . , s ↑ p (H r ) and X r ∈ R n×p , with H r X r = X r S r 7: Λ r ← X T r DX r ∈ R p×p (D ∈ N n×n , the diagonal degree matrix) 8: r * ← r : s ↓ 1 (r S r + (r -r )Λ r ) = (r -r )(1 + rr ), with line-search on r ∈ (1, r) 9:
r ← r * 10:

until convergence 11: δ ← multiplicity of s ↑ p (H r ) = 0 12: Set ζ p , ζ p-1 , • • • , ζ p-δ+1 to r 13: p ← p -δ return Vector ζ = (ζ 1 , ζ 2 , . . . , ζ k ) T .
Note that, although the subroutine compute_ζ only outputs the vector ζ, it can also be used to directly compute the informative eigenvectors {x p }.

Figure 7 provides a typical output of the computation of ζ and confirms the accuracy of the proposed algorithm. On the top line the algorithm is tested on a network created from the DC-SBM model, for which ρ(B)/s ↓ R,p (B) is a valid estimator for ζ p . The horizontal line indicates ρ(B) = √ cΦ which is the upper bound of ζ p . In our simulations, C is generated randomly (see the caption of Figure 4), for the first two points we obtained ν 3 < √ cΦ, invalidating Assumption 1. In this case we see that c/ν p > √ cΦ and the corresponding estimated value of ζ p saturates at √ cΦ. On the contrary, whenever Assumption 1 is verified, the estimate of ζ p is correct.

In the bottom line we compare the two methods to estimate ζ p on the three real networks that clearly show that these two methods are different on graphs not generated from the DC-SBM. Given that Subroutine 2 provides a direct computation of the eigenvalues of B, it should be preferred.

For illustrative purposes, we tested the execution time of Subroutine 2 for a SBM network with k = 4 classes of equal size, n = 10 5 , c in = 26.20 on the diagonal elements of C and c out = 4.60 for the off diagonal elements. The values of ζ p were computed to machine precision in approximately 40 seconds on a standard laptop using our Python implementation, while in approximately 5 seconds with the CoDeBetHe.jl Julia implementation. The complexity of this subroutine scales linearly with n, and thus can be applied to large networks but cubically with respect to k (and not quadratically as usual in the spectral clustering context) decreasing the computational efficiency when a rather large number of classes is present.

Projection of the Embedded Points on a Hyper-Sphere

The study performed so far identifies the presence of k informative eigenvalues and describes the content of the associated k eigenvectors X ≡ [x 1 , . . . , x k ] ∈ R n×k with x p the eigenvector [START_REF] Leskovec | SNAP Datasets: Stanford large network dataset collection[END_REF]. The value of k is estimated according to Subroutine 1. associated to the p-th smallest eigenvalue of H ζp (in particular, the vectors x p 's need not be orthogonal). The rows X 1,• , . . . , X n,• ∈ R k of the matrix X form a k-dimensional feature for every node, which are used in a last small-dimensional clustering step, usually employing the k-means algorithm. The fact that k-means is particularly efficient when the low-dimensional clusters are quite "isotropic" strongly motivates the need for the entries of the vectors x p not to be affected by the node degrees (which would otherwise spread the clusters unevenly).

Yet, to further tackle residual degree dependence, a classical method, prior to k-means, consists in normalizing all vectors X i,• to X i,• = 1 (this is Step 4 of Algorithm 0) . This method is motivated by the assumption that the degree dependence in each X ij is separable from the label dependence, a fact that is verified in sufficiently dense DC-SBM networks [START_REF] Jin | Fast community detection by score[END_REF] and to some extent also in sparser graphs [START_REF] Qin | Regularized spectral clustering under the degree-corrected stochastic blockmodel[END_REF]. Besides, under this normalization, the k-means algorithm is restricted to the unitary hypersphere, improving its convergence to the genuine solution, especially when in presence of many communities.

Algorithm 2 Community Detection in sparse and heterogeneous graphs 1: Input : adjacency matrix of a connected, undirected graph G 2: k ← estimate_number_of_classes (Subroutine 1)

3: ζ ← compute_ζ (Subroutine 2) 4: Initialize X ∈ R n× k 5: for p = 1 : k do 6: X •,p ← x p ∈ R n
, where x p is the eigenvector with eigenvalue s ↑ p (H ζp ) = 0. 7: Normalize the rows of X i,• ← X i,• / X i,• 8: Estimate community labels ˆ as output of k-class k-means on the rows of X. 9: return Estimated number k of communities and label vector ˆ .

As such, while our proposed algorithm naturally discards degree dependence in the entries of X under the DC-SBM setting, the reality of practical networks may disrupt this expected behavior and the projection of the vectors X i,• on the unit hyper-sphere both alleviates this deleterious effect and further improves the convergence of k-means. We thus adopt this normalization step in our final Algorithm 2 and will confirm its practical gains when clustering real graphs.

The total theoretical complexity of Alg. 2 is dominated by subroutines 1 and 2, that both run in O(n k3 ), as the k-means step costs only O(n k2 ). This is to compare with the usual complexity of spectral clustering algorithms on sparse graphs that are in O(n k2 ) (see for instance [START_REF] Tremblay | Approximating Spectral Clustering via Sampling: A Review[END_REF]. This additional cost comes with a better classification performance in many real-world graphs, as discussed in the next section. In practice, to give an order of magnitude of computation times 7 , running Algorithm 2 on a SBM 8 with k = 4 classes of equal size, n = 10 5 (resp. 10 6 ), c in = 26.20 on the diagonal elements of C and c out = 4.60 for the off diagonal elements, takes approximately 85 (resp., 2200) seconds on a laptop. If k is known in advance, this times drop to approximately 45 (resp., 700) seconds.

Algorithm and Performance on Real Networks

In this section we compare the performance of Algorithm 2 versus competing spectral methods on real-world networks that do not have ground-truth label assignment as well as LFR synthetic benchmark networks.

Measuring the quality of an inferred partition ˆ is in general not straightforward since communities are not uniquely defined. Two different scores are then adopted adopted. One is the modularity 9 (Newman and [START_REF] Mark | Finding and evaluating community structure in networks[END_REF][START_REF] Newman | Modularity and community structure in networks[END_REF], M:

M = 1 2|E| n i,j=1 A ij - d i d j 2|E| δ( ˆ i , ˆ j ). (29)
7. The laptop's RAM is 7.7Gb with Intel Core i7-6600U CPU @ 2.6GHz x 4 8. These times refer to our Python implementation, while the CoDeBetHe.jl implementation runs for n = 10 5 (resp., 10 6 ) in approximately 15 (resp., 450) seconds if k is not known and in approximately 6 (resp., 100) seconds if k is known a priori. 9. Note that the measure of the modularity is meaningful on assortative or disassortative networks but not on "hybrid" networks for which a more involved description would be needed. [START_REF] Lancichinetti | Benchmark graphs for testing community detection algorithms[END_REF] benchmark graphs as a function of the mixing parameter µ. Blue sharp diamonds are the result given by Algorithm 2, red diamonds the algorithm using the Bethe-Hessian of [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF], green squares are for the non-backtracking of [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF], yellow pentagons the adjacency matrix, purple stars the random walk Laplacian, black dots are the algorithm of [START_REF] Qin | Regularized spectral clustering under the degree-corrected stochastic blockmodel[END_REF]. The parameters of the LFR are n = 50.000, d max = 15 log(n), d = 5, τ 1 = -3, τ 2 = 1, min_class_size = n/8, max_class_size = n/8. Variances are taken over 8 runs of k-means.

High values of M correspond to good quality partitions. Alternatively, the partition quality is evaluated in terms of (normalized) posterior negative log-likelihood of the DC-SBM, L:

L = - 1 2|E|   (ij)∈E log θi θj Ĉˆ i , ˆ j n + (ij) / ∈E log 1 -θi θj Ĉˆ i , ˆ j n   , (30)
where θi = d i / d and Ĉab = i : ˆ i =a j : ˆ j =b A ij / i : ˆ i =a j : ˆ j =b θ i θ j . Good quality clustering correspond, in this case, to low values L.

Table 1 compares the results of different clustering algorithms 10 of 15 real-world networks of increasing size. For all networks, the number of communities, when not available, is estimated through kL (see Subroutine 1) and then the same value is used for all competing techniques (which in general do not provide their own dedicated estimator of k). The underlined numbers in the k column indicate instead that k is known. Furthermore, for all networks, community detection is performed only on the largest connected component of the graph and n, c, Φ, k refer to the characteristics of this dominant connected component.

Given the embedding X, 5 iterations of k-means are run and the best partition ˆ (in terms of k-means is kept) and the scores M, L are computed. To keep stochastic of k-means into account, this step is repeated for 8 times and Tables 1 reports the results of M, L in terms of mean ± standard deviation.

10. The choice of the spectral algorithm considered for comparison is based on two criterion: i) H √ cΦ , B and L sym τ are the state-of-the-art spectral methods for community detection in sparse graphs, while A and L rw are algorithms of great relevance in the literature; ii) these are the methods that are well explained by our unified framework. [START_REF] Leskovec | SNAP Datasets: Stanford large network dataset collection[END_REF][START_REF] Wayne | An information flow model for conflict and fission in small groups[END_REF][START_REF] Lusseau | The bottlenose dolphin community of doubtful sound features a large proportion of long-lasting associations[END_REF][START_REF] Girvan | Community structure in social and biological networks[END_REF][START_REF] Lada | The political blogosphere and the 2004 us election: divided they blog[END_REF]. The value of k is estimated according to Algorithm 2 only if it is not known. When k is known in advance, it appears in the table as an underlined number. All algorithms are performed on the largest connected component of the network. For each data set, the first row in the score L Equation (30), while the second is M, Equation ( 29). When the standard deviation is not indicated it means that it is below 0.01. The first score column of Table 1 indicates the output of the meta-algorithm Algorithm 1 for which ζ p is estimated from ρ(B)/s ↓ R,p (B). The second score column provides the output of Algorithm 2 in which the last step of projection on the hypersphere (Step 4 of Algorithm 0) is not performed. Comparing these two columns, it is clear that Algorithm 2 generally provides much better partitions both in terms of M and L, as a consequence of ρ(B)/s ↓ R,p (B) being an inappropriate estimator for ζ p in general. The following columns display the results scores obtained by Algorithm 2 with the projection step (highlighted in cyan), by clustering based on the leading eigenvectors of A, by the Bethe-Hessian as per [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF], by the non-backtracking matrix as per [START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF], by the random walk Laplacian as per [START_REF] Shi | Normalized cuts and image segmentation[END_REF] and by the symmetric normalized Laplacian of [START_REF] Qin | Regularized spectral clustering under the degree-corrected stochastic blockmodel[END_REF].

The algorithm of [START_REF] Shi | Normalized cuts and image segmentation[END_REF] based on L rw provides in certain cases very competitive partitions (e.g., in the Wikipedia data set) but is quite unreliable as it may dramatically fail in others (GNutella P2P and Polblogs). Algorithm 2 without the normalization step provides systematically good partitions, all comparable to those of L sym τ . This is an evidence that Algorithm 2 effectively produces a node embedding which is significantly resilient to degree heterogeneity. Finally, the improved version of Algorithm 2 including the projection step further improves the quality of the partition on most data sets, providing on all data sets but Wikipedia the highest reported modularity and lowest, or second lowest measured value of L.

To test Algorithm 2 we computed the vector ζ up to machine error precision, that is, the convergence stopping criterion in Subroutine 2 is met when r's update is below machine precision. Notably, due to the non-linearity of the k-means step, larger errors in the estimate of the values of ζ lead, within a certain range, to the same partition. The same classification precision can therefore be reached in fewer iterations, thus faster, if needed.

To further test the efficiency of Algorithm 2, we compared the clustering performances on some large size Lancichinetti-Fortunato-Radicchi networks [START_REF] Lancichinetti | Benchmark graphs for testing community detection algorithms[END_REF] that notoriously constitute a benchmark to evaluate clustering performance on synthetic graphs with ground truth labels that resemble real-world networks. In this case we have a clear definition of communities (provided by the generative model), hence the performance is evaluated in terms of overlap. Figure 8 shows the performance of the different spectral algorithms considered earlier as a function of the mixing parameters 1 -µ, the fraction of links connecting nodes in the same community. The result is in complete agreement with the observations in Table 1 andFigure 5 which all evidence: i) the importance of choosing a proper parametrization for H r ; ii) that Algorithm 2 is robust on real-world networks.

Let us mention that, comparing the columns "Alg 2 wp " and "H √ cΦ ", it appears evident that Algorithm 2 indeed achieves good quality embeddings on real-world graphs. In addition, it appeared from our simulations that the projection step on the unitary hypersphere of the rows of X significantly helps (both in terms of quality of the partitioning and of the variance of the score) the convergence of k-means at a low computational cost.

As a side comment, although formally not strictly comparable on even grounds, we evaluated the performances reached by Algorithm 2 against the popular greedy Louvain method [START_REF] Vincent D Blondel | Fast unfolding of communities in large networks[END_REF] using its scikit-network implementation [START_REF] Pedregosa | Scikit-learn: Machine learning in Python[END_REF]. The Louvain method comes with an estimate of the number of communities and relies on a different notion of communities than the one we used. Specifically, it is a hierarchical algorithm which looks for a partition maximizing the modularity, while Algorithm 2 relies on the DC-SBM assumption. Figure 9 shows that Algorithm 2 largely outperforms the Louvain method on a DC-SBM graphs, both in terms of overlap (left) and modularity (right), except below α c . There, the Louvain method reaches higher modularity, but this is likely an incidental artifact of the modularity optimization constraint [START_REF] Guimera | Modularity from fluctuations in random graphs and complex networks[END_REF].

The situation of real networks is less straightforward. We indeed observed that when the number of communities estimated by both algorithms is similar, they both produce node partitions with a similar modularity. This is however no longer the case when the estimated k are more distinct (e.g., on GNutella P2P). While we recommend our estimation method of k due to its interpretability (that we discussed all along the article), we acknowledge that the Louvain algorithm often provides very competitive outcomes, at a smaller computational cost but at the expense of any theoretical guarantee. Similar considerations can be made in relation to the more recently proposed Leiden method [START_REF] Vincent A Traag | From louvain to leiden: guaranteeing well-connected communities[END_REF].

Comments and Conclusions

This article provides a critical analysis of the current state-of-the-art spectral algorithms for community detection in sparse graphs. While these methods were mostly developed in parallel and from differing lines of arguments, we showed that they all relate and can be jointly optimized to tackle community detection in heterogeneous and sparse graphs. This in passing bridges ideas from statistics and physics. Notably, we discovered that the hyperparameters involved in these algorithms, so far chosen statistically, must instead be smartly selected, this selection being a function of the difficulty of the clustering problem. Besides this fundamental observation, we also showed that there is an efficient (with computational complexity linear in n) routine capable of estimating these optimal parameters, which comes from a deep relation between the non-backtracking matrix and the Bayes optimal solution.

The necessity in all these methods to adapt the choice of the graph representation matrix to the hardness of the problem at hand is a profound observation, likely not restricted to community detection on sparse graphs. Indeed, the classical spectral algorithms based on Laplacian matrices may all be interpreted as a continuous relaxation of an NP-hard discrete optimization problem, providing good (and sometimes provably optimal) solutions on dense but generally not on sparse networks or data. In light of this observation, it may appear that some common optimization problems, beyond clustering, have similarly been devised in the past to perform well in easy scenarios, but would require a task-hardness related update to perform better in harder settings.

Practically, examples of how the results of this article can be generalized to more involved scenarios include the dynamical community detection studied in [START_REF] Dall | Community detection in sparse time-evolving graphs with a dynamical bethe-hessian[END_REF] on time-evolving graphs and the i.i.d. sparsification (improving cost efficiency) of kernel matrices for spectral clustering [START_REF] Dall | Nishimori meets bethe: a spectral method for node classification in sparse weighted graphs[END_REF]. The latter work shows that the common kernel spectral clustering algorithms can be interpreted as a dense limit of a weighted Bethe-Hessian, that can be conveniently adopted in both the sparse and dense regimes. Beyond spectral clustering, we envision that a carefully parametrized Bethe-Hessian matrix can be exploited as a powerful alternative to the combinatorial graph Laplacian.

• we generate n nodes and, for each node, draw independently the class p to which it belongs with probability π p . Also, for each node i, we draw independently its intrinsic connectivity variable θ i ∈ [θ min , θ max ], according to an underlying distribution ν(θ) verifying E(θ) = 1 and E(θ 2 ) = Φ.

• For each pair of nodes (i, j), we generate an edge with probability min θ i θ j C i , j

n , 1 where i is the class of node i.

Theorem 1 (Percolation threshold in the DC-SBM) Consider a graph G generated according to the above DC-SBM procedure. Suppose that the constant vector is the Perron-Frobenius eigenvector of CΠ associated to eigenvalue c > 0. Then, for all large n with high probability, the graph G has a giant component if and only if cΦ > 1.

Proof In our proof we will use (Bollobás et al., 2007, Theorem 3.1) which requires to write the graph G as a graphon. Let us define the "node" variable x = θv ∈ S, with probability distribution µ(x), where • θ ∈ Θ = [θ min , θ max ] is a scalar that encodes the degree heterogeneity, distributed according to ν(•) and satisfying

Θ dν(θ) = 1 (normalization), E[θ] = Θ θ dν(θ) = 1 and E[θ 2 ] = Θ θ 2 dν(θ) = Φ.
• v ∈ {e 1 , . . . , e k } with e p ∈ R k satisfying e p, = δ p ; v is a k-dimensional random vector with law P(v = e p ) = π p ; that is, v encodes the class to which node x belongs.

• The probability density of x is denoted by µ(x) and is defined on S = ∪ k =1 S where S = [θ min e , θ max e ]. The measure µ(x) is equal to π ν(θ) if x ∈ S and it is indeed normalized:

S dµ(x) = k =1 π Θ dν(θ) = k =1 π = 1.
Given two such nodes, we next define the kernel κ(x, y) = x T Cy. The kernel κ satisfies the three conditions of (Bollobás et al., 2007, Definition 2.7) and is a so-called graphical kernel.

With these notations at hand, the generative procedure of the DC-SBM under consideration is equivalent to: drawing n independent values {x i } i=1,...,n from µ(•) and then generating the edges independently according to p ij = min κ(x i ,x j ) n

, 1 ; we thus fall precisely under Bollobas' framework.

In order to use the core argument of (Bollobás et al., 2007, Theorem 3.1), we still need to define the linear operator T κ , an operator on f : R k → R, endowed with the norm f 2 2 = S f 2 (x) dµ(x), as:

∀ x ∈ S, (T k f )(x) = S κ(x, y)f (y) dµ(y)
The square 2-norm of this operator, T κ 2 , reads

T κ 2 2 = sup f 2 ≤1 S (T κ f ) 2 (x) dµ(x). (31)
According to (Bollobás et al., 2007, Theorem 3.1), with high probability, a giant component appears in G if and only if T κ > 1. We thus are left to evaluating T κ for the kernel x T Cy: we will show that it equals cΦ, first by showing that T κ ≤ cΦ and then finding a particular function f for which the bound is attained.

Let f be a function from R k to R. We introduce its associated vector ω f ∈ R k as follows.

(T κ f )(x) = S κ(x, y)f (y) dµ(y) = S x T Cyf (y) dµ(y) = x T C k =1 π e Θ θ f (θe ) dν(θ) ≡ x T CΠω f , (32)
where ∀ , ω f, = Θ θ f (θe ) dν(θ). From Equation (32), one thus has:

S (T κ f ) 2 (x)dµ(x) = S ω T f ΠCxx T CΠω f dµ(x) = Φ ω T f ΠC k =1 π e e T CΠω f = Φ ω T f ΠCΠCΠω f . (33)
Injecting the result of Equation (33) into Equation ( 31) we obtain

T κ 2 2 = Φ sup f 2 ≤1 ω T f ΠCΠCΠω f = Φ sup f 2 ≤1 ω T f ΠCΠCΠω f ω T f Πω f • ω T f Πω f ≤ Φ sup f 2 ≤1 ω T f ΠCΠCΠω f ω T f Πω f • sup f 2 ≤1 ω T f Πω f (34)
Analyzing the first element, we can write for v f = Π 1/2 ω f :

sup f 2 ≤1 ω T f ΠCΠCΠω f ω T f Πω f = sup f 2 ≤1 v T f Π 1/2 CΠCΠ 1/2 v f v T f v f ≤ sup v∈R k v T Π 1/2 CΠCΠ 1/2 v v T v i.e.: sup f 2 ≤1 ω T f ΠCΠCΠω f ω T f Πω f ≤ s ↓ 1 (Π 1/2 CΠCΠ 1/2 ) = s ↓ 1 (CΠ) 2 = c 2 ,
as, by hypothesis, c is the Perron eigenvalue of CΠ, and as such is larger than the modulus of any other eigenvalue: s ↓ 1 (CΠ) 2 is indeed c 2 . We can therefore rewrite Equation (34) as

T κ 2 2 ≤ c 2 Φ sup f 2 ≤1 ω T f Πω f .
Analyzing the right-hand side term, we have

ω T f Πω f = k =1 π ω 2 f, = k =1 π Θ θf (θe ) dν(θ)
2 which, by Cauchy-Schwartz's inequality is bounded as

Θ θf (θe ) dν(θ) 2 ≤ Θ θ 2 dν(θ) Θ f 2 (θe ) dν(θ) ≤ Φ,
where in the last step we used the fact that the norm of f is less than or equal to one. We thus obtain that sup

f 2 ≤1 ω T f Πω f ≤ Φ and conclude that T κ 2 2 ≤ (cΦ) 2 . ( 35
)
We are thus left to showing that there exists a function f for which the bound is reached. Let us consider f (x) = x / √ Φ. It is easy to check that this function has unit norm

f 2 2 = S x √ Φ 2 dµ(x) = 1 Φ S θ 2 e 2 dµ(x) = 1 Φ S θ 2 dµ(x) = 1. (36)
Furthermore, observe that ω f = √ Φ1 k . Then, we have

T κ 2 2 = Φ sup f 2 ≤1 ω T f ΠCΠCΠω f ≥ Φω T f ΠCΠCΠω f = Φ 2 1 T k Π(CΠ) 2 1 k = (cΦ) 2 (37)
as, by hypothesis, we suppose that CΠ1 k = c1 k . Combining Equations ( 35) and (37), we find that T κ = cΦ so that, from (Bollobás et al., 2007, Theorem 3.1), we conclude that the percolation transition is at cΦ = 1 and a giant component exists if and only if cΦ > 1.

This appendix is structured as follows. First we enunciate two lemmas and one corollary to state that when an arbitrary connected graph has at most one loop, the spectrum of its associated non-backtracking matrix B does not have real eigenvalues that are larger than one in modulus.

We then proceed by arguing that with high probability, for large DC-SBM graphs, the small connected components do not have more than one loop, implying that all the real eigenvalues of B that are larger than one in modulus come from the giant component. Finally we show that this last statement, implies that the two definitions of ζ (j) p given in Equations (38,39) are indeed equivalent with high probability.

Lemma 1 (Spectrum of B on a tree) Let T (V, E) be a tree. Then, all the eigenvalues of its associated non-backtracking matrix B(T ), as defined in Equation (11), are equal to zero.

Lemma 2 (Spectrum of G plus a node) Let G be a graph obtained by adding one node and one edge to the graph G, i.e. for i / ∈ V and an arbitrary j ∈ V, V = V ∪ {i} and E = E ∪ (ij). Then, all the non-zero eigenvalues of B(G ) are equal to the non-zero eigenvalues of B(G), where B(G) is the non-backtracking matrix defined in Equation (11) of the graph G.

Thus concluding that any eigenvector g of B is associated to eigenvalue zero. Note that in other words, this means that B is nilpotent.

Proof [Lemma 2] Let i be the newly added node and j the node in V to which i is attached. The matrix B(G ) can be written by adding to the matrix B(G) two rows and two columns corresponding to the directed edges (ij) and (ji). We introduce the notation 1 (•j) ∈ R 2|E| . Its element-wise definition reads for all (yx) ∈ E d (the set of directed edges of G) as 1 Proof [Corollary 1] A connected graph with |E| = |V| can be obtained by adding an edge between any two nodes of a particular tree defined on V. The graph G thus contains only one loop. We now apply Lemma 2 "backwards", i.e., removing leaves from G without affecting the non-zero eigenvalues of B.

By iteratively removing all leaves, the graph G reduces to a loop. We are thus left to prove that ρ(B) = 1 on a loop. It is straightforward to check that on any d-regular graph (all nodes having d neighbours), the vector 1 2|E| is the Perron-Frobenius eigenvector of B with eigenvalue equal to d-1. A loop is a d-regular graph with d = 2, hence the result.

With the results of the former two Lemmas, we proceed to argue that all the eigenvalues of B that are larger than one in modulus come from the giant component. Since the graph is disconnected, S(B) = n CC j=1 S B (j) , i.e., each connected component contributes independently to the eigenvalues of B. The expected size of the small connected components grows as 11 log(n). We claim that this implies that the probability that a small connected component contains two or more loops is o n (1), thus tending to zero as n grows to infinity. Applying Lemmas 1 (for zero loops) and Corollary 1 (for one loop), we conclude that, with high probability all the real eigenvalues of B larger than one in modulus will come from the 11. Indeed, the kernel κ defined in the proof of Theorem 1 is irreducible in the sense of (Bollobás et al., 2007, Def. 2.10 ), as the nodes cannot be split into two parts that have no chance of being connected (all the entries of C and all πi are supposed strictly positive). One can thus apply (Bollobás et al., 2007, Thm. 3.12) stating that the small components are with high probability of order log n.

giant component. From the Ihara-Bass formula [START_REF] Terras | Zeta functions of graphs: a stroll through the garden[END_REF], for all the eigenvalues γ = ±1, γ ∈ S(B) ⇐⇒ det[H γ ] = 0. From the definition of ζ p≥2 > 1, then with high probability j is the giant component. We now conclude showing that this statement implies the equivalence between Equations ( 38) and (39).

Thanks to the properties we discussed in Section 2, for p > 1, then certainly ζ (j) p > 1, if it exists. We denote with j the giant component and we want to prove that for r, p > 1

s ↑ p (H (j) r ) = s ↑ p (H r ) (44)
Consider

ζ (j)
p > 1 such that s ↑ p (H ) = 0, then certainly there exists 12 q ≥ p such that s ↑ q (H ζ (j) p ) = 0 and so it is enough to show that p = q. We proceed with a proof by contradiction.

Suppose that p = q, then there exists j = j such that s ↑ 1 H p > 1 such that s ↑ 1 (H (j ) r ) = 0. From the Ihara-Bass formula, r is thus in S(B), so there is a real eigenvalue of B larger than one not coming from the giant component. This is in contradiction with what stated above, thus concluding our argument.

Appendix C. The Eigenvalues of T

We here want to show that s ↓ p (T ) = . The main result of Lemma 3 is to state that this eigenvalue is the p-th largest. Note that the condition s ↑ p (D -rA) < 0 could be loosened, but we don't need a stronger result in the following.

Lemma 4 Let D and A be the degree and adjacency matrices of any arbitrary graph. Let p be an integer ≥ 2. Suppose that there exists r p > 1 such that i/ ∀ r ≥ r p the eigenvalue s ↑ p (D -rA) is simple, and ii/ s ↑ p (D -r p A) < 0. Then:

∀ r ≥ r p , ∂ r s ↑ p (D -rA) < 0.

In words, once the function s ↑ p (D -rA) becomes negative, it is strictly decreasing.

The eigenvalue simplicity assumption of this lemma is technical and enables us to properly define and manipulate the derivative of the p-th smallest eigenvalue. In case of multiplicity, the tools involved are more complicated (see for instance [START_REF] Greenbaum | First-order perturbation theory for eigenvalues and eigenvectors[END_REF] and not included here. Note that according to Claim 1 this assumption is verified for a network generated from the DC-SBM because on a large random graph with independent entries all eigenvalues are simple with probability one.

Lemma 5 Let G(V, E) be a DC-SBM graph generated according to Equation (1), and D and A its degree and adjacency matrices respectively. Consider 2 ≤ p ≤ k and let ζ p be defined as per Equation (5). Let r ≥ ζ p . If s ↑ p (D -rA) < 0 is isolated with associated eigenvector x p , then x p is an eigenvector of L rw -s ↑ p (D-rA)

, whose corresponding eigenvalue is also isolated.

Otherwise, for 1 ≤ p ≤ n and r > 1, if s ↑ p (D -rA) belongs to the bulk of non informative eigenvalues of D -rA, then the corresponding eigenvalue of L rw .

Based on Lemmas 3-5 and Claim 1 we now proceed to the justification of Claim 2, that we recall here for convenience.

Claim 2 Consider the graph G(V, E), built on a sparse DC-SBM as per Equation (1) with k communities. Let ζ p for 2 ≤ p ≤ k be defined as per Equation (5) (imposing ζ 1 = 1) and τ ∈ R be such that ζ 2 p -1 ≤ τ ≤ cΦ -1. Then, under Assumption 1, for all large n with high probability, the p largest eigenvalues of the matrix L rw τ are isolated. In particular,

s ↓ p (L rw ζ 2 p -1 ) = 1 ζ p .
For simplicity, will use the notation s and it is isolated.

Figure 1 :

 1 Figure 1: For all the four plots n = 5000, k = 2, θ = 1 n if not otherwise specified. (Left) spectrum of the matrix L sym in the: A) dense regime (c in /n = 0.06, c out /n = 0.02); B) sparse regime (c in = 6, c out = 2). (Right) first eigenvectors x 1 vs second eigenvector x 2 of L sym in the dense regime for: C) θ = 1 n ; D) θ i ∼ [U(3, 10)] 3 .

  p≥2 exists: it is the giant component. In the following, abusing notation, we simply write ζ p instead of ζ (j) p to refer to the ζ's associated to this giant component. One has ζ 1 = 1 and, if it exists, ζ p≥2 verifies 1

Figure 2 :

 2 Figure 2: (Left): spectrum of the matrix B on the complex plane. The green dots represent the uninformative eigenvalues, the blue diamonds the real eigenvalues outside the bulk and the red squares the real eigenvalues inside the bulk. The lines indicate the theoretical values. For this network n = 5000, k = 3, θ i ∼ [U(3, 10)] 3 , c = 8. The off-diagonal elements of C are distributed according to C p>q ∼ N (c out , c out /k) with c out = 4 and the diagonal elements are then fixed to have CΠ1 k = c1 k . The element of π are distributed as π i ∼ N (1/k, 1/2k). The entries are then rescaled sot that Tr(Π) = 1. (Right): position of the first four real eigenvalues of B as a function of α. For this simulation n = 5000, k = 2, θ i = [U(3, 10)] 3 , c in : 12 → 7, c out : 0 → 5 with c = 6.

Figure 4 :

 4 Figure 4: Behavior of s ↑ p (H r ) as a function of r for 1 ≤ p ≤ k in solid lines, and for p = k + 1 in dotted line. The vertical lines are the theoretical values of c/ν p (dashed dotted line), of ρ(B) (dotted line) and of ν p Φ (dashed lines). For this simulation, n = 50.000, c = 6, c out = 2, k = 4, C p>q ∼ N (c out , 2c out /k) and π ∝ 1 k .

Algorithm 1 nc i d 2 i 3 :

 123 Community Detection with the Bethe-Hessian 1: Input : adjacency matrix of undirected graph G 2: Estimate c = 1 n i d i and Φ = 1 Detect number of classes k according to Remark 1 4: for 2 ≤ p ≤ k do 5:

Figure 5 :

 5 Figure 5: Overlap comparison as a function of the hardness α of the detection problem.Blue sharp diamonds are the result given by Algorithm 1, red diamonds the algorithm using the Bethe-Hessian of[START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF], green squares are for the non-backtracking of[START_REF] Krzakala | Spectral redemption in clustering sparse networks[END_REF], yellow pentagons the adjacency matrix, purple stars the random walk Laplacian, black dots are the algorithm of[START_REF] Qin | Regularized spectral clustering under the degree-corrected stochastic blockmodel[END_REF]. For both graphs n = 50.000,θ i ∼ [U(3, 10)] 4 , π ∝ 1 k , c = 5. (Left): k = 2, c out : 0.05 → 4. (Right): k = 5, π p ∼ N (k -1 , (4k) -1 ), C a>b ∼ N (c out , c out k -1 ), c out : 0.05 → 4. α := 2(c -c out )/ √ c and α c = 2/ √ Φ. Averages over 10 samples.

  Figure 6: (Left) Spectrum of ζ 4 L rw ζ 2 4 -1 , n = 5000, with c in = 4, c out = 0.5, θ i ∼ [U(3, 5)] 4 with C generated according Equation (26). (Right) Realization of G from the SBM with C defined as per Equation (26), with c in = 5 and c out = 1, n = 1500. The position of the nodes is assigned according to the true labels as N (µ i , I 2 ) and the colors are the output of Algorithm 1 on G.

Figure 7 :

 7 Figure 7: (Top: DC-SBM setting) Exact value of c/ν p (p = 2 on the left and p = 3 on the right) in blue, compared to the estimate of ζ p obtained using the dominant eigenvalues of B (green) and its direct computation using Subroutine 2 (red), for problems of different hardness. For this simulation, n = 50.000, k = 3, θ i ∼ [U(3, 10)] 3 , c = 10. For each point π i ∼ N (1/k, 1/2k) and the off-diagonal elements of C are distributed as ∼ N (c out , c out /k) and c out : 0.1 → 7.5. (Bottom: real networks) Estimate of the values of ζ p as function of p, computed as ρ(B)/s ↓ p (B) (green squares) vs the direct computation of ζ p as the output of Subroutine 2 (red diamonds) on the three real networks taken from[START_REF] Leskovec | SNAP Datasets: Stanford large network dataset collection[END_REF]. The value of k is estimated according to Subroutine 1.

Figure 9 :

 9 Figure 9: Comparison of the overlap (left) and modularity (right) resulting from the label estimation of Algorithm 2 (blue diamonds) and Louvain algorithm (orange crosses) applied on the giant component of G, as function of the detection hardness. For this simulation, n = 50.000, k = 2, π ∝ 1 2 , c = 8, c out = 0.05 → 7. Averages are taken over 10 samples.

Corollary 1 (

 1 Spectral radius of a connected graph with |E| = |V|) Let G(V, E) be a connected graph with |E| = |V|. Then ρ(B) = 1, where B is defined in Equation (11) and ρ(•) is the spectral radius. Proof [Lemma 1] Consider an arbitrary vector g ∈ R 2|E| and define g (m) ∈ R 2|E| as g notation d(ij, k ) = m indicates that there exists a non-backtracking path from the edge (ij) to the edge (k ) of length m. From a straightforward calculation one obtains (Bg (m) ) ij = ∈∂j\i (kq):d(j ,kq)=mg kq = g , for any two edges we have that d(ij, k ) ≤ n-1, so there is a value of m c ≤ n-1, which represents the maximal distance between any two directed edges, such that, for all vectors g,Bg (mc) = 0.(42)This relation comes from the fact that, by definition of m c , no two edges are at a distance equal to m c + 1, so {(k : d(ij, k ) = m c + 1} = ∅ for any edge (ij). Now, let us consider g to be an eigenvector of B, such that Bg = γg. Then we can write 0 = Bg (mc) = B 2 g (mc-1) = • • • = B mc-1 g = γ mc-1 g.

(

  •j) yx = δ xj . Similarly, we define 1 (j•) ∈ R 2|E| with 1(j•) yx = δ yj . Denote with M = (1 (j•) , 0 2|E| ) ∈ {0, 1} 2|E|×2 and M = (0 2|E| , 1 (•j) ) ∈ {0, 1} 2|E|×2, where 0 2|E| is the null vector of size 2|E|. The matrix B(G ) can be written as:B(G ) = B(G) M M T 0 2×2We now look for the non-zero eigenvalues of γ of B(G ). Recall that I n is the n × n identity matrix. Consider γ ∈ C * . Using a block matrix determinantal formula one has:det B(G ) -γI 2|E | = det(-γI 2 )det B(G) -γI 2|E| + 1 γ M M T . It is straightforward to check that M M T = 0 2|E|×2|E| , thus we have det B(G ) -γI 2|E | = γ 2 det B(G) -γI 2|E| .A non-zero eigenvalue of B(G ) (cancelling the determinant) is thus necessarily also an eigenvalue of B(G), ending the proof.

  are in the spectrum of B. We may conclude that, if ∃j such that ζ (j)

  it is easy to show that for r > d max -1, the matrix H r is positive definite. Consequently, there exists r ≥ ζ (j)

c

  for 1 ≤ p < k. Consider the following equivalent identities for the matrixT = ΠC c -Π1 k 1 T k : ΠC c -Π1 k 1 T k a p = s ↓ p (T )a p CΠ c -1 k 1 T k Π b p = s ↓ p (T )b p ; b p = Π -1 a p Since CΠ1 = c1 k and 1 T k Π1 k = 1, then b k = 1 k . By introducing d p = Π 1/2 b p , we can write: Π 1/2 CΠ 1/2 c -Π 1/2 b k b T k Π 1/2 d p = s ↓ p (T )d p Π 1/2 CΠ 1/2 c -d k d T k d p = s ↓ p (T )d p 12. As S(Hr) = ∪ n CC j=1 S H (j) r .The fact that the matrix L rw -s ↑ p (D-rA)has an eigenvalue equal to 1/r comes easily by construction of L rw -s ↑ p (D-rA)

  ↑ p (D -rA) = s p (r). The first part of Claim 2 asserts that for ζ 2 p -1 ≤ τ ≤ cΦ -1 the p largest eigenvalues of the matrix L rw τ are isolated. According to Claim 1, for ζ p ≤ r ≤ √ cΦ the eigenvalue s p (r) is isolated. Thanks to Lemma 4, since s p (ζ p ) < 0 we have: ∂ r s p (r) < 0 for all r in the interval ζ p ≤ r ≤ √ cΦ, implying s ↑ p (D -rA) < 0 for all r in the interval ζ p ≤ r ≤ √ cΦ. The hypotheses of Lemma 3 and Lemma 5 are therefore verified. We can then assert that the eigenvalue 1/r is the p-th largest of the matrix L rw -s ↑ p (D-rA)

Figure 10 :

 10 Figure 10: Plot of the f rt (r ) (in purple). The blue continuous line indicates the value of r t , the red dashed dotted line the theoretical value of ζ p = c/ν p , the two dashed black lines are the two roots of f rt (r ), the smaller of which is r t+1 . For this simulation n = 50 000, k = 2, c = 5, c out = 2.5, r = √ cΦ, θ i ∼ [U(3, 10)] 3 , π ∝ I k .

  that its smallest eigenvalue is null and that s ↑ p≥2 H

				(j) 1	> 0.
	Consequently, ζ	(j) 1 = 1 always exists and if ζ	(j) p≥2 exists, it is strictly superior to 1.
	2. If ζ p≥2 exists, then ζ (j) 2 , . . . , ζ (j)	(j) p-1 necessarily exist and are smaller or equal to ζ	(j) p . In
	fact, if ζ	(j) p exists, it means that at r = ζ	(j) p , 0 is the p-th smallest eigenvalue of H r : (j)
	the p -1 smallest are thus ≤ 0. Given that these p -1 smallest are ≥ 0 at r = 1 and
	by continuity of s ↑ p (H r ), they necessarily cross zero before ζ (j)	(j) p . Similarly, if ζ p≥2 does (j)
	not exist, then all ζ	(j)

  . . , ξ p,n grows linearly with d i (and thus becomes essentially negligible). Now, proceeding as in Equation (19), we can compute Var \i [(H ζp u p ) i |A, i ] from a direct calculation of E \i [(H ζp u p ) 2 i |A, i ] and E 2 \i [(H ζp u p ) i |A, i ], and we obtain Var \i [(H ζp u p

Table 1 :

 1 Modularity and log-likelihood comparison on real networks

	Vip	HepPh	Wikipedia	GNutella P2P	Politicians	Power grid	GrQc	Facebook	Tv	Polblogs	Mail	Football	Polbooks	Dolphins	Karate	Data set
	11565 11.6 4.4 53	11204 21.0 6.2 60	7066 28.3 5.1 22	6299 6.6 2.7 4	5908 14.1 3 62	4941 2.7 1.5 25	4158 6.5 2.8 29	4039 43.7 2.4 55	3892 8.9 3 41	1222 27,4 3 2	1133 9.6 1.9 21	115 10.7 1 12	105 8.4 1.4 3	62 5 1.3 2	34 4.6 1.7 2	n c Φ k
	0.34 ±0.01	2.60 ±0.03	0.39 ±0.02	1.75 ±0.06	0.20	1.93	0.20 ±0.01	3.32 ±0.01	0.52 ±0.02	1.78 ±0.04	0.32 ±0.01	3.63 ±0.02	0.50 ±0.02	2.40 ±0.04	0.52 ±0.01	1.17 ±0.04	0.55 ±0.02	2.01 ±0.03	0.43	1.50	0.45 ±0.01	1.95 ±0.01	0.60	0.83	0.51	1.17	0.38	1.25	0.37	0.86	Alg 1B
	0.58 ±0.15	2.48 ±0.13	0.37 ±0.06	1.98 ±0.05	0.23 ±0.01	1.88 ±0.01	0.26	3.28	0.83 ±0.01	1.60 ±0.05	0.92	2.61 ±0.01	0.68 ±0.09	2.22 ±0.11	0.76 ±0.02	1.05 ±0.06	0.70 ±0.06	1.89 ±0.09	0.43	1.50	0.49 ±0.01	1.89 ±0.02	0.60	0.83	0.50	1.16	0.38	1.25	0.37	0.86	Alg 2wp
	0.60 ±0.01	2.22 ±0.01	0.57	1.52 ±0.01	0.26	1.88	0.40	3.21	0.84 ±0.01	1.43 ±0.01	0.93	2.57 ±0.01	0.80	1.89 ±0.01	0.76	0.88 ±0.01	0.84	1.50 ±0.01	0.43	1.50	0.52 ±0.01	1.88 ±0.01	0.60	0.83	0.51	1.17	0.38	1.25	0.37	0.86	Alg 2
	0.27 ±0.01	2.72 ±0.02	0.42 ±0.04	1.74 ±0.08	0.15	1.99±0.01	0.19±0.01	3.33	0.54 ±0.02	1.94 ±0.06	0.17 ±0.01	3.84 ±0.01	0.49 ±0.04	2.39 ±0.06	0.44 ±0.03	1.13 ±0.03	0.53 ±0.04	1.99 ±0.09	0.25	1.62	0.31	2.13 ±0.01	0.60	0.83	0.46	1.26 ±0.01	0.22 ±0.01	1.38 ±0.01	0.37	0.86	A
	0.32 ±0.01	2.64 ±0.03	0.39 ±0.02	1.79 ±0.04	0.17	1.97	0.14	3.37	0.50 ±0.02	1.90 ±0.02	0.35 ±0.03	3.63 ±0.03	0.50 ±0.02	2.40 ±0.04	0.47 ±0.01	1.21 ±0.03	0.53 ±0.02	2.02 ±0.05	0.27	1.63	0.40	2.03 ±0.01	0.60	0.83	0.50	1.21	0.32 ±0.03	1.30 ±0.03	0.37	0.86	H √ cΦ
	0.30 ±0.01	2.67 ±0.03	0.44 ±0.03	1.68 ±0.05	0.14±0.01	1.96	0.13±0.01	3.38±0.01	0.47 ±0.01	1.92 ±0.02	0.28 ±0.01	3.71 ±0.01	0.50 ±0.01	2.39 ±0.01	0.47 ±0.02	1.21 ±0.03	0.53 ±0.02	2.05 ±0.05	0.24	1.65	0.36 ±0.01	2.05 ±0.01	0.60	0.83	0.47 ±0.01	1.28 ±0.01	0.22	1.38	0.37	0.86	B
	0.50 ±0.12	2.65 ±0.01	0.25 ±0.09	2.12 ±0.07	0.38 ±0.02	1.94 ±0.03	0.00	3.41	0.78 ±0.07	1.77 ±0.11	0.92	2.61 ±0.01	0.58 ±0.10	2.44 ±0.09	0.75 ±0.05	1.07 ±0.14	0.59 ±0.12	2.06 ±0.15	0.00	1.73	0.50 ±0.02	1.99 ±0.03	0.60	0.83	0.50	1.16	0.38	1.25	0.36	0.98	L rw
	0.54 ±0.01	2.18 ±0.01	0.51 ±0.01	1.44 ±0.01	0.27	1.88	0.35	3.27	0.73 ±0.01	1.40 ±0.01	0.86 ±0.01	2.80 ±0.02	0.78	1.90 ±0.02	0.59 ±0.01	0.87 ±0.01	0.79 ±0.01	1.47 ±0.01	0.43	1.50	0.51	1.87 ±0.01	0.60	0.83	0.51	1.22	0.38	1.25	0.37	0.86	L sym τ

in this case of two equal-sized classes

Recall that νp is the p-th largest eigenvalue of CΠ.

Note that s ↑ 1 (T ) = 0.

Acknowledgments

Couillet's work is supported by the IDEX GSTATS DataScienceChair and the MIAI Large-DATA Chair at University GrenobleAlpes. Tremblay's work is partly supported by the CNRS PEPS I3A (Project RW4SPEC) as well as the French National Research Agency in the frameworks of the ANR GraVa (ANR-18-CE40-0005), the "Investissements d'avenir" program (ANR-15-IDEX-02) and the LabEx PERSYVAL (ANR-11-LABX-0025-01). The authors thank Lennart Gulikers and Simon Barthelmé for fruitful discussions and the anonymous reviewers for the constructive comments that helped to improve the article.

Appendix A. Existence of a Giant Component in the Sparse DC-SBM

In this appendix we state and prove Theorem 1 which provides a necessary and sufficient condition for the existence of a giant component in a multi-class (with unequal sizes) sparse DC-SBM. This theorem is an application of (Bollobás et al., 2007, Theorem 3.1) to the present DC-SBM setting. Note in particular that the number of nodes in each class p is deterministic and equal to nπ p . In Bollobas' setting, however, the class of each node i is drawn independently and is equal to p with probability π p : nπ p is thus only the expected number of nodes in class p. Asymptotically, we will see that this of course does not affect the result. To be precise, the generative procedure under analysis in the theorem is the following:

• letting k be the number of classes, and n the number of nodes, we consider a symmetric matrix C ∈ R k×k where C > 0 encodes the affinity between classes and ; a diagonal matrix Π = diag(π 1 , . . . , π k ), where 0 < π < 1 is the expected proportion of nodes in class and Tr(Π) = 1.

Appendix B. Definition of ζ on the Sparse DC-SBM

In this Appendix we show that for a graph generated from the sparse DC-SBM model, defined in Equation ( 1), the two following definitions of ζ

p (where j is the giant component) are indeed equivalent with high probability for p ≥ 2 if ζ (j) p exists:

We will further discuss that, with high probability, ζ

> 1 exists only for the giant component.

Given that the d k are eigenvectors of a symmetric matrix, then d T p d q = δ pq , so, for p < k

From this last equation the result follows directly.

Appendix D. The Excited States of the Ising Hamiltonian on G

In this section we give an explicit connection between the problem described in Section 2.3.2 and its statistical physics analogue. Consider the graph G obtained from the k-class DC-SBM model as per Equation (1). To study the stability of the configurations of σ as a function of r, one needs to compute the free energy F , defined in Equation ( 46) (which we recall is a function of the moments of the Boltzmann distribution µ(σ)):

To evaluate F , a common way to proceed is to exploit the Bethe approximation that is exact on a tree and is thus relevant on sparse graphs. By denoting with µ ij (•) and µ i (•) the edge and the node marginals of µ(•), the Bethe free energy is defined as:

By a direct computation, one can verify that

where D KL (• •) is the Kullback-Leibler divergence which is always greater or equal to zero. In order to find the best estimate of the exact free energy, one needs to find the minimum of F Bethe with respect to m, χ. At large r one sees that F Bethe has a unique minimum in m = 0, χ = 1 r 1 2|E| , whereas for r small enough, there appear multiple minima and m = 0 is a saddle point. To quantitatively study this behavior, one needs to consider the Hessian matrix of the Bethe free energy at the paramagnetic point, i.e.:

The eigenvectors corresponding to the negative eigenvalues of the Bethe-Hessian matrix H r represent the directions along which the paramagnetic phase becomes unstable and they will have a non-trivial structure. We now detail an approach inspired from [START_REF] Janusz | Hołyst. Ising model on two connected barabasi-albert networks[END_REF] but adapted to the Bethe approximation to study the stability of the paramagnetic phase. The following equations hold:

It is possible to verify that by plugging the ansatz of Equations (50,51), one recovers the factorization of Equation ( 47). We define the following change of variable ξ ji (σ i ) = e h ji σ i , resulting in

ξ ji (-1) . This leads to the following expressions:

We next linearize Equations ( 52) and ( 53) at the paramagnetic point, i.e., for h → 0 in the high temperature regime (large r). For mathematical convenience we here assume the regime where c 1 and r of the order of c or larger. We will comment in the end on how the results we obtain still hold in the sparse regime in which c = O n (1).

Further developing the second term of the summation j∈V

we obtain

Equation ( 56) gives a self consistent relation between the magnetizations of the stable (or metastable) configurations. We study the magnetizations taking an average over the realizations of A as proposed in [START_REF] Janusz | Hołyst. Ising model on two connected barabasi-albert networks[END_REF])

Defining

We thus obtain the relation controlling the stability of the solutions:

The transition temperatures are hence given by r = ν p Φ, as expected. For r > cΦ Equation (59) only admits the trivial solution s = 0 which represents the paramagnetic configuration. At r = cΦ there is a paramagnetic-ferromagnetic transition, consistent with the result of [START_REF] Leone | Ferromagnetic ordering in graphs with arbitrary degree distribution[END_REF]. For smaller values of r there appear other directions (aligned to the eigenvectors of CΠ) along which the Bethe free energy has local minima.

The main result of this appendix is thus to argue that (i) there are only k directions along which the paramagnetic phase can get unstable, hence H r can have at most k negative eigenvalues, so that s ↑ k+1 (H r ) ≥ 0 for r > 1; (ii) the directions of instability precisely correspond to the eigenvectors of the matrix CΠ.

To conclude, we comment the case c = O n (1). We just detailed an approach to understand and interpret the problem at hand with tools taken from statistical physics. For mathematical convenience we had to assume the average degree c to be rather large and we were able to predict the exact positions of the transition temperatures, which coincide with the largest eigenvalues of the non-backtracking matrix. Our approach in determining the positions of these eigenvalues is non-rigorous and does not intend to substitute the existing proofs, but rather to give a physical intuition of the problem at hand. In particular, the spectrum of B has been rigorously studied in (Gulikers et al., 2017b;[START_REF] Bordenave | Non-backtracking spectrum of random graphs: community detection and non-regular ramanujan graphs[END_REF] in the regime c = O n (1). The conclusion above is not altered: the matrix B has k isolated eigenvalues, so only the k smallest eigenvalues of the Bethe-Hessian matrix can become negative, in correspondence to the emergence of new local minima in the Bethe free energy profile. According to our earlier argument, we expect these minima (represented by the eigenvectors of H r ) to be correlated with the class structure of A.

Appendix E. Supporting Arguments for Claim 2

To support Claim 2, we will rely on the following three intermediary results that will be proved subsequently. We now proceed to the proof of the three lemmas. For the sake of clarity we here enunciate the Courant-Fischer theorem that will be of fundamental use in the sequel.

Theorem 2 (Courant-Fischer, see for instance [START_REF] Bhatia | Matrix analysis[END_REF] Let M ∈ C n×n be a Hermitian matrix and U a vector subspace of C n . Then, Recalling that L rw τ and L sym τ share the same spectrum, one has in particular:

Thus, proving the lemma amounts to proving that s ↑ p ( Mr ) = 0, which we do now by first proving that s ↑ p ( Mr ) ≤ 0 and then that s ↑ p ( Mr ) ≥ 0. Denote by X = (x 1 | . . . |x p ) ∈ R n×p the matrix concatenating the eigenvectors associated to the p smallest eigenvalues of M r . One has:

Now define the vector space V as V = span D 1/2 -sp(r) X . Since s p (r) is strictly negative by hypothesis, D -sp(r) 0 (is positive definite), which in turn implies that dim(V ) = p. By

-sp(r) Xu. From the Courant-Fischer theorem, we can write:

where for the last step we exploited s ↑ q (M r ) ≤ 0 and D -sp(r) 0. We thus conclude that s ↑ p ( Mr ) ≤ 0.

To prove that the equality holds, we exploit the second relation of the Courant-Fischer theorem. We define X = (x p | . . . |x n ) ∈ R n×(n-p+1) the matrix concatenating the eigenvectors associated to the n -p + 1 largest eigenvalues of M r . For q ≥ p, s ↑ q (M r ) ≥ 0. We further define

X , satisfying dim(W ) = n -p + 1. We can write:

As a consequence, s ↑ p ( Mr ) ≥ 0. Combining both inequalities, we obtain that s ↑ p ( Mr ) = 0. The fact that the eigenvectors are shared comes from the following. Let x p (r) be the eigenvector of D -rA associated to s p (r), i.e., [D -rA]x p (r) = s p (r)x p (r). This can be re-written as L rw -sp(r) x p (r) = 1 r x p (r).

Proof [Lemma 4] Once again, let us write s p (r) = s ↑ p (D -rA) to lighten notations. Let p, r p be an integer ≥ 2 and a scalar > 1 (if they exist) such that s p (r p ) < 0. Let r ≥ r p . As s p (r) is an eigenvalue, ∃ x p (r) with ||x p (r)|| 2 = 1 such that:

which implies in particular that x p (r) T Ax p (r) = 1 r x p (r) T D -sp(r) x p (r). As we suppose s p (r) to be simple, we can apply the eigenvalue perturbation theorem (see, for instance [START_REF] Greenbaum | First-order perturbation theory for eigenvalues and eigenvectors[END_REF]):

In Equation ( 60), D -sp(rp) 0 (as s p (r p ) < 0 by hypothesis), consequently ∂ r s p (r)| r=rp < 0.

We now want to show that for all r > r p , ∂ r s p (r) < 0. We proceed with a proof by contradiction.

Suppose that there exists a value r > r p such that ∂ r s p (r)| r=r ≥ 0. From Equation (60) it follows that a necessary (but not sufficient) condition to be verified is that s p (r )≥0. From a continuity argument on the function s p (r), and as s p (r p ) < 0 and s p (r ) ≥ 0 with r > r p , there exists r ∈ (r p , r ) such that ∂ r s p (r)| r=r > 0 and s p (r ) < 0. Invoking once again Equation ( 60), no such r can exist, invalidating the hypothesis we made by absurd. We thus conclude that

finishing the proof.

Proof [Lemma 5] Consider the values {ζ p } 2≤p≤k as defined in Eq.( 5). As G is generated from a DC-SBM we know from Claim 1 that these [START_REF] Saade | Spectral clustering of graphs with the bethe hessian[END_REF]. We thus have s p (ζ p ) < 0, ∀ 2 ≤ p ≤ k + 1 and it is isolated by hypothesis. We can apply Lemma 4, for r p = ζ p : for all r ≥ ζ p , ∂ r s p (r) < 0 and s p (r) < 0. Consequently, there exists a unique value of r ≥ r, satisfying s p+1 (r ) = s p (r). From lemma 3, we know that:

• as s p (r) < 0: x p (r), the eigenvector of D -rA associated to s p (r), is also the eigenvector of L rw -sp(r) associated to its p-th largest eigenvalue, 1/r.

• as s p+1 (r ) < 0: x p+1 (r ), the eigenvector of D -r A associated to s p+1 (r ), is also the eigenvector of L rw -s p+1 (r ) = L rw -sp(r) associated to its (p + 1)-th largest eigenvalue, 1/r . Thus, the p-th (resp. (p + 1)-th) largest eigenvalue of L rw -sp(r) is 1/r (resp. 1/r ). Consider r ≥ ζ p . By hypothesis, s p (r) is an isolated eigenvalue, that is, we can write

for some constant κ = O n (1) independent of n, representing the average value of ∂ x s p (x) on the integration interval. The constant κ = O n (1) because for any r, r = O n (1), we have that

So, if s p (r) is isolated, the eigenvalue 1/r of L rw τ is isolated as well. On the other hand, if s p (r) is in the bulk, then ∃ q : |s p (r) -s q (r)| = o n (1). By an argument of continuity on s p (r), one can analogously define an r satisfying s p (r) = s q (r), concluding that in this case |r -r | = o n (1) and so eigenvalues in the bulk are mapped into eigenvalues in the bulk.

Appendix F. Fast Estimate of the ζ p

We here describe a fast algorithm to estimate the values of the ζ p using the Courant-Fisher theorem.

Let us denote by X T rt = (x T 1 (r t ), . . . , x T p (r t )) T ∈ R n×p where x p (r t ) is the eigenvector of H rt corresponding to the p-th smallest eigenvalue, s ↑ p (H rt ), while S rt = diag s ↑ 1 (H rt ), . . . , s ↑ p (H rt ) . For another value r = r t , and applying the Courant-Fischer theorem, we can write

i.e.:

We can further simplify the earlier expression by exploiting the identity:

We thus obtain

We now study the function f rt (r ) for r t ∈ (ζ p , ρ(B)) and define r t+1 ∈ (ζ p , r t ) as the solution (if it exists) to

The idea is to iteratively approach ζ p from the right and substitute r t ← r t+1 . If r t+1 ∈ (ζ p , r t ) as defined above exists, then {r t } t≥0 is a lower-bounded decreasing sequence: it thus converges to a limit r ∞ (potentially different from ζ p ). Exploiting (Rellich and Berkowitz, 1969, Theorem in I.4), denoting with A, B two Hermitian matrices and λ a scalar, s ↓ 1 (A + λB) is a convex function of λ. As a consequence, the function f rt (r ) is convex and has either no root or two roots. Since lim |r |→∞ f rt (r ) = +∞, it is enough to find a value of r for which f rt (r ) < 0 to prove that this function has two roots. With a straightforward computation, one can verify that f rt (r t ) < 0, so f rt (r ) has two roots, satisfying r t+1 < r t and r + > r t . By construction, since s ↑ p (H r ) is negative in the considered interval and s ↑ p (H r ) ≤ f rt (r )/r t , if f rt (r t+1 ) = 0, then ζ p ≤ r t+1 . Consequently, f rt (r t+1 ) = 0 has a unique solution satisfying ζ p ≤ r t+1 < r t and the algorithm converges to r ∞ = lim t→∞ r t ≥ ζ p . We are left to prove that r ∞ = ζ p . By convergence of r t , we have r t+1 -r t = o t (1). Plugging this relation solution into Equations (62,63), we obtain s ↓ 1 (r t+1 S r t+1 ) = r t+1 s ↑ p (H r t+1 ) = o t (1).