A Scalable Community Detection Algorithm for Large Graphs Using Stochastic Block Models

1 Introduction

Community structures, in which nodes belonging to the same community are more densely connected to each other than externally, prevail in real graphs. Finding those structures can be beneficial in many fields, such as finding protein complexes in biological networks and topical or disciplinary groups in collaborative networks [Fortunato and Castellano, 2012; Kemp et al., 2006; Ans´otegui et al., 2012].

In this work, we consider non-overlapping community structures. Many community detection algorithms handle such a problem [Newman, 2004; Raghavan et al., 2007; Blondel et al., 2008; Fortunato and Castellano, 2012; Liu et al., 2013; Wickramaarachchi et al., 2014; Staudt and Meyerhenke, 2015]. However, they come along with limitations for large graphs, for example, in handling community heterogeneity [Nadler and Galun, 2006; Fortunato and Barthelemy, 2007; Xiang and Hu, 2012]. Alternatively, model-based methods can produce more reliable and accurate results when the model assumptions are in accordance with the real graphs. Stochastic block models (SBMs) are among the important probabilistic tools describing the connectivity relationship between pairs of nodes [Holland et al., 1983], and have received considerable attention both in theoretical [Celisse et al., 2012] and application domains [Daudin et al., 2008].

Many algorithms have been proposed to infer the parameters of SBMs, such as Bayesian estimation [Hofman and Wiggins, 2008] and nuclear norm minimization [Chen et al., 2012]. Bayesian estimation defines the prior distributions on latent variables (community labels of nodes) and maximizes the posterior distribution when a graph is given [Hofman and Wiggins, 2008]. Nuclear norm minimization of a matrix minimizes the sum of its singular values, and their algorithm is verified on graphs containing one thousand nodes.

In this work, we propose a multi-stage likelihood maximization algorithm based on SBMs, which has three advantages:

• Speed: We devise an algorithm based on coordinate descent and use approximations to simplify computations. Our algorithm runs in linear time with respect to the number of edges.
• Scalability: We propose a parallel algorithm based on a message-passing, which tends to be the most portable and performance-consistent parallel programming model for a variety of memory structures. We overlap the communication and computation in the algorithm, so that for large graphs, it can achieve significant speedup proportional to the number of processors. To the best of our knowledge, it is the first parallel algorithm for inferring SBMs.
• Quality: The algorithm can produce high-quality results. In the initialization, it considers each node as one community and then employs a multi-stage iterative strategy to construct larger communities gradually. It outperforms traditional community detection algorithms empirically as measured by normalized mutual information (NMI) [Danon et al., 2005] with respect to the

ground-truth. Experiments on real-world graphs show that our algorithm can produce more meaningful communities with good quality.

2 Stochastic BlockModel

In this work, we develop a community detection algorithm based on a stochastic block model (SBM). We define N as the number of nodes and M as the number of undirected and unweighted edges connecting those nodes. Each node belongs to one of K blocks, and we use Z ∈ {0, 1}N×K to represent the block labels. That is, Zir = 1 means node i belongs to block r and each row of Z contains only one nonzero entry. We also define a matrix B ∈ [0, 1]K×K where Brk(r ̸= k) represents the probability of connections between nodes drawn from block r and k, respectively. If r = k, Brk represents the probability of connections inside the block.

Using the matrices B and Z, we define a probability matrix Θ = ZBZ T . Then the adjacency matrix W of a sample network can be drawn from the following model:

for i, j ∈ {1, 2, · · · ,N} and i ̸= j, indicating that Wij is a sample from the Bernoulli distribution with success rate Θij . Typically, the adjacency matrix W is available from the data set. Our primary purpose is to estimate Z.

3 Methodology

For a specific model, the likelihood is a function of the model parameters, describing the probability of obtaining the observed data with these parameters. The maximum likelihood estimator is the setting of parameters that maximizes the likelihood function.

As defined in Eq. (1), if onlyW is given, the log-likelihood function is

