Spectral community detection in heterogeneous large networks

In this article, we study spectral methods for community detection based on $ \alpha$-parametrized normalized modularity matrix hereafter called $ {\bf L}_\alpha $ in heterogeneous graph models. We show, in a regime where community detection is not asymptotically trivial, that $ {\bf L}_\alpha $ can be well approximated by a more tractable random matrix which falls in the family of spiked random matrices. The analysis of this equivalent spiked random matrix allows us to improve spectral methods for community detection and assess their performances in the regime under study. In particular, we prove the existence of an optimal value $ \alpha_{\rm opt} $ of the parameter $ \alpha $ for which the detection of communities is best ensured and we provide an on-line estimation of $ \alpha_{\rm opt} $ only based on the knowledge of the graph adjacency matrix. Unlike classical spectral methods for community detection where clustering is performed on the eigenvectors associated with extreme eigenvalues, we show through our theoretical analysis that a regularization should instead be performed on those eigenvectors prior to clustering in heterogeneous graphs. Finally, through a deeper study of the regularized eigenvectors used for clustering, we assess the performances of our new algorithm for community detection. Numerical simulations in the course of the article show that our methods outperform state-of-the-art spectral methods on dense heterogeneous graphs.


Introduction and Motivations
The advent of the big data era is creating an unprecedented need for automating large network analysis. Community detection is among the most important tasks in automated network mining [1]. Given a network graph, detecting communities consists in retrieving hidden clusters of nodes based on some similarity metric (the edges are dense inside communities and sparse across communities). While quite simple to define, community detection is usually not an easy task and many methods arising from different fields have been proposed to carry it out. The most important of them are statistical inference, modularity maximization and graph partitioning methods. Statistical inference methods consist in fitting the observed network to a structured network model and infer its parameters (among which the assignment of the nodes to the communities) [2,3]. Modularity maximization algorithms rely on the modularity metric which quantifies the subdivision of networks into communities [1]. 1 The maximization of this quantity over all possible partitions in the graph gives the best possible subdivision of this graph in the modularity measure sense.
However, this is generally an NP-hard problem and many approximation methods have been proposed based on some polynomial-time heuristics: greedy methods [4], simulated annealing [5], extremal optimization [6] and spectral methods [7]. Spectral algorithms consist in retrieving the communities from the eigenvectors associated with the extreme eigenvalues of some matrix representation of the graph structure (adjacency matrix, modularity matrix, Laplacian matrix). By relaxing the modularity optimization problem from discrete values of the community memberships to continuous values, it is shown that approximate modularity maximization and even statistical inference methods can be performed via a low dimensional clustering of the entries of the dominant eigenvectors of the Laplacian matrix [8,9] in polynomial time. We focus in this article on the latter methods. Precisely, spectral methods for community detection generally follow successively those steps 1. Compute the, say, eigenvectors corresponding to the extreme (largest or smallest) eigenvalues of one of the matrix representations of the network (adjacency, modularity, Laplacian).
2. Stack those eigenvectors column-wise in a matrix W ∈ R n× with n the number of nodes, or correspondingly the size of the matrix representation of the network. 3. Take each row of W as a (feature) vector in a -dimensional (feature) space. 4. Cluster those n vectors in K groups using a standard classification algorithm e.g., k-means or expectation maximization (EM). The EM algorithm, for example, aims to roughly identify clusters at first before to sequentially update the individual cluster means and covariances until convergence.
Real world networks are in general sparse in that the number of connections of each node (degree) scales in O(1) when the number of nodes n grows large. When the degrees scale instead like O(log n) or O(n), the network is said to be dense. The standard spectral algorithms based on the network matrix (adjacency, modularity, Laplacian) of sparse graphs are generally suboptimal in the sense that they fail to detect the communities down to the transition point where the detection is theoretically feasible [10]. New operators (non-backtracking [10], Bethe Hessian [11]) based on statistical physics have recently been proposed and have been shown to perform well down to the aforementioned sparse regime. We focus however in this article on dense networks for which spectral methods are often optimal [10].
Most of the works proposing statistical analysis of the performance of community detection (for dense as well as sparse networks) consider the basic Stochastic Block Model (SBM) as a model for networks decomposable into communities. Denoting G a K-class graph of n vertices with communities C 1 , . . . , C K with g i the group assignment of node i, the SBM assumes an adjacency matrix A ∈ {0, 1} n×n , with A ij independent Bernoulli random variables with parameter P gigj where P ab represents the probability that any node of class C a is connected to any node of class C b .
The main limitation of this model is that it is only suited to homogeneous graphs where all nodes have the same average degree in each community (besides, class sizes are often taken equal). A more realistic model, the Degree-Corrected SBM (DCSBM), was proposed in [12,13] to account for degree heterogeneity inside communities. For the same graph G defined above, by letting q i , 1 ≤ i ≤ n, be some intrinsic weights which affect the probability for node i to connect to any other network node, the adjacency matrix A ∈ {0, 1} n×n of the graph generated by the DCSBM is such that A ij are independent Bernoulli random variables with parameter q i q j C gigj , where C gigj is a class-wise correction factor.
The main motivation of this work arises from the observation that classical spectral algorithms based on the adjacency matrix (modularity, Laplacian) may drastically fail to detect the genuine communities in some synthetic graphs generated using the DCSBM. The same observation is made even for the aforementioned and very competitive Bethe Hessian (BH) method. 2 To illustrate those limitations of spectral methods under the DCSBM, the two graphs of Figure 1 provide 2D representations of dominant eigenvector 1 versus eigenvector 2 for the standard modularity matrix and the BH matrix, when half the nodes connect with low weight q (1) and half the nodes with high weight q (2) . For both methods, it is clear that k-means or EM alike would erroneously induce the detection of extra communities and even a confusion of genuine communities in the BH approach.
We have come to understand that those extra communities are produced by some biases created by the heterogeneity of the intrinsic weights q i 's; intuitively, nodes sharing the same intrinsic connection weights tend to create their own sub-cluster inside each community, thereby forming additionnal sub-communities inside the genuine communities.
In order to understand the aforementioned limitations and the different mechanisms into play when using spectral methods based on matrices derived from the adjacency matrix A (such as the modularity matrix), we study here a generalized version of the normalized modularity matrix, 2 The Bethe Hessian (BH) spectral method [11] is based on the union of the eigenvectors associated to the negative eigenvalues of H(rc) and H(−rc) respectively where H(r) = (r 2 − 1)In − rA + D with In the identity matrix of size n, D a diagonal matrix containing the degrees on the diagonal and rc = i d 2 given for α ∈ A ⊂ R, by where d is the vector of degrees (d i = n j=1 A ij ), D is the diagonal matrix of degrees (containing d on the main diagonal) and m = 1 2 d T 1 n is the number of edges in the network. In particular, L 0 is the modularity matrix, L 1 2 is a modularity equivalent to the Laplacian matrix and L 1 was studied in [14,15].
In the dense DCSBM model where q i = O(1) (with growing n), when the correction factors C gigj differ by O(1), the classification is trivial as asymptotically vanishing error rates are easily guaranteed. We thus place ourselves in the more interesting regime where C gigj = O (1) individually but the C gigj 's differ by O(n − 1 2 ). We study the dominant eigenvalues and associated eigenvectors (used for classification) of L α for large dimensional dense graphs following the DCSBM in the aforementioned "non-trivial" regime.
In a nutshell, our main findings are as follows • We show that, as n → ∞, L α can be arbitrarily well approximated by a theoretically tractable random matrixL α which falls in the family of so-called spiked random matrix models and which allows for a thorough understanding of spectral methods based on L α .
Those random matrices generally exhibit a phase transition beyond which useful information can be extracted from the eigenvectors associated to outlying eigenvalues (and below which nothing can be said). In our context, this phase transition corresponds to a community detectability threshold, common in community detection algorithms analysis. We characterize exactly this phase transition for each value of α.
• We prove the existence of an optimal value α opt of α for which the aforementioned phase transition is maximally achievable.
• We provide a consistent estimatorα opt of α opt based on d alone.
• A thorough analysis of the spiked random matrix model then shows that, to achieve consistent clustering in the DCSBM model, the dominant eigenvectors used for clustering should be pre-multiplied by D α−1 prior to the low dimensional classification (step (4) of the spectral algorithm described previously). 3 • A deeper study of those regularized eigenvectors allows us to improve the initial setting of the EM algorithm (in the step (4) of the spectral algorithm described above) in comparison with a random setting.
find the theoretical clustering error rate of spectral community detection in the regime under study.
• Numerical simulations (throughout the article) show that our methods outperform state-ofthe-art techniques both on synthetic graphs and on real world networks.
Notations: Vectors (matrices) are denoted by lowercase (uppercase) boldface letters. {v a } n a=1 is the column vector v with (scalar or vector) entries v a and {V ab } n a,b=1 is the matrix V with (scalar or matrix) entries V ab . For a vector v, the operator D(v) = D ({v a } n a=1 ) is the diagonal matrix having the scalars v a down its diagonal and for a matrix V, D(V) is the vector containing the diagonal entries of V. The vector 1 n ∈ R n stands for the column vector filled with ones. The Dirac measure at x is δ x . The vector j a is the canonical vector of class C a defined by (j a ) i = δ i∈Ca and J = [j 1 , . . . , j K ] ∈ {0, 1} n×K . The set C + is {z ∈ C, [z] > 0}.

