A scalable community detection algorithm for large graphs using stochastic block models

Abstract

Community detection in graphs is widely used in social and biological networks, and the stochastic block model is a powerful probabilistic tool for describing graphs with community structures. However, in the era of “big data”, traditional inference algorithms for such a model are increasingly limited due to their high time complexity and poor scalability. In this paper, we propose a multi-stage maximum likelihood approach to recover the latent parameters of the stochastic block model, in time linear with respect to the number of edges. We also propose a parallel algorithm based on message passing. Our algorithm can overlap communication and computation, providing speedup without compromising accuracy as the number of processors grows. For example, to process a real-world graph with about 1.3 million nodes and 10 million edges, our algorithm requires about 6 seconds on 64 cores of a contemporary commodity Linux cluster. Experiments demonstrate that the algorithm can produce high quality results on both benchmark and real-world graphs. An example of finding more meaningful communities is illustrated consequently in comparison with a popular modularity maximization algorithm.

Keywords

Stochastic block model parallel computing community detection

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, disciplinary groups in collaborative networks, etc. [35, 47, 17, 51, 26, 2, 46].

In this work, we consider non-overlapping community structures. Many community detection algorithms handle such a problem [37, 40, 6, 17, 34, 49, 43, 32]. However, they come along with limitations for large graphs, for example, in handling community heterogeneity [36, 16, 50].

Alternatively, model-based methods can produce more reliable and accurate results when the model assumptions are in accordance with real graphs. Stochastic block models (SBMs) are among the important probabilistic tools describing the connectivity relationship between pairs of nodes [23], and have received considerable attention both in theoretical [7] and application domains [14].

Many algorithms have been proposed to infer the parameters of SBMs, such as Bayesian estimation [22] and nuclear norm minimization [8]. Bayesian estimation defines the prior distributions on latent variables (community labels of nodes) and maximizes the posterior distribution when a graph is given [22]. Nuclear norm minimization solves a relaxed problem to find the community indication matrix.

We can also utilize parallel computers to detect communities on large graphs. Parallel computing has been proposed for decades [19, 3]. It can decrease the execution time of applications by executing them on multiple processors [15].

Parallel computers are classified according to the programming model natural to their memory structure, whether shared, or distributed, or hybrid. In a shared memory model [12], each processor can access the entire system memory, while each processor in a distributed memory model [21] can access only its local memory and requires an interprocessor network to communicate data to and from other local memories. In the hybrid model [39], subsets of processors share memory and these subsets are connected by a communication network. The hardware supporting these programming models need not fully resemble them. The convenient global addressing of shared memory can be implemented over distributed memory hardware and the efficient controlled access of distributed memory over shared memory hardware. However, there is a general association of distributed memory protocols with the systems of the largest memory capacities. Therefore, in our algorithm, we choose the message passing model as the most performance-portable to the hardware resources that will be required for graph applications of the future.

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 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) [13] with respect to the ground-truth. Experiments on real-world graphs show that our algorithm can produce more meaningful communities with good quality.

2. Related work

Many community detection algorithms have been proposed in the past decades. Some algorithms rely on a quality function, for example, modularity. Modularity is the fraction of inter-community connections minus the expected fraction when edges are randomly placed [38]. Newman [37] optimizes the modularity in a hierarchical manner. A typical spectral clustering approach [45] can also detect communities. Li and Schuurmans [33] use an iterative rounding strategy for the spectral method to improve the quality of the results. The Louvain method [6] improves the modularity hierarchically. Nevertheless, they may have some disadvantages. For example, modularity-based algorithms may have resolution limit [16], which means they may be unable to find small communities without modifying their objective function [44]. In some cases, other algorithms can give high quality results when maximizing the modularity, but the time complexity could be high. For example, spectral clustering has a typical time complexity no less than $O(NK^{2})$ [48], where ${N}$ is the number of nodes and ${K}$ is the number of communities. Zhang and Newman [52] developed a multiway spectral algorithm that can divide graphs into many communities and achieves good quality in terms of modularity, but the time complexity is still high. Liu et al. adopted spectral clustering on a reduced graph, and the run time is quadratic to the size of the reduced graph [34]. Block modularity maximization [5] although reduces the run time from cubic to quadratic in terms of the number of clusters, is still expensive and not suitable for large graphs with many communities.

On the other hand, using SBMs is more reliable under particular assumptions. Theoretically, the threshold of exact community recovery can been identified [1]. In practice, many algorithms have been proposed to infer the parameters of SBMs. Bayesian estimation defines the prior distributions on latent variables (community labels of nodes) and maximizes the posterior distribution when a graph is given [22]. Nuclear norm minimization of a matrix minimizes the sum of its singular values, and their algorithm is verified on graphs containing one thousand nodes [8]. Spectral algorithms also demonstrate the capability of inferencing parameters of SBMs [9]. However due to the high time complexity, the run time of those algorithms typically becomes enormous when the graph size grows to millions or billions, which prohibits the application on big graphs.

Recent years have seen work on parallel implementations for non-SBM community detection algorithm. For example, Riedy et al. proposed an parallel agglomerative approach for community detection that supports both maximizing modularity and minimizing conductance [41]. Bhowmick and Srinivasan used a shared memory interface for parallelizing the Louvain method [4]. Staudt and Meyerhenke developed a parallel framework for multiple algorithms, including label propagation and the Louvain method with refinement [43]. Gregori et al. devised an parallel algorithm for k-clique community detection [20]. Parallel algorithms for spectral clustering are also proposed [24]. Those implementations are mostly on shared-memory systems, which constrain the maximum speedup due to the hardware limit. In addition, none of these implementations can suitably address the problem in SBMs.

In this work, we propose a community shrinking and expanding approach to overcome the limitation of the popular coordinate descent method for stochastic blockmodels. It is also scalable on parallel computers. Our fast and scalable algorithm fills the gap of processing very large graphs using stochastic block models.

3. Stochastic block model

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\in\{0,1\}^{N\times K}$ to represent the block labels. That is, $Z_{ir}=1$ means node $i$ belongs to block $r$ and each row of $Z$ contains only one nonzero entry. We also define a matrix $B\in[0,1]^{K\times K}$ where $B_{rk}(r\neq k)$ represents the probability of connections between nodes drawn from block $r$ and $k$ , respectively. If $r=k$ , $B_{rk}$ represents the probability of connections inside the block.