It is very time consuming to maximize such a likelihood function directly through traditional optimization methods (for example, branch-and-bound) for large graphs in which there are at least NK unknown variables.

For the sake of speed and scalability, we propose a fast algorithm that updates B and Z in turn to maximize the objective function L(B,Z|W), and use a multi-stage framework to help the solution be close to the global optimum. We also develop a parallel implementation to make the model more scalable.

3.1 Estimation Algorithm

Given W, the maximum likelihood estimates of (B,Z) are defined as

subject to 0 ≤ Bij ≤ 1, Zij ∈ {0, 1}, ∑j Zij = 1. Roughly speaking, we solve the above optimization problem by alternatively updating B and Z and using a community shrinking and expanding strategy to improve accuracy.

We first describe the alternating updating procedure. When Z is fixed and B is considered as unknown, without loss of generality, let β = Brk. If other entries are fixed, L(B,Z|W) = s log(β) + ˆ s log(1 - β) + C, where C is a constant, s = ∑ij Wij(Z:rZ T :k )ij and ˆ s = ∑ij (1 - Wij )(Z:rZ T :k )ij. Taking the derivative of L with respect to β,

and setting the derivative to be zero, we have

As the inter-block connection probabilities are small, we use a representative scalar value to replace the off-diagonal entries of B, which can be computed by counting all the intercommunity edges. Thus, the total time complexity of updating B is O(N) + O(K) + O(M) = O(M).

Theorem 1 For fixed Z, the objective function L(B,Z|W) achieves its global maximum if entries of B are updated according to Eq. (4).

Proof When Z is fixed, because each entry of B can optimize the objective function independently, after updating all the entries by Eq. (4), the resulting B is a stationary point of the objective function and each entry satisfies the constraints.

By taking the second derivative of the objective function, we have

Therefore, when Z is given, the objective function at β (determined by Eq. (4)) is a global maximum. As each entry of B is irrelevant to each other, we can update B by Eq. (4) sequentially.

When B is fixed, we use the block coordinate descent method to update Z row by row. When updating the first Z in ZBZ T in Eq. (2), the algorithm keeps the second Z as its previous estimate Z (t-1) . Then, the likelihood function can be locally maximized by setting all the elements in the row to

0 but Zirmax = 1, where rmax is chosen by

The time complexity for updating one node by Eq. (6) is O(NK) and for all the nodes it is O(N 2K).

To reduce the time complexity, we propose a faster method. Let N (c) be a column vector containing the number of nodes in each community, which is updated once the row Zi: is changed. Let N (d) be a column vector containing the number of nodes connected to node i in each community. Thus, we have a new update rule:

The time complexity is O(Mi)+O(K2 ), where Mi is the degree of node i used in computing N (d) . If we use a scalar to represent the inter-block connection probabilities as before, we can replace the summation in Eq. (7) by two multiplications. We also note that a node will take only a block label from its neighbors, and the time to enumerate over different choices of r becomes O(Mi) for node i in Eq. (7). Thus, the time complexity of computing labels for all the node is reduced to ∑Ni=1 O(Mi) = O(M).

Wij)+(2Wij -1)(Z (2)BZ (1) T )ij ]. Therefore, the objective function can be rewritten as L(B, Z|W) := G(B, Z, Z).

If entries in Z are continuous variables, then the above optimization can find a global maximum.

Theorem 2 min Z(t) -G(B, Z (t) , Z (t-1) ) under the constraints ∑r Zir = 1 and 0 ≤ Zir ≤ 1 is a convex optimization problem when B and Z (t-1) are given.

We omit the proof because it is direct.

However, in our model, Z is a Boolean matrix. Thus, it will probably converge to a local optimum, in which nodes from several true communities may become one single temporary community. Here we refer to the community that a subset of nodes indeed belongs to as the true community for those nodes, and the community determined by an algorithm at each iteration is called the temporary community.

The local optimum problem can be avoided by limiting each temporary community to contain only one true community in the initialization. The next subsection describes the detail.

Community Shrinking and Expanding