Preliminaries
This section describes the network model under study, which is based on the DCSBM defined in the previous section, and provides preliminary technical results.
We consider an n-node random graph with K classes C 1 , . . . , C K of sizes |C k | = n k . Each node is characterized by an intrinsic connexion weight q i which affect the probability that this node gets attached to another node in the graph. A null model would consider that the existence of an edge between i and j has probability q i q j . In order to take into account the membership of the nodes to some group, we define C ∈ R K×K as a matrix of class weights C ab , independent of the q i 's, affecting the connection probability between nodes in C a and nodes in C b .
As such, following [13], the adjacency matrix A of the graph generated from a DCSBM model has independent entries (up to symmetry) which are Bernoulli random variables with parameter P ij = q i q j C gigj ∈ (0, 1) where g i is the group assignment of node i. We set A ii = 0 for all i. In the dense regime under consideration, q i = O(1) and C gigj = O(1) as n → ∞. For convenience of exposition and without loss of generality, we discard the nodes having no neighbor and we assume that node indices are sorted by clusters i.e, nodes 1 to n 1 constitutes C 1 , nodes n 1 + 1 to n 1 + n 2 form C 2 and so on.
The matrix under study is given by for some α ∈ A, a compact subset of R, where d = A1 n , D = D(d) and m = 1 2 d T 1 n . As stated in the introduction, we are mainly interested in a dense network regime where clustering is not asymptotically trivial. This regime is ensured by the following growth rate conditions. Assumption 1. As n → ∞, K remains fixed and, for all i, j ∈ {1, . . . , n} random variables with measure µ having compact support in (0, 1).