Using the matrices $B$ and $Z$ , we define a probability matrix $\Theta=\textit{ZBZ}^{T}$ . Then the adjacency matrix $W$ of a sample network can be drawn from the following model:

$\displaystyle\Pr(W_{ij})=\begin{cases}\Theta_{ij}&\text{ if }\;W_{ij}=1,\\ 1-\Theta_{ij}&\text{ if }\;W_{ij}=0,\\ \end{cases}$ (1)

for $i,j\in\{1,2,\cdots,N\}$ and $i\neq j$ , indicating that $W_{ij}$ is a sample from the Bernoulli distribution with success rate $\Theta_{ij}$ .

For example, we can define a probability matrix

$\displaystyle B=\begin{bmatrix}0.2&0.01\\ 0.01&0.2\end{bmatrix}_{2\times 2},$

meaning that nodes inside each community are connected at probability 0.2 and nodes from different communities are connected at probability 0.01. We also define a matrix

$\displaystyle Z=\begin{bmatrix}1&&0\\ 1&&0\\ \vdots\\ 1&&0\\ 0&&1\\ 0&&1\\ \vdots\\ 0&&1\end{bmatrix}_{100\times 2},$ (2)

indicating the first half and the second half of the nodes are from two different groups. Thus, we can have the probability matrix

$\displaystyle\Theta=\begin{bmatrix}0.2&\ldots&0.2&0.01&\ldots&0.01\\ \vdots&&\vdots&\vdots&&\vdots\\ 0.2&\ldots&0.2&0.01&\ldots&0.01\\ 0.01&\ldots&0.01&0.2&\ldots&0.2\\ \vdots&&\vdots&\vdots&&\vdots\\ 0.01&\ldots&0.01&0.2&\ldots&0.2\end{bmatrix}_{100\times 100}.$ (3)

In this case, a possible adjacency matrix $W$ can be generated from $\Theta$ , by considering each edge indicator $W_{ij}$ (for the edge between node $i$ and $j$ ) as a binomial random variable with parameter $\Theta_{ij}$ . We can use an image to represent the matrix $W$ , in which non-zero entries (corresponding to edges) are plotted in blue, as shown in Fig. 1. The corresponding matrix $W$ here is of size $100\times 100$ with 1042 non-zero entries.

Figure 1.

One adjacency matrix generated by the model.

Typically in practice, the adjacency matrix $W$ is available from the data set. Our primary purpose is to estimate $Z$ .

4. Model analysis: An ideal example

In this section, through a simplified example, we illustrate how our algorithm solves the SBM problem. In particular, we assume the diagonal part of the ground truth $B$ is $p$ , and the off-diagonal part is $q$ . Thus, each pair of nodes from one community is connected with probability $p$ , and those from two different communities are connected with probability $q$ . As mentioned earlier, $N$ is the number of nodes in the graph, and $K$ is the number of communities. We define $n$ as the average number of nodes in each community.

4.1 Self-assembling process

We define 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 as the temporary community.

We initialize $Z$ by putting nodes almost evenly and randomly into different temporary communities. Thus, we can define $p_{R}^{(0)}\approx\frac{1}{K^{2}}$ as the expected fraction of nodes belonging to the same true community over the total number of nodes, in a temporary community at time 0. Due to the randomness, there is very likely a small bias that some temporary community should have a larger number of nodes from the same true community.

At each iteration, each node will associate with a particular temporary community where it has a tight connection. In this case, the temporary community with more nodes from the same true community will take an advantage because the probability of the intra-community connection is higher than other temporary communities. Due to this bias, the association will result in

$\displaystyle p_{R}^{(t)}=\frac{p_{R}^{(t-1)}}{\frac{1}{K}}\times\left({\frac{% 1}{K}-p_{R}^{(t-1)}}\right)+p_{R}^{(t-1)}=2p_{R}^{(t-1)}-(p_{R}^{(t-1)})^{2}\times K$ (4)

for the fraction of nodes in a temporary community and from the same true community. In Eq. (4), $\frac{1}{K}-p_{R}^{(t-1)}$ is the fraction of nodes that has not yet been correctly clustered, and $\frac{p_{R}^{(t-1)}}{\frac{1}{K}}$ is the fraction of wrongly clustered nodes that can be fixed in the $t$ th iteration. As the maximum $p_{R}^{(t)}$ is $\frac{1}{K}$ , the condition $2>p_{R}^{(t-1)}K$ is always true. Therefore, this equation indicates at the early stage, nodes from the same true community in other temporary communities will tend to assemble to a one temporary community.

In practice, we generate many more temporary communities than the number of true communities at the initialization, to avoid two or more true communities assembling at the same temporary community.

4.2 Minimum gap between

p

and

q

The gap between $p$ and $q$ represents the difference in edge densities within and across communities. A smaller value of $p-q$ in SBMs causes more challenges in the identification of the true $Z$ . Thus, a parametric inference algorithm is more powerful if it allows a smaller gap between $p$ and $q$ . In this section, we try to find the lower bound of $p-q$ .

Let $\Pr(Z_{ik}=1)$ be the probability that node $i$ of true community $k$ can be correctly assigned via maximum likelihood. We define the destination temporary community as the temporary community that node $i$ will be assigned to after the current iteration. When correctly assigned, the destination temporary community should have more nodes connected to $i$ than any other temporary community. We define $X$ and $Y$ as collections of binary random variables, and $X_{j}=1$ means node $i$ has a connection to node $j$ in set $X$ . We use $\sum_{j}X_{j}$ and $\sum_{j}Y_{j}$ to represent the number of nodes connected to $i$ in the destination temporary community and one other temporary community respectively.

Then we will demonstrate that node $i$ can be correctly identified at least in certain probability $p_{I}$ , where $p_{I}$ is the probability for correct identification between the true community and all the other communities. Because there are $K-1$ communities other than the true one, then, the probability of correct identification from the true community and one of the other communities is $p_{I}^{\frac{1}{K-1}}$ .

Using an auxiliary variable $d$ to represent the cutting point of $\sum_{j}X_{j}$ and $\sum_{j}Y_{j}$ , we have