A community shrinking and expanding approach works as follows. First, randomly initialize the nodes into αK temporary communities where α should be large enough so that the community size is tiny. Second, run the inference algorithm while gradually merge communities and reduce the number of communities to K. When K is unknown, the algorithm proceeds until no more merging is possible.

This approach reduces the "collision" probability in the initialization significantly, where a "collision" is the situation that two true communities both have most of the nodes in one temporary community. The ratio of the "non-collision" probability after shrinking to the probability without shrinking is:

The lower bound is not sensitive to the choice of α when α is large enough. For example, the lower bounds at α = 10 and α = 100 are 0.9KK K! and 0.99KK K! , respectively, both of which are very significant. In practice, we can choose αK = N.

Initializing an extra number of communities may introduce some tiny communities in the final result. The community expanding is devised by considering the likelihood of the intercommunity edges. When they are more likely to be in certain communities, the corresponding nodes are merged. In detail, for any two communities r and k, the number of edges between them is

which is out of nrk = ∑i Zir × ∑i Zik, the maximum possible connections. Therefore, the log likelihood of the intercommunity edges crk belonging to the true community r can be represented by Lp, while the likelihood of not belonging to the community is Lq.

where p = (crr + crk + ckr + ckk)/(nrr + nrk + nkr + nkk) indicating the internal density after merging, and q = Brk. If Lp > Lq, we merge r and k into one community, otherwise leave them separated.

The time complexity for each merging is O(K2 ). However, if we only consider the pairs of communities that have at least one edge between them, the time complexity becomes O(M). Community merging runs at the end of each stage after updating Z and B. Algorithm 1 is the pseudocode of the serial algorithm.

Convergence Rate

We consider an ideal stochastic block model, which comprises K true communities with the same properties (community size, internal density, etc.). We define p(0) R ≈ 1 K2 as the expected ratio of nodes belonging to the same true community over the total number of nodes over the total number of nodes, in a temporary community at iteration zero.

Theorem 3 For an ideal stochastic block model, ML-SBM with random initialization converges quadratically as

The result can be proved by utilizing the concept of temporary communities.

Algorithm 1 Serial Algorithm (ML-SBM)

set community number to αK initialize B, Z (t) into many tiny communities Z (t-1) = Z (t) repeat repeat for i = 1 : N do compute N (d) for node i update Z (t) i: according to Eq. (7) Z (t-1) i: = Z (t) i: update N (c) end for update B according to Eq. (4) until Z (t) does not change anymore merge communities according to Eq. (8) and (9) until no more merging is possible

3.2 Parallelism

The parallelism is designed as follows. Let SI be the set of integers from 1 to N, partitioned into non-overlapping and approximately equal-size subsets, each of which is represented by SIb for the bth processor. For each processor, we also use local variables sb, ˆ sb and βb, which have the similar definitions as s, ˆ s and β, but taking only the ith row (i ∈ SIb) of W into account.

We choose a message-passing model which is applicable on a variety of memory structures, whether shared, or distributed, or hybrid. We use non-blocking MPI communications to improve the communication-computation overlap. This has two-fold benefits: if the hardware allows, it utilizes the network bandwidth during computation, and it maintains a convergence rate for Z similar to the serial version by transmitting up-to-date community labels. In order to reduce the overall communication data, changes of Zi: are sent to processor b only if i has neighboring nodes in SIb. Because the communication bandwidth will be higher as the message size is larger, we buffer the change information of Z until computing f rows of Z finishes. Here f is chosen to be inversely proportional to the number of MPI ranks, and proportional to N, maintaining a consistent convergence rate.

The algorithm uses similar partitioning and communicating strategies when computing B. Algorithm 2 presents the pseudocode of the parallel algorithm.

In parallel computing, the rows of W are distributed on different processors. If there are totally g processors (equal to the number of MPI ranks) and data are evenly distributed, the computation time of each processor for Eq. (4) and Eq. (7) is O( Mg ). When the computations for B are partitioned on each processor, reducing the number of communities also requires O( Mg ) for each processor. The total message length for Z and B during the iterations is O(g(N + K)). If the bandwidth is constant, the speedup of the parallel algorithm is