ni
n → c i > 0 and we will denote c = {c k } K k=1 (in particular, c T 1 K = 1).
Before delving into the main technical results, we provide a uniform consistent estimator of the (a priori unknown) intrinsic weight q i which shall be used in the course of the article. .
Note that q i can be retrieved from the empirical graph degrees irrespective of the class matrix C, which is a direct (and important) consequence of Assumption 1-(1).
The first goal of the article is to study deeply the eigenstructure of L α . As can be observed, L α has non independent entries as D (and d) depend on A, and it does not follow a standard random matrix model. Our strategy is to approximate L α by a more tractable random matrix which asymptotically preserves the eigenvalue distribution and isolated eigenvectors of L α . Before providing the complete proof in Section 5.1, let us give the main steps of the approximation. As per our model and assumptions, the random variable A ij is Bernoulli distributed with parameter P ij = . We may thus write . The normalized adjacency matrix is thus 4 where q (i) = [q n1+...+ni−1+1 , . . . , q n1+...+ni ] T ∈ R ni (n 0 = 0) and X = {X ij } n i,j=1 . Note that the right-hand side of (4) is composed of a dominant (in terms of operator norm) matrix A d, and of smaller order terms. From there, we may then provide a Taylor By grouping all those expansions consistently following the structure of Equation (2) and by only keeping non-vanishing operator norm terms, we obtain the corresponding approximate of L α as follows Theorem 1. Let Assumption 1 hold and let L α be given by (2). Then, for D q = D(q), as n → ∞, L α −L α → 0 in operator norm, almost surely, wherẽ As far as the spectral analysis is concerned,L α is asymptotically equivalent to L α , as they asymptotically share the same set of eigenvalues and isolated eigenvectors. 5 Thus, for large enough n, the spectral analysis of L α can be performed through that ofL α . Interestingly,L α is an additive spiked random matrix [16] as it is the sum of a standard full rank random matrix n − 1 2 D −α q XD −α q (symmetric matrix having independent entries of zero mean and O(n −1 ) variances) and a low rank matrix UΛU T . As shown in Figure 2, the spectrum (eigenvalue distribution) of spiked random matrices is generally composed of (one or several) bulks of concentrated eigenvalues and, when a phase transition is met, of additionnal eigenvalues which isolate from the aforementioned bulks. The eigenvectors associated to those isolated eigenvalues contain important information related to the low rank matrix up to some noise; the more those eigenvalues isolate from the bulks, the lesser noise is contained in the corresponding eigenvectors. More specifically, the eigenvectors of the spiked random matrix become more correlated to the eigenvectors of the low rank matrix as the isolated eigenvalues are far away from the phase transition threshold.
From Theorem 1, we see that the low rank matrix UΛU T contains the matrix D 1−α q J; so in our case, when the phase transition is met, the eigenvectors ofL α will be correlated to some extent to D 1−α q J. But, for a consistent clustering, one expects instead the vectors used for classification to be correlated to the canonical vectors j a , 1 ≤ a ≤ K. An important outcome of this first preliminary result is thus that the eigenvectors of L α should be pre-multiplied by D α−1 prior to the classification. 6 This first result helps correcting the biases (creation of artificial classes) introduced by the degree heterogeneity when using classical spectral methods (as observed earlier in Figure 1). As shown in Figure 3, which assumes the same setting as Figure 1, when the aforementioned eigenvector regularization is performed prior to EM or k-means classification, the genuine communities are correctly recovered.
As is classical in the analysis of spiked random matrices, our next task is to study the isolated eigenvalues of L α and their associated eigenvectors. This will in particular allow us to i) evidence the phase transition phenomenon discussed earlier which corresponds here to a community detectability threshold and ii) evaluate the per-class average means and average covariances of 6 As far as the eigenvectors are concerned, we may freely replace Dq (unknown in practice) by D (which can be computed from the observed graph) since, from Lemma 1, the vector of degrees d is, up to the scale factor   the eigenvectors used for clustering, thereby leading to the clustering performances. This is the objective of the next sections.

Eigenvalues
In this section, we are interested in the localization of the eigenvalues of L α . Since L α is asymptotically equivalent to a spiked random matrix, its eigenvalues are expected to be asymptotically the same as the eigenvalues of the full rank "noise" matrix n − 1 2 D −α q XD −α q which are essentially concentrated in bulks, but possibly for finitely many of them which can isolate from those bulks when some eigenvalues of the low rank matrix UΛU T are sufficiently large [17].
To study the eigenvalues and eigenvectors of the spiked random matrix L α , we follow standard random matrix approaches [18] and [19]. We will first determine the support of the limiting eigenvalue distribution of where most of the eigenvalues of L α concentrate (in bulks). Then, we will find the plausible isolated eigenvalues of L α only induced by the low rank matrix UΛU T and which lie outside the aforementioned main support.
For future use, we further define for z,z ∈ C \ S α e α ab;2 (z,z) = and e α ab;3 (z,z) = . We have in particular Figure 2 for illustration.
Remark 2 (Semi-circle law). For homogenous graphs where ∀i, which is the Stieltjes transform of the well known semi-circle probability measure 7 with support S α = [−2, 2] and density given byπ By Proposition 1, we can define for ρ / ∈ S α and for large n almost surely the resolvent Q α ρ = (n − 1 2 D −α q XD −α q − ρI n ) −1 . As per Theorem 1, since the hypothetically isolated eigenvalues of L α are to be found outside the limiting support S α , we thus need to find those ρ's at a non-vanishing . This leads to solving, for large n, det(I K+1 + U T Q α ρ UΛ) = 0. We then show that U T Q α ρ UΛ converges almost surely to a deterministic matrix and we have the following result Theorem 3 (Isolated Eigenvalues). Let Assumption 1 hold and, for z ∈ C \ S α (given in Note that the limiting value E α 0 (z) does not depend on α in this case.
We shall denote byρ such eigenvalues when they exist.
• 1 + θ α (ρ) = 0: the left and right eigenvectors associated to the zero eigenvalues of G α ρ are respectively orthogonal to the right and left eigenvectors associated to the non-zero eigenvalues. So, by letting V l , V r be matrices containing in columns the respectively left and right As we show in Appendix Appendix D, for 1 + θ α (ρ) = 0, the eigenvectors associated to the aforementioned isolated eigenvaluesρ will not contain information about the classes. This case is thus of no interest for clustering. It is nevertheless important from a practical viewpoint to note that, even in the absence of communities, spurious isolated eigenvalues may be found that may deceive the experimenter in suggesting the presence of node clusters. From now on, we will only consider the isolated eigenvalues ρ for which 1 + θ α (ρ) = 0.
Since it is more convenient to work with symmetric matrices (having identical left and right eigenvectors), the following remark will be useful in what follows.
Remark 4 (Informative eigenvectors). The next three statements are equivalent. For ρ a limiting isolated eigenvalue of L α such that 1 + θ α (ρ) = 0, is a set of left/right eigenvectors of G α ρ associated to zero eigenvalues.
As L α is asymptotically equivalent (through Theorem 1) to a spiked random matrix, the eigenvectors of L α associated to isolated eigenvalues are expected to correlate (to some extent) to the eigenvectors of UΛU T (defined in Theorem 1) and thus to D 1−α q J. Clustering based on the eigenvectors of L α should then be possible when they are associated to such isolated eigenvalues.
Based on this fact and with the help of Theorems 2 and 3, we can exactly characterize the phase transition threshold beyond which community detection is indeed possible. .

Remark 5 (Rank ofM and maximum number of isolated eigenvalues). From Corollary 1,
there is a one-to-one mapping between isolated eigenvalues ρ of L α and non zero eigenvalues of M = D(c) − cc T M. As 1 T KM = 0,M has a maximum of K − 1 non zero eigenvalues which means that at most K − 1 eigenvalues of L α can be found at macroscopic distance from S α (excluding the case 1 + θ α (ρ) = 0). Thus, at most K − 1 eigenvectors of L α can be used in the first step of the spectral algorithm described in the introduction.