$\displaystyle p_{I}^{\frac{1}{K-1}}\geqslant\Pr(\sum_{j}X_{j}>d)\times\Pr(\sum% _{j}Y_{j}<d)=\big{[}1-\Pr(\sum_{j}X_{j}<d)\big{]}\big{[}1-\Pr(\sum_{j}Y_{j}>d)% \big{]}=\Big{[}1-\Pr(\sum_{j}X_{j}<pn-(pn-d))\Big{]}\times\Big{[}1-\Pr(\sum_{j% }Y_{j}>qn+(-qn+d))\Big{]}\geqslant\Big{[}1-e^{-\frac{(pn-d)^{2}}{2pn}}\Big{]}% \Big{[}1-e^{-\frac{2(-qn+d)^{2}}{n}}\Big{]}=\gamma^{\frac{1}{K-1}},$ (5)

where $n$ is the expected number of nodes within the community and $\gamma$ is a lower bound of probability $p_{I}$ . The last inequality is obtained from Chernoff’s inequality [10]

$\displaystyle Pr\left({\sum_{j}X_{j}<pn+\delta}\right)\leqslant e^{-\frac{% \delta^{2}}{2pn}}$

and Hoeffding’s inequality

$\displaystyle Pr\left({\sum_{j}Y_{j}>qn+\delta}\right)\leqslant e^{-\frac{2% \delta^{2}}{n}}$

for Bernoulli random variables.

By choosing an appropriate $d$ , the lower bound of $p_{I}^{\frac{1}{K-1}}$ can be well approximated. If we choose $d$ to have the two factors in Eq. (5) equal, we obtain the following expressions

$\displaystyle pn-d=2\sqrt{p}(d-qn),$ $\displaystyle d=\frac{pn+2\sqrt{p}(qn)}{1+2\sqrt{p}},$ $\displaystyle pn-d=pn-\frac{pn+2\sqrt{p}(qn)}{1+2\sqrt{p}}=\frac{pn+2\sqrt{p}(% pn)-pn-2\sqrt{p}(qn)}{1+2\sqrt{p}}=\frac{2\sqrt{p}(pn-qn)}{1+2\sqrt{p}}.$

Hence,

$\displaystyle\frac{(pn-d)^{2}}{2pn}=\frac{\frac{4p(pn-qn)^{2}}{(1+2\sqrt{p})^{% 2}}}{2pn}=\frac{2n(p-q)^{2}}{(1+2\sqrt{p})^{2}}.$ (6)

Because two factors in Eq. (5) are the same, we have

$\displaystyle p_{I}^{\frac{1}{K-1}}\geqslant\Big{[}1-e^{-\frac{2n(p-q)^{2}}{(1% +2\sqrt{p})^{2}}}\Big{]}^{2}$ (7) $\displaystyle=\gamma^{\frac{1}{K-1}}.$ (8)

As the lower bound of the probability $p_{I}$ , $\gamma$ is expected to be as large as possible. Then we can have $1-\gamma\leqslant 1-c_{0}N^{-10}$ for some constant $c_{0}$ to meet the definition of high probability, where $N$ is the total number of nodes. According to Eq. (8), we have

$\displaystyle\frac{2n(p-q)^{2}}{(1+2\sqrt{p})^{2}}=-\log(1-\gamma^{\frac{1}{2(% K-1)}})\geqslant-\log[1-(c_{0}N^{-10})^{\frac{1}{2(K-1)}}]=(c_{0}N^{-10})^{% \frac{1}{2(K-1)}}+\frac{1}{2}(c_{0}N^{-10})^{\frac{2}{2(K-1)}}+\frac{1}{3}(c_{% 0}N^{-10})^{\frac{3}{2(K-1)}}+\cdots=\Omega(N^{-\frac{5}{K-1}}),$

where the expansion is based on the Taylor series.

Therefore, as $K$ is no less than two and $n\approx\frac{N}{K}$ , we can have the following approximation

$\displaystyle p-q=\sqrt{(1+2\sqrt{p})^{2}\Omega(N^{[-\frac{5}{K-1}-1]}K)}% \approx\Omega\Big{(}\sqrt{\frac{pK}{N}}\Big{)},$ (9)

Thus, the lower bound of $p-q$ by our algorithm can be as small as $\Omega\Big{(}\sqrt{\frac{pK}{N}}\Big{)}$ . Compared to many other methods summarized in [8], the lower bound by our approach is relatively small.

4.3 Graphs with heterogeneous partitions

It is possible that the sizes of the partitions are different, and we analyze the consequent self-assembling process under this assumptions.

First, we define $N_{1},N_{2},...,N_{K-1},N_{K}$ in ascending order as the size of the $K$ communities. Without loss of generality, we assume the $k$ th partition corresponding to $N_{k}$ will emerge in the $k$ th temporary community.

Then, when temporary communities are initialized to have the same size, nodes from different communities are assigned to temporary communities in equal probability. To have the self-assembling process still work, we shall have the following expression, so that nodes from community 1 will not assemble to temporary community $K$ .

$\displaystyle\frac{qN_{K}}{K}<\frac{pN_{1}}{K},$ (10)

where $k\in\{1,2,...,K\}$ . Under the above condition, nodes from other communities are expected to assemble to correct temporary communities because of the following inequality.

$\displaystyle\frac{qN_{k}}{K}<\frac{qN_{K}}{K}<\frac{pN_{1}}{K}<\frac{pN_{k}}{% K}.$ (11)

During the self-assembling procedure, similar to Eq. (4), we have the following expression for community $k$

$\displaystyle p_{R}^{(t)}=\frac{p_{R}^{(t-1)}}{\frac{N_{k}}{N}}\times\left({% \frac{N_{k}}{N}-p_{R}^{(t-1)}}\right)+p_{R}^{(t-1)}=2p_{R}^{(t-1)}-(p_{R}^{(t-% 1)})^{2}\times\frac{N}{N_{k}^{(t-1)}}.$

Hence, the convergence rate is still quadratic.

At the end of the iteration, to obtain the correct result, the following expression is expected

$\displaystyle qN_{K}<pN_{1},$ (12)

which is in accordance with the assembling condition at the initialization in Eq. (10).

For undirected networks, $B$ is symmetric. When the probability matrix contains a variety of values, rather than only $p$ and $q$ , the following inequality is required as a necessary condition for converging to the correct result.

$\displaystyle B_{ij}N_{i}<B_{jj}N_{j},\;\text{ for }\;i>j.$ (13)

5. 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 only $W$ is given, the log-likelihood function is