where T is a linear function, such that the numerator and de-

Algorithm 2 Parallel Algorithm (Par-SBM)

set community number to αK initialize B, Z (t) into many tiny communities Z (t-1) = Z (t) repeat repeat do in parallel for i ∈ SIb do compute N (d) for node i update Z (t) i: according to Eq. (7) Z (t-1) i: = Z (t) i: synchronize changes of Zi: at frequency f end for update sb and ˆ sb for each βb entries synchronize s = ∑b sb and ˆ s = ∑b ˆ sb update B according to Eq. (4) update N (c) until Z (t) does not change anymore merge communities according to Eq. (8) and (9) until no more merging is possible

nominator indicate the running times of the serial version and parallel version, respectively. Therefore, the speedup will increase towards g if the number of edges increases given the same N and K.

4 Experimental Results

4.1 The Serial Algorithm

In this section, we compare our algorithm (ML) with other serial algorithms in MATLAB. The competitors include the Louvain method (LV) [Blondel et al., 2008], a Bayesian inference algorithm (BI) [Hofman and Wiggins, 2008], spectral clustering (SC) [Von Luxburg, 2007], label propagation (LP) [Raghavan et al., 2007], and modularity optimization (MM) [Newman, 2004]. We use the LFR benchmark generator [Lancichinetti and Fortunato, 2009] to create graphs and run all of them for comparison. Each generated graph is of size 1000, average degree 30, maximum degree 50, exponent of degree distribution -2, exponent of community size distribution -2, minimum community size 20, and maximum size 100. µ is the proportion of inter-community connections to the intra-community ones. The accuracy is evaluated by normalized mutual information (NMI) [Danon et al., 2005], which compares the similarity between the computed community structure and the ground-truth one. The larger the NMI is, the more similar the two schemes are. If two schemes are identical, NMI is 1. The results are averaged on five runs.

From Figs. 1(a) and 1(b) we can see that our algorithm is relatively fast and has the best accuracy. When µ = 0, most of the algorithms can find the correct result, and our algorithm is the second fastest one. As µ increases to 0.7, our algorithm needs a longer time to iterate, but its accuracy is still the best. Only SC and LV can achieve a similar accuracy over a similar amount of time.

However, for large graphs, SC is not scalable, and LV is not accurate. Our algorithm is the method of choice for both scalability and accuracy. Figs. 1(c) and 1(d) show the results on graphs similar as the previous experiment but the average degree is 10 and µ = 0.1. On large graphs, compared with other algorithms, ours consistently produces more accurate results with respect to the ground-truth setting, and the running time growth rate is smaller than others as the graph size grows. The slow growth of the running time of our algorithm is in accordance with the theoretical time complexity which is linear in the number of edges.

4.2 The Parallel Algorithm

In this section, we exploit distributed memory parallelism for larger graphs. When the input graph has more than one million nodes, typical algorithms for stochastic block models are not able to infer the parameters within a reasonable amount of time (for example, 24 hours). Therefore, we compare our algorithm (PS+number of processors) to a popular community detection algorithm, the Louvain method (LV) [Blondel et al., 2008] and a fast graph partitioning algorithm Metis (MT) [Karypis and Kumar, 1998]. The Louvain method implemented in C can also generate level one structures as byproducts containing the finest communities (LVF). Those communities are merged into larger ones in LV. Parallel algorithms are often extended from serial versions [Wickramaarachchi et al., 2014], but the results are typically invariant. We use several metrics to compare the results, such as NMI, modularity [Clauset et al., 2004], and conductance (the smaller, the better) [Leskovec et al., 2009]. The results on those quantities are averaged on five runs.