Application 1: Optimal α
In this section, we find the α for which the community detectability threshold is maximally achieved. This, in turn, shall allow to extract some non-trivial information about the classes from the extreme eigenvectors.
From Corollary 1, sinceM does not depend on α, the smaller τ α the more likely the detectability condition λ(M) > τ α is met. We then seek α for which τ α is minimal. We may thus define α opt ≡ argmin α∈A {τ α } .
Retrieving α opt has a tremendous practical advantage as it optimizes the detection of barely detectable communities (and, as shall be seen through simulations, greatly improves the performance). The estimation of α opt however requires the knowledge of E α 2 (x) for each α ∈ A. The estimation of E α 2 (x) can be done numerically by solving the fixed point equation defined in Theorem 2 provided µ is known. Thanks to Equation (3), µ can be consistently estimated from the empirical degrees d i 's. We thus have all the ingredients to estimate α opt .
Remark 6 (Estimation of S α ). To estimate numerically the right edge S α of the support S α = To this end, we evaluate S α by an iterative dichotomic search in intervals of the type [l, r] for which E α 2 (l) is undefined (and thus the system of fixed point equations defining E α 2 (x) in Theorem 2 does not converge for x = l) and E α 2 (r) is defined (the system of fixed point equations defining E α 2 (x) in Theorem 2 converges for x = r), starting from e.g., l = 0 and r quite large. Since for the points x where E α 2 (x) is not defined, the fixed point algorithm solving the equations defining E α 2 (x) run indefinitely, one must fix a given number of iterations after which the algorithm is stopped and x is then considered to be in S α .
The aforementionned importance of choosing α =α opt along with the need to pre-multiply the dominant eigenvectors of L α by D α−1 before classification, as discussed after exposing Theorem 1, naturally bring us to a novel heterogeneous-graph improved community detection method.
Algorithm 1 below summarizes the main steps of our improved spectral algorithm. Note that the same algorithm can be applied for any α by replacingα opt by the corresponding value of α and skipping step 1.
Algorithm 1: Improved spectral algorithm 1: Evaluate α =α opt = argmin α∈(0,1) lim x↓Ŝ αÊ α 2 (x) as per Lemma 2. 2: Retrieve the eigenvectors corresponding to the largest eigenvalues (which are found away from the bulks in the spectrum of Let r 1 , . . . , r n ∈ R be the rows of W. Cluster r i ∈ R , 1 ≤ i ≤ n in one of the K groups using any low-dimensional classification algorithm (e.g., k-means or EM). The label assigned to r i then corresponds to the label of node i. 1 αopt Figure 4: Ratio between the largest eigenvalue λ of Lα and the right edge of the support S α , as a function of the largest eigenvalue ofM, In the following, we will restrict ourselves to α ∈ A = [0, 1] for the numerical simulations.
To illustrate the importance of the choice of α opt , Figure 4 presents the theoretical (asymptotic) ratio between the largest eigenvalue of L α and the right edge of the limiting support S α (S α + ) with respect to the amplitude of the eigenvalues ofM. Intuitively, the larger this ratio the better the clustering performance as the eigenvector associated to this largest eigenvalue contains less noise.
Although α opt only ensures in theory to have the best isolation of the eigenvalues only in "worst cases scenarios"(i.e., when λ(M) is only slighty larger than τ αopt ), Figure 4 shows that taking α = α opt provides the largest gap λi(Lα) S α + for all values of λ(M). This suggests (again, without any theoretical support) better performances with α = α opt in all cases (for any value of M).
In the sequel, to compare the different algorithms, we will use the performance evaluation measure known as the overlap to ground truth communities, defined in [10] as where g i andĝ i are the true and estimated labels of node i, respectively. Note that this definition implicitly suggests that all communities are of equal proportions and is therefore not fully compatible with our more general present setting which allows for unbalanced classes. Figure 5 subsequently shows the overlap performance under the setting of Figure 4. It is worth mentioning that the empirically observed phase transitions closely match the theoretical ones (drawn in circles and the same as in Figure 4). We then present in Figure 6 an example where the BH 9 corresponds in practice to the number of eigenvalues of Lα which isolate from the bulks. From Remark 5, we can have up to K − 1 of them.  algorithm fails due to strongly heterogeneous node degrees. Assuming nodes connect with either low q (1) = 0.1 or high q (2) > q (1) intrinsic weights, we observe a sudden drop of the BH overlap for large q (2) − q (1) . This phenomenon is consistent with the fact, observed earlier in Figure 1, that BH creates artificial communities out of nodes with the same q i parameter. This is a practical demonstration of the need for a proper eigenvector normalization to avoid degree biases. This observation has recently led [20] to consider a regularization for the non-backtracking operator on which the BH method is based.  In Figure 7, we consider a more realistic synthetic graph where the q i 's assume a power law of This mismatch is likely due to the fact that our theoretical results in this article require P ij = O (1) which is not always the case in this scenario.  We finally confront the overlap performance on real world benchmarks in Table 1. The best overlap score for each benchmark is set in boldface and quasi-equal scores are shown in italic.
Our approach largely outperforms the BH method on some benchmarks and has competitive performances on others. However note that, for so small network sizes, the performance achieved by Lα opt may be quite unsatisfactory. 10 10 We should note here that the scores for the BH are different from the ones found in the article [11] since here we are running a consistent algorithm (EM in the last step of the spectral algorithm) for all K while the authors of [11] have instead used the signs of the eigenvector for networks with two communities and k-means algorithm for those with more than two communities.
In order to assess the performance of Algorithm 1, we now need to investigate closely the content of the eigenvectors of L α used for clustering (and of their pre-multiplied by D α−1 versions).
These regularized eigenvectors happen to be shapped like noisy "plateaus" (step functions), each plateau characterizing a class. The properties of those noisy plateaus are significant to assess the performance of spectral clustering. The objective of the next section is to characterize exactly those quantities.

Eigenvectors
In this section, in order to fully characterize the performances of Algorithm 1, we study in depth the normalized eigenvectorsv α i used for the classification in the algorithm (step 3 of Algorithm 1). The eigenvectors corresponding to the eigenvalues for which 1 + θ α (ρ) = 0 are still not considered since we recall they are of no interest for the classification. For technical reasons, we restrict As one can see in Figure 3, the different clusters of points (rows of W in Algorithm 1) have different dispersions (variances) in the DCSBM model under consideration. The most appropriate algorithm to use in step 4 of Algorithm 1 is the expectation maximization (EM) method. EM considers each point r i ∈ R arising from [v α 1 , . . . ,v α ] as a mixture of K Gaussian random vectors with means ν a EM and covariances Σ a EM ∈ R × , a ∈ {1, . . . , K}. Starting from initial means and covariances, they are sequentially updated until convergence. To identify ν a EM , Σ a EM and thus understand the performance of Algorithm 1, we may writev α i 11 as the "noisy plateaus" vector where w a i ∈ R n is a random vector orthogonal to j a , of norm √ n a and supported on the indices of C a and with D a = D(j a ). The vector ν a = (ν a i ) i=1 ∈ R and the matrix Σ a = (σ a ij ) i,j=1 ∈ R × represent respectively the empirical means and empirical covariances of the points r i (defined in Algorithm 1) belonging to class C a . Thus, provided that EM converges to the correct solution, (ν a EM ) i and (Σ a EM ) ij shall converge asymptotically to the limiting values of ν a i ∈ R and σ a ij respectively. 11 Recall that the graph nodes were assumed labeled by class, and thus the entries ofv α i are similarly sorted by class.
Clearly, for small values of Σ a compared to ν a , clustering the vectorsv α i shall lead to good performances. We shall investigate here the links between Σ a , ν a and the model parameters M, µ and α.
Technically, the standard tools used in spiked random matrix analysis do not allow for an immediate assessment of the quantities ν a i and σ a ij . As a workaround, we follow the approach used in [25] which relies on the possibility to estimate bilinear forms of the type a T u α i (u α i ) T b for given vectors a, b ∈ R n and unit multiplicity eigenvectors u α i of L α as we have from Cauchy formula, as n → ∞ almost surely, and for a given matrix D where Γ ρ is a positively oriented contour circling around the limiting eigenvalue ρ of λ i (L α ) associated to the eigenvector u α i of L α . Similar to the eigenvalue estimation step in previous sections, the estimation of those quantities then consists in relating L α toL α and then to G α z . Precisely, • we will estimate the ν a i 's by obtaining an estimator of the K × K matrix 1 n the diagonal entries of which allow to estimate |ν a i | while the off-diagonal entries are used to decide on the signs of the ν a i 's (up to a convention in the sign of u α i ).
• Similarly, we may first estimate the more involved object 1 n can be retrieved by dividing any entry e, f of the former quantity by non-vanishing quantities ν e i ν f i . For the eigenvectors u α i used for clustering, there is always at least one index f such that ν f i is non zero (otherwise, this eigenvector is of no use for clustering).
All calculus done, we obtain the following limit for the empirical means ν a i 's.
of unit multiplicity, mapped to eigenpair (ρ, u α i ) of L α as defined in Corollary 1, under the conditions of Assumption 1 and for ν a i defined in (9), we have almost surely as n → ∞, And for the empirical covariances σ a ij 's, we have the following limit 2 mapped respectively to (ρ 1 , u α i ) and (ρ 2 , u α j ) eigenpairs of L α and for σ a ij defined in (10), we have almost surely as n → ∞, where χ α (ρ) is defined in Theorem 4.
We thus obtain from Theorems 4 and 5 along with Remark 7, Corollary 2 (Means and covariances for K = 2 classes). For a = 1, 2 for ρ the unique isolated eigenvalue of L α (if it exists).
Remark 8 (Case q i = q 0 ). Here, since ∀i, q i = q 0 , all the e ij (ρ)'s are completely explicit and only depend on q 0 , ρ and α while for q i = q 0 , the e ij (ρ)'s can only be found by solving a fixed point equation numerically. Furthermore, e 00;2 (ρ 1 , ρ 2 ) = e 00 (ρ 1 )e 00 (ρ 2 ) and thus This indicates that the normalized eigenvectorsv α i ∈ R n , 1 ≤ i ≤ used for the classification (Step 4 of Algorithm 1) are asymptotically uncorrelated. Thus, the classification can be performed by treating each eigenvector independently, rather than jointly. In addition, we have σ a,∞ ii = σ b,∞ ii for 1 ≤ a = b ≤ K meaning that the nodes can only be discriminated based on the means ν a i 's. The performances of EM highly depend on the chosen starting parameters; a first natural choice is to set them randomly, which as we shall see leads to poor performances especially in cases where the clusters are not easily separable.
Since the theoretical limiting means ν a,∞ and covariances Σ a,∞ are respectively the limiting values of ν a EM and covariances Σ a EM provided EM converges to the correct solution, we may set as initial parameters of EM our findings ν a,∞ (Theorem 4) and Σ a,∞ (Theorem 5) for a ∈ {1, . . . , K} if those can be estimated. It will be confirmed later via simulations that this new setting of initial parameters is much better than a random initial setting of EM, and is close to the theoretical performances evaluated in Section 3.4.2. Unfortunately, in most scenarios, the many unknowns prevent such an estimation. Nonetheless, from Corollary 2, provided the class proportions (or the sizes of each class) are (more or less) known, we can consistently estimate ν ∞ and Σ ∞ in a 2-class scenario. To show the effect of our setting of initial parameters of EM based on the findings ν ∞ and Σ ∞ , Figure 8 compares the empirical performances of our new spectral algorithm based on the regularized eigenvector of L 0.5 for different initial settings of the EM parameters i) random setting (Random EM) ii) our theoretical setting (by assuming that the class proportions are known) and iii) the ground truth setting (oracle EM where we set the initial points to the empirically evaluated means and covariances of each cluster based on ground truth) versus the theoretical ones (evaluated in Section 3.4.2). Below the transition point, no consistent clustering can be achieved for large n using the eigenvectors associated to highest eigenvalues since the clusters are not separable and our theoretical limiting means and covariances are not defined since there is no isolated eigenvalues in that case. We have thus initialized EM at random in this non interesting regime (as for Random EM). The EM algorithm may in that regime set all the nodes to the same cluster, which will then result to a classification rate close to the proportion of the nodes in the cluster of largest size. This explains the important mismatch between theory and practice (the retrieval of the theoretical plot is discussed in the subsequent section) below the transition threshold, when c 1 and c 2 are not equal