$\displaystyle L(B,Z|W)=\sum_{i\neq j}\log\Pr(W_{ij})=\sum_{i\neq j}\log[(1-W_{% ij})+(2W_{ij}-1)\Theta_{ij}].$

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 $N K$ 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.

5.1 Estimation algorithm

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

$\displaystyle\mathop{\text{argmax}}\limits_{B,Z}\;\Big{\{}L(B,Z|W)$ $\displaystyle=\sum\limits_{i\neq j}\log[(1-W_{ij})+(2W_{ij}-1)\Theta_{ij}]$ $\displaystyle=\sum\limits_{i\neq j}\log[(1-W_{ij})+(2W_{ij}-1)(\textit{ZBZ}^{T% })_{ij}]\Big{\}},$ (14)

subject to $0\leqslant B_{ij}\leqslant 1$ , $Z_{ij}\in\{0,1\}$ , $\sum_{j}Z_{ij}=1$ . 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 $\beta=B_{rk}$ . If other entries are fixed, $L(B,Z|W)=s\log(\beta)+\hat{s}\log(1-\beta)+C$ , where $C$ is a constant, $s=\sum_{ij}W_{ij}(Z_{:r}Z_{:k}^{T})_{ij}$ and $\hat{s}=\sum_{ij}(1-W_{ij})(Z_{:r}Z_{:k}^{T})_{ij}$ .

Taking the derivative of $L$ with respect to $\beta$ ,

$\displaystyle\frac{\partial L}{\partial\beta}=\frac{s}{\beta}-\frac{\hat{s}}{1% -\beta}$ (15)

and setting the derivative to be zero, we have

$\displaystyle\beta=\frac{s}{s+\hat{s}}.$ (16)

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 inter-community edges.