First, we run our algorithm on benchmark graphs (LFR1e6) generated by LFR benchmark [Lancichinetti and Fortunato, 2009]. The parameters are similar to the serial case except the graph size is increased to one million, and we change two parameters (the mixing parameter µ and maximum community size mc) to have more variations. Fig. 2 shows the results, in which GT is the ground truth given by the generator. Generally, the increase of µ and mc will decrease the community detection accuracy, but the results generated by our algorithm are always most close to the ground truth according to NMI. The results with our algorithm are also very close to the ground truth in numbers of communities (except for MT which is predefined) and modularity, and are usually better than others in modularity. In Fig. 2(b), as the problem becomes more difficult from left to right, our algorithm spends more time to maintain the quality, while the other algorithms finish within the same amount of time while sacrificing the quality. In addition, the performance of our algorithm is stable as the running results using 8 processors (P8) and 64 processors (P64) are fairly close.

We also compare algorithms on real-world graphs: citPatents [Leskovec et al., 2005] and dblp-Coauthor [Ley, 2002]. Nodes in cit-Patents are patents in U.S. patent dataset and edges represent the citation relationships between the two patents. dblp-Coauthor is a coauthor graph extracted from publications in DBLP database. It is also worth noting that LVF has higher NMI and more communities than LV in benchmark test. It indicates that LV suffers from the resolution limit problem [Fortunato and Barthelemy, 2007] that prefers unrealistically large communities. A similar problem happens on LV for real-world data as it generates many communities containing more than 104 nodes (Fig. 3). Especially, for example, for dblp-Coauthor network, the biggest community by LV contains about 105 nodes (10% of the whole graph). It is not realistic for collaboration networks and is orders of magnitude larger than a reasonable community size for human interaction [Leskovec et al., 2009], so LV is omitted from the remaining comparisons. Alternatively, choosing low level structures such as LVF can be a remedy [Good et al., 2010], although it still contains a few extremely large communities. On the other hand, MT tends to find even-size communities, which is not desirable for community detection either. Table 2 shows the quality of results by different algorithms.

4.3 Microscopic Example on dblp-Coauthor Network

In this section, we analyze the community detection results of a run by our Par-SBM (PS) and the Louvain method [Blondel et al., 2008] using random initialization. For the Louvain method, we pick the Level One structure with finest communities (LVF) and the Level Two structure denoted by LV2. Communities become larger in higher level structures by the Louvain method. Taking the Conference Chair as an example, we consider community containing author "Michael Wooldridge," and use SPS, SLV F , and SLV 2 to represent the nodes in the community found by Par-SBM and the Louvain method, respectively. The community size by our algorithm is comparable to the one by LVF, but much smaller than the one by LV2, as |SPS| = 255 and |SLV F | = 166, while |SLV 2| = 26294.

All the algorithms are able to identify blue nodes in Fig. 4 as a separate community, while our Par-SBM provides the clearest detail. SLV 2 contains one hundred times more nodes than SPS. Among all the communities found by Par-SBM, the red, green, and purple communities in Fig. 4(a) are the top-three contributors to SLV 2, owning 2.96% of total nodes in SLV 2. The three communities can be represented by three nodes respectively: purple by Dr. Michael Wooldridge (multi-agent systems), green by Dr. Bruce W. Porter (knowledge systems), and red by Dr. Christopher W. Geib (Reports of the AAAI Conference Workshops). Nodes in red forms a strong community because the authors (although from different areas of AI) are fully connected by the workshop reports annually. However, some nodes in the aforementioned three communities are inappropriately separated by LVF as shown in Fig. 4(b).

5 Concluding Remarks

In this paper, we have proposed a fast algorithm for clustering nodes in large graphs using stochastic block models. We have adopted alternative updates and coordinate descent methods to infer the latent parameters, and used community shrinking and expanding to improve accuracy. Our algorithm has linear time complexity in the number of edges, and can scale to multi-processor systems.

Compared with other community detection algorithms, our SBM-based method can produce high quality results on benchmark graphs and find interesting communities in the real-world graphs.

This work can boost the application of SBMs on big data and bring new insights by overcoming the accuracy or running time shortages of traditional algorithms. The code is available at https://github.com/pdacorn/par-sbm.