Performance analysis
In this section, we analyze the performance of our new spectral community detection algorithm in terms of probability of misclassification.
From our theoretical class means and class covariance findings, we may compute the theoretical misclassification error probability by defining the proper decision regions which separate the different clusters in -dimensional space. For simplicity, we focus here on K = 2 so that = 1, but the generalization to K > 2 is straightforward.

Concluding Remarks
In this article, we have studied a family of graph affinity matrices L α ∝ 1 √ n D −α A − dd T 2m D −α which generalize the matrices (modularity, Laplacian) used for spectral community detection in dense networks. The main difficulty for the study of those random matrices comes from the dependency between their entries. We tackle this difficulty by establishing the approximation L α −L α → 0 using similar techniques as in [26] whereL α belongs to the class of so called "spiked" random matrices for which the study of eigenvalues and eigenvectors is classical. The study of eigenvalues and eigenvectors of L α used for the classification is thus performed usingL α .
We go further than the observation of [14] and [20] which state that it is important to use the eigenvectors of L 1 rather than the classically used L 0 for the classification when the network has heterogeneous degree distribution to avoid some important misclassifications induced by degree biases. We saw in Figure 3 for example that the eigenvectors of L 1 correct the degree biases but we show that it is better to use instead the eigenvectors of L 0 premultiplied by D −1 (see Figures 4-7 for example). Better still, we show that there exists an optimal α called α opt for which taking the eigenvectors of L αopt pre-multiplied by D αopt−1 ensures best performance (or to be more precise best asymptotic cluster detectability).
We generalize the study in [27] concerning the evaluation of per-class means of the entries of the unique eigenvector used for the classification, which was limited to the symmetric stochastic block The results and methods in this article are all based on the strong assumption that the classwise correction factors C gigj differ by O(n − 1 2 ) since ∀i, j ∈ {1, . . . , n}, C gigj = 1 + Mg i g j √ n . When this condition is not fulfilled, our results may not be applicable. However, this is a more interesting assumption, challenging in practice since when the C ab 's differ by O(1), one can easily achieve a vanishing misclassification error rate asymptotically and most algorithms perform similarly well in detecting the communities. Besides, our analyses assume very large networks (n → ∞ with convergence rates of the order n − 1 2 ) and thus the performances of our algorithm on small size graphs can be quite poor as observed in Table 1.
Our broad study of spectral methods for community detection is so far limited to dense networks. Real world networks being sparse in general, it would be interesting to extend our study to sparse graph models. The Non Backtacking and Bethe Hessian (BH) methods developped under homogeneous graph models, are currently state-of-the-art spectral methods which allow for non-trivial performances of community detection in sparse graphs. For heterogeneous graphs, the "flow" matrix (the Non Backtracking matrix normalized by the degrees) proposed in [20] is shown to have a much better behavior than the two aforementionned matrices in sparse heterogeneous graphs. This follows the same spirit as what is done in this article where we normalize the modu-larity matrix by the degrees for dense heterogeneous graphs. The use of the symmetric BH matrix (of size n the number of nodes) is highly preferred in terms of computational time and memory than the non-symmetric, high-dimensional Non Backtracking matrix (of size 2m with m the number of edges). Following then [20] and our current work, one may wonder about the existence of an equivalent α-normalized BH matrix and its associated performances. This is, we believe, a promising avenue of future investigation.
Pushing further the applicability reach of the present study, note that one of the remaining issues of spectral clustering methods especially for large graphs is the expensive computational complexity of the eigenvectors of the large random matrices representing those graphs. This computation burden is reduced when using power methods and is thus less expensive for sparse graphs (with many zeros in the corresponding matrix). This suggests potential computational gains incurred by smartly removing some edges in the graph to make it sparse prior to eigenvectors computation, which will of course reduce the performance to the benefit of the computational cost.
Our mathematical framework may allow us to study the tradeoff between complexity/cost and performances of spectral methods for community detection on such subsampled large dense graphs.