Thus, the total time complexity of updating $B$ is ${\cal O}(N)+{\cal O}(K)+{\cal O}(M)={\cal 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. (16).

Proof..

When $Z$ is fixed, because each entry of $B$ can optimize the objective function independently, after updating all the entries by Eq. (16), 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

$\displaystyle\frac{\partial^{2}L}{\partial\beta^{2}}\bigg{|}_{\beta=\frac{s}{s% +\hat{s}}}=-\frac{s}{\beta^{2}}-\frac{\hat{s}}{(1-\beta)^{2}}\bigg{|}_{\beta=% \frac{s}{s+\hat{s}}}<0.$ (17)

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

When $B$ is fixed, we use the block coordinate descent method to update $Z$ row by row. When updating the first $Z$ in $\textit{ZBZ}^{T}$ in Eq. (14), 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 $Z_{ir_{\max}}=1$ , where $r_{\max}$ is chosen by

$\displaystyle r_{\text{max}}=\mathop{\text{argmax}}\limits_{r}\;\sum\limits_{j% \neq i}\log\big{[}(1-W_{ij})+(2W_{ij}-1)(B[Z^{(t-1)}]^{T})_{rj}\big{]}.$ (18)

The time complexity for updating one node by Eq. (18) is ${\cal O}(NK)$ and for all the nodes it is ${\cal O}(N^{2}K)$ .

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 $Z_{i:}$ 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:

$\displaystyle r_{\max}=\mathop{\text{argmax}}\limits_{r}\;\sum\limits_{k=1}^{K% }\big{[}N^{(d)}_{k}\log(B_{rk})+(N^{(c)}_{k}-N^{(d)}_{k})\log(1-B_{rk})\big{]}.$ (19)

The time complexity is ${\cal O}(M_{i})+{\cal O}(K^{2})$ , where $M_{i}$ 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. (19) 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 ${\cal O}(M_{i})$ for node $i$ in Eq. (19). Thus, the time complexity of computing labels for all the node is reduced to $\sum_{i=1}^{N}{\cal O}(M_{i})={\cal O}(M)$ .

For fixed $B$ , we define $G(B,Z^{(2)},Z^{(1)})=\sum_{i\neq j}\log[(1-W_{ij})+(2W_{ij}-1)(Z^{(2)}B{Z^{(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. Optimization $\min\limits_{Z^{(t)}}\;-G(B,Z^{(t)},Z^{(t-1)})$ under the constraints $\sum_{r}Z_{ir}=1$ and $0\leqslant Z_{ir}\leqslant 1$ is a convex optimization problem when $B$ and $Z^{(t-1)}$ are given.

Proof..

First we reshape the matrix $Z$ into a column vector $z$ where $Z_{ij}$ corresponds to $z_{K(i-1)+j}$ .

Let $H\in\mathbb{R}^{NK\times NK}$ be the Hessian matrix of the objective function $G(B,Z^{(t)},Z^{(t-1)})$ with respect to $z^{(t)}$ . Let $M=\log[(1-W)+(2W-1){\odot}Z^{(t)}B{Z^{(t-1)}}^{T}]$ , which is the matrix in $G(B,Z^{(t)},Z^{(t-1)})$ before the final summation, and we only need to consider the off-diagonal part of the matrix in this prove. Here ${\odot}$ denotes the element-wise product of two matrices.

If $[Z^{(t-1)}_{i:}]^{T}$ are constants, the second derivative of the $(i,j)$ th term ( $i\neq j$ ) in $M$ with respect to variables other than $[Z^{(t)}_{i:}]^{T}$ is 0, while with respect to the variables in $[Z^{(t)}_{i:}]^{T}$ , it is as follows

$\displaystyle\frac{\partial^{2}-M_{ij}}{(\partial{[Z^{(t)}_{i:}]^{T}})^{2}}=% \frac{\partial^{2}-\log[(1-W_{ij})+(2W_{ij}-1)Z^{(t)}_{i:}B[Z^{(t-1)}_{j:}]^{T% }]}{(\partial{[Z^{(t)}_{i:}]^{T}})^{2}}=\frac{\partial-\frac{(2W_{ij}-1)B[Z^{(% t-1)}_{j:}]^{T}}{(1-W_{ij})+(2W_{ij}-1)Z^{(t)}_{i:}B[Z^{(t-1)}_{j:}]^{T}}}{% \partial[Z^{(t)}_{i:}]^{T}}=\frac{[(2W_{ij}-1)]^{2}B[Z^{(t-1)}_{j:}]^{T}[Z^{(t% -1)}_{j:}]B^{T}}{[(1-W_{ij})+(2W_{ij}-1)Z^{(t)}_{i:}B[Z^{(t-1)}_{j:}]^{T}]^{2}},$

where the square operation is element-wise. It is positive semi-definite because it is a covariance matrix multiplied with a positive scalar.

Therefore, for any column vector $x$ , we have

$\displaystyle-x^{T}Hx=-x^{T}\frac{\partial^{2}\sum_{i\neq j}M_{ij}}{\partial{z% }^{2}}x=\sum_{i\neq j}x^{T}\frac{\partial^{2}-M_{ij}}{\partial{z}^{2}}x% \geqslant 0.$ (20)

Hence, the objective function is convex. As obviously the constraints are convex as well, this problem is a convex optimization. ∎

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. 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.

5.1.1 Community shrinking and expanding

A community shrinking and expanding approach works as follows. First, randomly initialize the nodes into $\alpha K$ temporary communities where $\alpha$ 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 fraction of the “non-collision” probability after shrinking to the probability without shrinking is:

$\displaystyle\prod_{k=0}^{K-1}\frac{1-\frac{k}{\alpha K}}{1-\frac{k}{K}}=\prod% _{k=0}^{K-1}\left[{1+\frac{(1-\frac{1}{\alpha})k}{K-k}}\right]\geqslant\left({% 1-\frac{1}{\alpha}}\right)\frac{K^{K}}{K!}.$

The lower bound is not sensitive to the choice of $\alpha$ when $\alpha$ is large enough. For example, the lower bounds at $\alpha=10$ and $\alpha=100$ are $\frac{0.9K^{K}}{K!}$ and $\frac{0.99K^{K}}{K!}$ , respectively, both of which are very significant. In practice, we can choose $\alpha 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 inter-community 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

$\displaystyle c_{rk}=\sum_{i\neq j}Z_{ir}W_{ij}Z_{jk},$

which is out of $n_{rk}=\sum_{i}Z_{ir}\times\sum_{i}Z_{ik}$ , the maximum possible connections. Therefore, the log likelihood of the inter-community edges $c_{rk}$ belonging to the true community $r$ can be represented by $L_{p}$ , while the likelihood of not belonging to the community is $L_{q}$ .

$\displaystyle L_{p}=c_{rk}\log(p)+(n_{rk}-c_{rk})\log(1-p),$ (21) $\displaystyle L_{q}=c_{rk}\log(q)+(n_{rk}-c_{rk})\log(1-q),$ (22)

where $p=(c_{rr}+c_{rk}+c_{kr}+c_{kk})/(n_{rr}+n_{rk}+n_{kr}+n_{kk})$ indicating the internal density after merging, and $q=B_{rk}$ . If $L_{p}>L_{q}$ , we merge $r$ and $k$ into one community, otherwise leave them separated.

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

Serial Algorithm (ML-SBM)set community number to $\alpha K$ initialize $B$ , $Z^{(t)}$ into many tiny communities

$Z^{(t-1)}=Z^{(t)}$

$i=1:N$ compute $N^{(d)}$ for node $i$ update $Z_{i:}^{(t)}$ according to Eq. (19) $Z_{i:}^{(t-1)}=Z_{i:}^{(t)}$ update $N^{(c)}$ update $B$ according to Eq. (16) $Z^{(t)}$ does not change anymore merge communities according to Eqs (21) and (22)

no more merging is possible

In Fig. 2, we use the graph in Fig. 1 as an example to illustrate the process. When two nodes at both ends of an edge belong to the same community, we assign the edge with a specific color corresponding to the community. Otherwise, the edge is in background color (white) and invisible.

Figure 2.

Changes of community labels in each iteration.

At the beginning, we shrink the temporary communities as initialization, so at Iteration 1 of Stage 1, there are many small communities (temporary communities) and the edge colors are very diverse. At the end of Stage 1, due to the self assembling process, some temporary communities become a bit larger. After each stage, we expand the communities by merging temporary communities that very likely belong to the same community. Therefore, after the transition of each stage, we can see much larger communities which are slightly refined in the stage. Eventually, the algorithm converges and the two true communities are identified.

5.1.2 Convergence rate

In this part, we consider an ideal stochastic block model. We assume that there are $K$ true communities, each containing the same number of nodes. We initialize $Z$ by putting nodes almost evenly and randomly into different temporary communities, with naturally a small bias that some temporary community should have a larger number of nodes from the same true community. Namely, we can define $p_{R}^{(0)}\approx\frac{1}{K^{2}}$ as the expected fraction of nodes belonging to the same true community 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

$\displaystyle\lim\limits_{t\rightarrow\infty}\frac{|p_{R}^{(t)}-\frac{1}{K}|}{% |p_{R}^{(t-1)}-\frac{1}{K}|^{2}}=K.$

Proof..

When the sizes of communities are equal, each community has $\frac{1}{K}$ nodes. Let $p_{R}^{(t-1)}=\frac{1}{K}+\theta$ where $\theta$ is a bias item. According to the analysis of the self-assembling process, by Eq. (4), we have the following expression:

$\displaystyle\lim\limits_{t\rightarrow\infty}\frac{|p_{R}^{(t)}-\frac{1}{K}|}{% |p_{R}^{(t-1)}-\frac{1}{K}|^{2}}=\frac{|2(\frac{1}{K}+\theta)-(\frac{1}{K}+% \theta)^{2}\times K-\frac{1}{K}|}{|\theta|^{2}}=K,$

which means that the convergence is quadratic. ∎

5.2 Parallelism

The parallelism is designed as follows. Let $S_{I}$ be the set of integers from 1 to $N$ , partitioned into non-overlapping and approximately equal-size subsets, each of which is represented by $S_{Ib}$ for the $b$ th processor. For each processor, we also use local variables $s_{b}$ , $\hat{s}_{b}$ and $\beta_{b}$ , which have the similar definitions as $s$ , $\hat{s}$ and $\beta$ , but taking only the $i$ th row $(i\in S_{Ib})$ of $W$ into account.

[b] Parallel Algorithm (Par-SBM)

set community number to $\alpha K$ initialize $B$ , $Z^{(t)}$ into many tiny communities

$Z^{(t-1)}=Z^{(t)}$

do in parallel $i\in S_{Ib}$ compute $N^{(d)}$ for node $i$ update $Z_{i:}^{(t)}$ according to Eq. (19) $Z_{i:}^{(t-1)}=Z_{i:}^{(t)}$

synchronize changes of $Z_{i:}$ at frequency $f$

update $s_{b}$ and $\hat{s}_{b}$ for each $\beta_{b}$ entries synchronize $s=\sum_{b}s_{b}$ and $\hat{s}=\sum_{b}\hat{s}_{b}$ update $B$ according to Eq. (16) update $N^{(c)}$ $Z^{(t)}$ does not change anymore merge communities according to Eqs (21) and (22)

no more merging is possible

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 $Z_{i:}$ are sent to processor $b$ only if $i$ has neighboring nodes in $S_{Ib}$ . 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 communication strategies when computing $B$ . Algorithm 5.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 Eqs (16) and (19) is ${\cal O}\left({\frac{M}{g}}\right)$ . When the computations for $B$ are partitioned on each processor, reducing the number of communities also requires ${\cal O}\left({\frac{M}{g}}\right)$ for each processor. The total message length for $Z$ and $B$ during the iterations is ${\cal O}\left({g(N+K)}\right)$ . If the bandwidth is constant, the speedup of the parallel algorithm is

$\displaystyle\frac{T(M)}{T(\frac{M}{g})+{\cal O}(N+K)},$ (23)

where $T$ is a linear function, such that the numerator and denominator indicate the run 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$ .

6. Experimental results

6.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) [6], a Bayesian inference algorithm (BI) [22], spectral clustering (SC) [45], label propagation (LP) [40], and modularity optimization (MM) [37]. We use the LFR benchmark generator [28] 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. $\mu$ is the proportion of inter-community connections to the intra-community ones.

The accuracy is evaluated by normalized mutual information (NMI) [13], 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 Fig. 3a and b we can see that our algorithm is relatively fast and has the best accuracy. When $\mu=0$ , most of the algorithms can find the correct result, and our algorithm is the second fastest one. As $\mu$ increases to 0.7, our algorithm needs 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. Figure 3c and d shows the results on graphs similar as the previous experiment but the average degree is 10 and $\mu=0.1$ . On large graphs, compared with other algorithms, ours consistently produces more accurate results with respect to the ground-truth setting, and the run time growth rate is smaller than others as the graph size grows. The slow growth of the run time of our algorithm is in accordance with the theoretical time complexity which is linear in the number of edges.

Figure 3.

Comparison to popular algorithms in MATLAB. (a) and (b) represent the variation of average accuracy and run time over different choices of $\mu$ . (c) and (d) represent the variation over different choices of graph size $N$ .

6.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) [6] and a fast graph partitioning algorithm Metis (MT) [25]. 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 [49], but the results are typically invariant. We use several metrics to compare the results, such as NMI, modularity [11], and conductance (the smaller, the better) [30]. The results on those quantities are averaged on five runs.

First, we run our algorithm on benchmark graphs (LFR-1e6) generated by LFR benchmark [28]. 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 $\mu$ and maximum community size $m c$ ) to have more variations. Figure 4 shows the results, in which GT is the ground truth given by the generator. Generally, the increase of $\mu$ and $m c$ 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. 4b, 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 as before while sacrificing the quality. In addition, the performance of our algorithm is stable as the run results using 8 processors (P8) and 64 processors (P64) are fairly close.

Table 1
Properties of the test graphs

Graph name	Node number	Edge number
LFR-1e6	1,000,000	$\sim$ 5,000,000
cit-Patents	3,774,768	16,518,948
dblp-Coauthor	1,314,050	10,724,828

Figure 4.

Comparison on benchmark data.

Table 2

Comparison on real-world data

Graph	Algo	Cond	Mod	Time
cit-Patents	PS1	0.418	0.631	1514
	PS8	0.416	0.638	284
	PS64	0.415	0.639	51
	LVF	0.594	0.560	103
	MT	0.488	0.487	243
dblp-Coauthor	PS1	0.114	0.631	186
	PS8	0.113	0.632	33
	PS64	0.112	0.632	6
	LVF	0.397	0.561	22
	MT	0.621	0.253	73

Figure 5.

Number of communities versus community size.

We also compare algorithms on real-world graphs: cit-Patents [29] and dblp-Coauthor [31]. 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 [16] that prefers unrealistically large communities. A similar problem happens on LV for real-world data as it generates many communities containing more than $10^{4}$ nodes (Fig. 5). Especially, for example, for dblp-Coauthor network, the biggest community by LV contains about $10^{5}$ 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 [30], so LV is omitted from the remaining comparisons. Alternatively, choosing low level structures such as LVF can be a remedy [18], 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, where run time is in seconds.

6.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 [6] 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. For example, we consider community containing author “Michael_Wooldridge,” and use $S_{PS}$ , $S_{\textit{LVF}}$ , and $S_{LV2}$ 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 $|S_{PS}|=255$ and $|S_{\textit{LVF}}|=166$ , while $|S_{LV2}|=26294$ .

All the algorithms are able to identify blue nodes in Fig. 6 as a separate community, while our Par-SBM provides the clearest detail. $S_{LV2}$ contains one hundred times more nodes than $S_{PS}$ . Among all the communities found by Par-SBM, the red, green, and purple communities in Fig. 6a are the top-three contributors to $S_{LV2}$ , owning 2.96% of total nodes in $S_{LV2}$ . 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 form 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. 6b.

Figure 6.

Subgraph of the dblp-Coauthor network, where colors represent the community labels determined by different algorithms.

Figure 7.

Comparison to the state-of-the-art algorithm by nuclear norm minimization. The curves represent the average smallest $p$ required by the algorithms to successfully recover the communities.

In addition, Par-SBM can find relevant nodes more than LVF does. For example, Par-SBM uniquely includes “H._H._Pham” and “Ali_M._Aseere” into the purple community in Fig. 6a. The corresponding authors have “agent” in the title of all their publications.

6.4 Comparing to nuclear norm minimization

In 2012, a nuclear norm minimization algorithm (NN) [8] is proposed to solve such SBM problem. It naturally comes with time complexity ${\cal O}(N^{3})$ due to the SVD decomposition, but it can find the global optimum theoretically and has the advantage of clustering sparse graphs. As their source code is not available, we run our serial algorithm with exactly the same setting as it is in [8] and compare the accuracy.

Table 3
Properties of the test graphs

Graph name	Node number	Edge number
Bitcoin	13,086,528	44,032,115
soc-LiveJournal1	4,847,571	68,993,773
Wikipedia	5,716,808	130,160,392

Figure 8.

Performance evaluation for graphs.

The test graph contains 1000 nodes with 5 communities; over different $q$ , we take the smallest $p$ that succeeds, and average over 20 trials; success means the misclassification for pairs of nodes is less than 0.1%.

Figure 7 compares the accuracy of our algorithm (ML) with NN, which is the best reported in [8]. This figure shows that our algorithm is very competitive to the nuclear norm minimization (NN), and the minimum possible gap between $p$ and $q$ for the two algorithms are similar. This fact indicates that our algorithm can find a local optimum that is very close to a global one with much less time complexity.

6.5 Strong scaling on large graphs

We use four different graphs. The first graph (LFR-1e6) is generated by LFR benchmark [28] using the same setting as in Fig. 3c and d except that the graph size is much larger. The second graph (Bitcoin) is a Bitcoin transaction network [27], where each node corresponds to one Bitcoin address, and each edge indicates that there is at least one transaction between two addresses. The third graph (soc-LiveJournal1 [30]) is obtained from Stanford Large Network Dataset Collection. The fourth graph (Wikipedia) contains most of the pages and the interlinks on Wikipedia in 2009.

The graph properties are listed in Column “Node Number” and “Edge Number” in Table 3. These graph data are preprocessed and partitioned on each processor so that different processors will have a similar number of nodes and edges.

We run our algorithm on these graphs using different number of processors, but the accuracy is quite stable as expected. For example, when running on LFR-1e6 with ground-truth cluster information, the resulted NMI by our algorithm is in the range of 0.997 $\pm$ 0.002, indicating that the detected clusters are very similar to the ground truth, regardless the number of processors.

Figure 8 demonstrates the run time and the speedup. Our algorithm on all the graphs can achieve a significant speedup as the number of processors grows. The Wikipedia graph has the highest average degree in all the test graphs, and Fig. 8c shows that our algorithm can achieve the biggest speedup of 48x on 128 cores. This phenomenon is in accordance with Eq. (23).

We also compare the speedup performance with a recently proposed parallel community detection algorithm [41, 42] running on shared-memory systems. Their community detection algorithm achieves 5.1x speedup on 16 cores using Intel Xeon X5570 CPU and 9.4x speedup using 32 cores on XMT2 [42] for soc-LiveJournal1 data.

Using the same data set and the same Intel Xeon CPU, our algorithm achieves an 8.2x and 12.8x speedup on 16 and 32 cores respectively, as shown in Fig. 8d. Furthermore, 25.7x speedup can be achieved on 128 cores.

7. 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 run time shortages of traditional algorithms. The code is available at https://github.com/pdacorn/par-sbm.

Footnotes

Acknowledgments

Research reported in this publication was supported by King Abdullah University of Science and Technology (KAUST).

References

Abbe

Bandeira

A.S.

and Hall

, Exact recovery in the stochastic block model, IEEE Transactions on Information Theory 62(1) (2016), 471–487.

Carlos

Jesús

G.-C.

and Jordi

, The community structure of SAT formulas, In Theory and Applications of Satisfiability Testing–SAT 2012, Springer, 2012, 410–423.

Asanovic

Bodik

Catanzaro

B.C.

Gebis

J.J.

Husbands

Keutzer

Patterson

D.A.

Plishker

W.L.

Shalf

Williams

S.W.

et al., The landscape of parallel computing research: A view from Berkeley. Technical report, Technical Report UCB/EECS-2006-183, EECS Department, University of California, Berkeley, 2006.

Bhowmick

and Srinivasan

, A template for parallelizing the louvain method for modularity maximization, In Dynamics On and Of Complex Networks, Springer, 2 (2013), 111–124.

Biesan

Anthony

and desJardins

, Block modeling in large social networks with many clusters, In AAAI Fall Symposium Series, 2012.

Blondel

V.D.

Guillaume

J.-L.

Lambiotte

and Lefebvre

, Fast unfolding of communities in large networks, Journal of Statistical Mechanics: Theory and Experiment 2008(10) (2008), P10008.

Celisse

Daudin

J.-J.

and Pierre

, Consistency of maximum-likelihood and variational estimators in the stochastic block model, Electronic Journal of Statistics 6 (2012), 1847–1899.

Chen

Sanghavi

and Xu

, Clustering sparse graphs, In Advances in Neural Information Processing Systems 25, 2012, 2213–2221.

Chin

Rao

and Vu

, Stochastic block model and community detection in the sparse graphs: A spectral algorithm with optimal rate of recovery, arXiv preprint arXiv:1501.05021 2(4) (2015).

10.

Chung

and Lu

, Complex Graphs and Networks (Cbms Regional Conference Series in Mathematics), American Mathematical Society, Boston, MA, USA, 2006.

11.

Clauset

Newman

M.E.J.

and Moore

, Finding community structure in very large networks, Physical Review E 70(6) (2004), 066111.

12.

Dagum

and Menon

, Openmp: an industry standard api for shared-memory programming, IEEE Computational Science & Engineering 5(1) (1998), 46–55.

13.

Danon

Diaz-Guilera

Duch

and Arenas

, Comparing community structure identification, Journal of Statistical Mechanics: Theory and Experiment 2005(09) (2005), P09008.

14.

Daudin

J.J.

Picard

and Robin

, A mixture model for random graphs, Statistics and Computing 18(2) (2008), 173–183.

15.

Diaz

Munoz-Caro

and Nino

, A survey of parallel programming models and tools in the multi and many-core era, IEEE Transactions on Parallel and Distributed Systems 23(8) (2012), 1369–1386.

16.

Fortunato

and Barthelemy

, Resolution limit in community detection, Proceedings of the National Academy of Sciences 104(1) (2007), 36–41.

17.

Fortunato

and Castellano

, Community structure in graphs, In Computational Complexity, Springer, 2012, 490–512.

18.

Good

B.H.

de Montjoye

Y.-A.

and Clauset

, Performance of modularity maximization in practical contexts, Physical Review E 81(4) (2010), 046106.

19.

Grama

Gupta

Karypis

and Kumar

, Introduction to Parallel Computing, Second Edition, Pearson Education, 2003.

20.

Gregori

Lenzini

and Mainardi

, Parallel k-clique community detection on large-scale networks, IEEE Transactions on Parallel and Distributed Systems 24(8) (2013), 1651–1660.

21.

Gropp

Lusk

and Skjellum

, Using MPI: portable parallel programming with the message-passing interface, volume 1. MIT Press, 1999.

22.

Hofman

J.M.

and Wiggins

C.H.

, Bayesian approach to network modularity, Physical Review Letters 100(25) (2008), 258701.

23.

Holland

P.W.

Laskey

K.B.

and Leinhardt

, Stochastic blockmodels: First steps, Social Networks 5(2) (1983), 109–137.

24.

Kang

Meeder

Papalexakis

E.E.

and Faloutsos

, Heigen: Spectral analysis for billion-scale graphs, IEEE Transactions on Knowledge and Data Engineering 26(2) (2014), 350–362.

25.

Karypis

and Kumar

, A fast and high quality multilevel scheme for partitioning irregular graphs, SIAM Journal on Scientific Computing 20(1) (1998), 359–392.

26.

Kemp

Tenenbaum

J.B.

Griffiths

T.L.

Yamada

and Ueda

, Learning systems of concepts with an infinite relational model, In AAAI 3 (2006), 5.

27.

Kondor

Psfai

Csabai

and Vattay

, Do the rich get richer? an empirical analysis of the bitcoin transaction network, PLoS ONE 9(2) (2014), e86197.

28.

Lancichinetti

and Fortunato

, Benchmarks for testing community detection algorithms on directed and weighted graphs with overlapping communities, Physical Review E 80(1) (2009), 016118.

29.

Leskovec

Kleinberg

and Faloutsos

, Graphs over time: densification laws, shrinking diameters and possible explanations, in: Proceedings of the 11th ACM International Conference on Knowledge Discovery in Data Mining, ACM, 2005, pp. 177–187.

30.

Leskovec

Lang

K.J.

Dasgupta

and Mahoney

M.W.

, Community structure in large networks: Natural cluster sizes and the absence of large well-defined clusters, Internet Mathematics 6(1) (2009), 29–123.

31.

Ley

, The DBLP computer science bibliography: Evolution, research issues, perspectives, in: Proceedings of the International Symposium on String Processing and Information Retrieval, 2002, pp. 1–10.

32.

Pan

Xiao

, and Huang

, Community discovery and importance analysis in social network, Intelligent Data Analysis 18(3) (2014), 495–510.

33.

and Schuurmans

, Modular community detection in networks, In Twenty-Second International Joint Conference on Artificial Intelligence, 2011.

34.

Liu

Wang

Danilevsky

and Han

, Large-scale spectral clustering on graphs, in: Proceedings of the 23rd International Joint Conference on Artificial Intelligence, AAAI Press, 2013, pp. 1486–1492.

35.

Mitchell

, Complex systems: Network thinking, Artificial Intelligence 170(18) (2006), 1194–1212.

36.

Nadler

and Galun

, Fundamental limitations of spectral clustering, In: Advances in Neural Information Processing Systems, 2006, pp. 1017–1024.

37.

Newman

M.E.J.

, Fast algorithm for detecting community structure in networks, Physical Review E 69(6) (2004), 066133.

38.

Newman

M.E.J.

, Modularity and community structure in networks, Proceedings of the National Academy of Sciences 103(23) (2006), 8577–8582.

39.

Rabenseifner

Hager

and Jost

, Hybrid mpi/openmp parallel programming on clusters of multi-core smp nodes, in: 17th Euromicro International Conference on Parallel, Distributed and Network-based Processing, IEEE, 2009, pp. 427–436.

40.

Raghavan

U.N.

Albert

and Kumara

, Near linear time algorithm to detect community structures in large-scale networks, Physical Review E 76(3) (2007), 036106.

41.

Riedy

E.J.

Meyerhenke

Ediger

and Bader

D.A.

, Parallel community detection for massive graphs, in: Parallel Processing and Applied Mathematics, Springer, 2012, pp. 286–296.

42.

Riedy

Bader

D.A.

and Meyerhenke

, Scalable multi-threaded community detection in social networks, in: the 26th International Parallel and Distributed Processing Symposium Workshops & PhD Forum (IPDPSW), IEEE, 2012, pp. 1619–1628.

43.

Staudt

and Meyerhenke

, Engineering parallel algorithms for community detection in massive networks, IEEE Transactions on Parallel and Distributed Systems PP(99) (2015), 1–1.

44.

Traag

V.A.

Van Dooren

and Nesterov

, Narrow scope for resolution-limit-free community detection, Physical Review E 84(1) (2011), 016114.

45.

Luxburg

U.V.

, A tutorial on spectral clustering, Statistics and Computing 17(4) (2007), 395–416.

46.

Wang

Zhang

Tsang

Chen

and Zhang

, Defragging subgraph features for graph classification, in: Proceedings of the 24th ACM International on Conference on Information and Knowledge Management, CIKM ’15, New York, NY, USA, 2015, pp. 1687–1690. ACM.

47.

Wang

and Zhao

, A tale of three social networks: User activity comparisons across facebook, twitter, and foursquare, IEEE Internet Computing 18(2) (Mar 2014), 10–15.

48.

White

and Smyth

, A spectral clustering approach to finding communities in graphs, in: Proceedings of the Fifth SIAM International Conference on Data Mining, vol. 119, 2005, pp. 274. SIAM.

49.

Wickramaarachchi

Frincu

Small

and Prasanna

, Fast parallel algorithm for unfolding of communities in large graphs, in: 18th IEEE High Performance Extreme Computing Conference (HPEC 14), 2014, pp. 1–6.

50.

Xiang

and Hu

, Limitation of multi-resolution methods in community detection, Physica A: Statistical Mechanics and its Applications 391(20) (2012), 4995–5003.

51.

Yang

Zhao

Zheng

V.W.

and Li

X.L.

, An aggressive graph-based selective sampling algorithm for classification, in: Data Mining (ICDM), 2015 IEEE International Conference on, Nov 2015, pp. 509–518.

52.

Zhang

and Newman

M.E.J.

, Multiway spectral community detection in networks, Physical Review E 92(5) (2015), 052808.

A scalable community detection algorithm for large graphs using stochastic block models

Abstract

Keywords

1. Introduction

2. Related work

3. Stochastic block model

4.1 Self-assembling process

5.1 Estimation algorithm

Proof..

Proof..

Proof..

5.2 Parallelism

6.1 The serial algorithm

Table 1 Properties of the test graphs

Table 3 Properties of the test graphs

7. Concluding remarks

Footnotes

Acknowledgments

References

Table 1
Properties of the test graphs

Table 3
Properties of the test graphs