Preliminaries
The random matrix under study L α = (2m) α 1 √ n D −α A − dd T 2m D −α is not a classically studied matrix in random matrix theory. We will thus first find an approximate tractable random matrixL α which asymptotically preserves the eigenvalue distribution and the extreme eigenvectors of L α (Section 5.1). Then, the empirical distribution of the eigenvalues of L α is studied in Section 5.2 along with the exact localization of those eigenvalues in Section 5.3. Finally, a thorough study of the eigenvectors associated to the aforementioned eigenvalues is investigated in Section 5.4.
We follow here the proof technique of [25]. In the sequel, we will make some approximations of random variables in the asymptotic regime where n → ∞. For the sake of random variables comparisons, we give the following stochastic definitions. For x ≡ x n a random variable and Most of the proofs here are classical in random matrix theory (see e.g., [17]) but require certain controls inherent to our model. The goal of the article not being an exhaustive development of the proofs techniques, we will admit a number of technical results already studied in the literature.
However, we will exhaustively develop the calculus to obtain our final results which are not trivial.

Random equivalent for
D −α has non independent entries and is not a classical random matrix model. The idea is thus to approximate L α by a more tractable random matrix modelL α in such a way that they share asymptotically the same set of outlying eigenvalues/eigenvectors which are of interest in our clustering scenario. We recall that the entries A ij of the adjacency matrix were defined from the DCSBM model as independent Bernoulli random variables with parameter q i q j 1 + ; one may thus write where X ij , 1 ≤ i, j ≤ n, are independent (up to the symmetry) zero mean random variables of . We can then write the normalized adjacency matrix as follows where 12 q (i) = [q n1+...+ni−1+1 , . . . , q n1+...+ni ] T ∈ R ni (n 0 = 0) , X = {X ij } n i,j=1 and D q = D(q). The idea of the proof is to write all the terms of L α based on Equation (12), since all those terms depend on A. To this end, we will evaluate successively d = A1 n , D = D(d), dd T and 2m = d T 1 n .
It will appear that D and d T 1 n are composed of dominant terms (with higher operator norm) and vanishing terms (with smaller operator norm); we may then proceed to writing a Taylor expansion of D −α and (2m) α = (d T 1 n ) α for any α around their dominant terms to finally retrieve a Taylor expansion of L α .
Let us start by developing the degree vector d = A1 n . We have . (13) Let us then write the expansions of d T 1 n , (d T 1 n ) α , dd T and dd T (d T 1n) respectively. From (13), we obtain Thus for any α, proceeding to a 1 st order Taylor expansion, we may write Besides, from (13) we have Keeping in mind that we shall only need terms with non vanishing operator norms asymptotically, to have terms with spectral norms of order at least O(1). We get from multiplying (16) and (15) By subtracting (17) from (12), we obtain It then remains to evaluate D −α . From (13), we may write D = D(d) as .
The right hand side of D (in brackets) having a leading term in O(1) and residual terms in O(n − 1 2 ), the Taylor expansion of the (−α)-power of D is then retrieved By combining the expressions (15), (18) and (19), we obtain a Taylor approximation of L α as follows The three following arguments allow to complete the proof • 1 n = J1 K and D q 1 n = q. −→ m µ c where we recall that m µ = tµ(dt).
• As X is a symmetric random matrix having independent entries of zero mean and finite variance, from the law of large numbers, we have 1 Using those three arguments, L α may be further rewritten By rearranging the terms of (20), we obtain the expected result This proves Theorem 1.

Limiting spectral distribution of L α
It follows from Theorem 1 thatL α = D −α q X √ n D −α q + UΛU T is equivalent to an additive spiked random matrix [28] where with the difference that the deterministic part UΛU T is not independent of the random part (an issue that we solve here) and U is not composed of orthonormal vectors. Let us then studyX = D −α q X √ n D −α q (having entriesX ij with zero mean and variance σ 2 ij /n with σ 2 ij = q i q j (1 − q i q j ) + O(n − 1 2 )) and show that its e.s.d.π α converges weakly toπ α with Stieljes transform E α 0 (z) = (t − z) −1 dπ α (t) for z ∈ C + defined in Theorem 2. This will imply (By Weyl interlacing formula) that the empirical spectral measure π α ≡ 1 n n i=1 δ λi(Lα) (with λ i (L α ) eigenvalues ofL α ) will also converge toπ α .
The matrixX is a classical random matrix model in RMT already studied in similar cases. It is well known for those random matrix models (having entries with given means, variances and bounded first order moments) that the law of theX ij 's does not change the results on the limiting law of the e.s.d.π α : this property is kwown as universality (e.g., [29]). For technical reasons, we can thus assume that theX ij 's are Gaussian random variables with the same means and variances in order to use standard Gaussian calculus, introduced in [30]. The objective of the proof is to find the deterministic limit E α 0 (z) for the random quantity 1 n tr X − zI n −1 which is the Stieljes transform of the e.s.d.π α . Deterministic equivalents for the Stieljes transform of empirical spectral measures associated with centered and symmetric random matrix models with a variance profile have already been studied in for example [31,32]. We give in Appendix Appendix B an exhaustive development of the Gaussian calculus to obtain E α 0 (z). The final result is as follows. From Lemma 3, we get directly 1 from which we get E α 0 (z) = e α 00 (z) where From this, we have that E α 0 (z) does not depend on n, so that 1 n tr Q a.s.
Remark 9 (Convergence of the e i 's.). Similar results to Lemma 3 have been derived for example in [32] and the fixed point algorithm (21) which consists of iterating the e i 's is shown to converge. Since the calcultation of the e ab 's is an intermediary step of (21) from (22), the fixed point algorithm (23) also converges. Similarly to [32], when none of the (D −α q ) ii 's is isolated, the random matrixX does not produce isolated eigenvalues outside the support S α ofπ α . Here, for large n, this property is verified since from Assumption 1, the q i 's are i.i.d. arising from a law with compact support (the probability that a (D −α q ) ii gets isolated tends to 0 asymptotically). This gives Proposition 1 which we will not prove here; similar proofs are provided for example in [33].
From the analyticity of the Stieljes transform outside its support, Lemma 3 extends naturally to C \ S α . This proves Theorem 2.

Isolated eigenvalues of L α
In the previous section, we have shown that the e.s.d. of L α converges weakly to the limiting law of the eigenvalues ofX since they only differ by a finite rank matrix. We shall have in addition isolated eigenvalues of L α induced by the aforementionned low rank matrix. We are interested here in the localization of eigenvalues of L α isolated from the support S α of the limiting law of its e.s.d. According to Proposition 1, there is almost surely no eigenvalue ofX at non-vanishing distance from S α asymptotically as n → ∞ and hence the plausible isolated eigenvalues of L α are only due to the matrix UΛU T . We follow classical random matrix approaches used for the study of the spectrum of spiked random matrices [18,28]. From Theorem 1, the eigenvalues λ of L α falling at non-vanishing distance from the limiting support S α solve for large n, 0 = det(L α − λI n ) almost surely for λ / ∈ S α . Since L α −L α a.s. a.s.
We may then just solve 0 = det(D −α q X √ n D −α q + UΛU T − λI n ). Now, as from Proposition 1, the random matrixX does not have eigenvalues at non-vanishing distance from S α asymptotically, for λ / ∈ S α , we can thus factor and cancel out det(X − λI n ) from the previous determinant equation, so that we are left to solve where Q α λ = (X − λI n ) −1 . As we will show next, the matrix I K+1 + U T Q α λ UΛ converges to a deterministic matrix, almost surely for large n. By the argument principle (similar to e.g., [28]), the roots of I K+1 + U T Q α λ UΛ are asymptotically those of the limiting matrix, with same multiplicity and it suffices to study the latter.
We then proceed to retrieving a limit for I K+1 + U T Q α λ UΛ. From Theorem 1, we have The entries (1, 2), (2, 1) and (2, 2) of U T Q α λ U are random as they contain the random matrix X but tend to be deterministic in the limit. In fact, using the resolvent identity, we have that λ D α q 1 n and the entry (2, 2) is equal to n (q T 1n) 2 1 T n X1 n + λ1 T n D 2α q 1 n + λ 2 1 T n D α q Q α λ D α q 1 n . Now, we can freely use Lemma 3 to evaluate the limits of the entries of U T Q α λ U since all the terms are of the form a T Q α λ b with a and b deterministic vectors. From Lemma 3, the entries (1, 1), (1, 2) and (2, 2) converge almost surely respectively to 1 qQ α λ D α q 1 n for large n. Now, using the fact that for any bounded continuous function f , from the law of large numbers, After some algebra, we obtain 1 where the e ij 's are given in Theorem 2. Similarly for the terms (1, 2) and (2, 2), we obtain respectively with v µ = q 2α µ(dq) and where we have also used the fact that 1 n 1 T The limit of I K+1 + U T Q α λ UΛ is then obtained as Using the Schur complement formula for the determinant of block matrices, we have that the determinant of the RHS matrix is zero whenever In particular, when S α = [S α − , S α + ] is composed of a single connected component (as when S α is the support of the semi-circle law as well as most cases met in practice), then isolated eigenvalues of L α may only be found beyond S α , for some non-zero eigenvalue ofM. For the asymptotic spectrum of L α , S α − = −S α + as one can show that for any z ∈ R \ S α , E α 2 (−z) = −E α 2 (z) so that both previous conditions reduce to | |> lim x↓S α . This proves Corollary 1.
where in (27) This proves Remark 4 concerning the eigenvectors of G λ associated to vanishing eigenvalues.

Informative eigenvectors
In this section, we are interested in characterizing the content of the vectorsv α i = (which are used for the classification in the Step 4 of our algorithm) where u α i is the eigenvector of L α corresponding to a unit multiplicity isolated eigenvalue converging to ρ which is such that 1 + θ α (ρ) = 0 (the case 1 + θ α (ρ) = 0 is treated in Appendix Appendix D). By writinḡ where w a i ∈ R n is a random vector orthogonal to j a of norm √ n a , supported on the indices of C a with identically distributed entries, our goal is to estimate the ν a i 's and σ a ij 's the limits of which are asymptotically the limiting class means (ν a EM ) i 's and class covariances (σ a EM ) ij 's that the EM algorithm will find after convergence (whenever it converges to the correct solution) in Step 4 of Algorithm 1.
In Section 3.3, we have shown that the estimation of the ν a i 's requires the evaluation of with G α will produce zero in the limit z → ρ.
where the subscript denotes the first derivative with respect to z. Using the fact that v r,ρ is Next, we need to estimate the denominator term (u α i ) T D 2(α−1) u α i of 1 For u α i an eigenvector of L α associated to an isolated eigenvalue converging to ρ asymptotically, we have As in the previous section, by applying Woodburry idendity, this is equivalent to evaluating where R 3 (z) is a matrix having no residue in the considered contour.
Again here, we just need the top left entry of U T Q α z D 2(α−1) Q α z UΛ which is given from Theorem 1 by e α 00 (z) with e α ij;2 defined in Equation (6). This given, we can show that where χ α (z) is defined in Theorem 4. We thus have By applying l'Hopital rule to evaluate this limit as in the previous section, we obtain . Finally, We recall that one goal of this section is to estimate From Equation (38), the former quantity is easily retrieved and we have This proves Theorem (4).
In Section 3.3, we have shown that to estimate the σ a ij 's, we need to evaluate the more involved object 1 n Similarly to what was done previously for the estimation of 1 where Γ ρ1 and Γ ρ2 are two positively oriented contours circling around some limiting isolated eigenvalues ρ 1 and ρ 2 respectively. We will use the same technique as in the proof of Theorem 4 to evaluate this integrand. Namely, by applying the Woodburry identity to each of the inverse in the integrand, we get 1 2πi 2 where we have used the fact that the cross-terms 1 n J T D α−1 Q α zi D α−1 J, i = 1, 2 will vanish asymptotically as the latter do not have poles in the considered contours. By using the iden- Most of those quantities have been evaluated in the proof of Theorem 4. We thus obtain where R 4 (z 1 , z 2 ) has no poles in the considered contours. It is then sufficient to evaluate the top left entry of each of the matrices The former quantity has already been evaluated in the previous section but for (z, z) replaced by (z 1 , z 2 ) and the diagonal matrix between Q α z1 and Q α z2 being here D α−1 D a D α−1 instead of D 2(α−1) . We thus have Finally, we are left to evaluate We can then perform a residue calculus similar to what was done in the proof of Theorem 4. Additionnaly, we use the fact that the eigenvectors v ρ1 and v ρ2 corresponding to distinct eigenvalues ρ 1 and ρ 2 of the symmetric matrix D(c) −→ e α 00 (ρ i )e α 00 (ρ j ) e α 00;2 (ρ i , ρ i ) + χ α (ρ i ) e α 00;2 (ρ j , ρ j ) + χ α (ρ j ) We are thus now ready to evaluate the σ a ij 's. By definition, The first right hand side term is estimated by dividing 1 = 0 for any couple of indexes (e, f ) such that the aforementioned quantities are non zeros. Indeed from Theorem 4, we have Theorem (5)  We need to prove that ∞ n=1 P (max 1≤i≤n |q i −q i | > η) < ∞ for any η > 0 so that we can conclude from the first Borel Cantelli lemma (Theorem 4.3 in [34]) that P (lim sup n max 1≤i≤n |q i −q i | > η) = 0 from which Lemma 1 unfolds. We have that Let us treat for instance the term P(q i −q i > η) in the following. Since A ij = q i q j +q i q j Mg i g j √ n +X ij with X ij a zero mean random variable, we have Since A, B and C tend to zero in the limit n → ∞, we will next use the fact that P(q i − q i > η) ≤ P(A > η/3) + P(B > η/3) + P(C > η/3) and show that all those individual probabilities vanish asymptotically. Since the term C is deterministic and tends to zero in the limit n → ∞, we have P(C > η/3) = 0 for all large n. Let us then control P(A > η/3) and P(B > η/3). We have with σ 2 = lim sup n max 1≤i≤n q i ( j q j ) − q 2 i ( j q 2 j ) and where in the last inequality of (A.2), we have used Bernstein's inequality (Theorem 3 in [35]) since the A ij 's are independent Bernoulli random variables with variance σ 2 ij = q i q j (1 − q i q j ) + O(n − 1 2 ). For the term B we have where in the inequality (1) we have used the fact that n −1 n j=1 A ij ≤ 1; in the inequality (2) ψ > 0 is any constant smaller than m 2 µ and in the inequality (3) we have used n −2 i,j A ij > √ ψ and the fact that the probability of the intersection between two events is always smaller than the probability of one of those events. It then remains to control B 1 and B 2 . For B 2 we have where the inequality follows from Bernstein's inequality with the similar arguments as previously.
Finally for the term B 1 we have where in the inequality (1)  We then first compute E(XQ). By writingX il = σ il √ n Z il where Z il is a random variable with zero mean and unit variance, we thus have By applying Stein's Lemma (Lemma 5 in Section Appendix A.2), we have where E il is the matrix with all entries equal to 0 but the entry (i, l) which is equal to 1. Using simple algebra, we have since E(y) is equal to 0 in that case. Using Nash Poincaré inequality (Lemma 6 in Section Appendix A.2), one can show that Var(x) = O 1 n 2 [32]. Additionally, ∀i, j and z ∈ C + , |Q ij | ≤ 1 | (z)| . This finally implies that E Q ij tr ΣiQ By using standard techniques [32], one can show that the unique solution e i (z) to e i (z) = Similarly, for any vectors a, b of bounded norms, we may write We also have that i,j (a * Ξ −1 ) i E Q Σi n Q ij b j = O(n −1 ). This can be proved similarly to the proof of Lemma 7. Hence, which is independent of the class information (class proportions or inter-class affinities). This concludes the proof.