Hierarchical community detection via rank-2 symmetric nonnegative matrix factorization

Background Community discovery is an important task for revealing structures in large networks. The massive size of contemporary social networks poses a tremendous challenge to the scalability of traditional graph clustering algorithms and the evaluation of discovered communities. Methods We propose a divide-and-conquer strategy to discover hierarchical community structure, nonoverlapping within each level. Our algorithm is based on the highly efficient rank-2 symmetric nonnegative matrix factorization. We solve several implementation challenges to boost its efficiency on modern computer architectures, specifically for very sparse adjacency matrices that represent a wide range of social networks. Conclusions Empirical results have shown that our algorithm has competitive overall efficiency and leading performance in minimizing the average normalized cut, and that the nonoverlapping communities found by our algorithm recover the ground-truth communities better than state-of-the-art algorithms for overlapping community detection. In addition, we present a new dataset of the DBLP computer science bibliography network with richer meta-data and verifiable ground-truth knowledge, which can foster future research in community finding and interpretation of communities in large networks.

Earlier work in this area mainly focused on finding linkage-based communities on small networks. Various measures have been defined and optimized for the difference between intraconnections and interconnections among network nodes, such as normalized cut, conductance, and others, using a variety of methods [4]. Manual examinations (such as visualization) of the results have typically been used in those studies, which may provide some valuable information and insights. Community detection on large networks is more challenging due to the following reasons: (1) many algorithms suitable for small networks are often not scalable; (2) there is a dearth of ground-truth communities defined for large networks and even for the datasets with ground truth, where the quality of the ground truth is often questionable; and (3) examining the results is nearly impossible. It has not been until recently [4] that ground-truth communities have been defined and studied in several real-world large-scale networks.
This paper introduces a scalable algorithm based on rank-2 symmetric nonnegative matrix factorization (rank-2 SymNMF) for large-scale hierarchical community detection. After summarizing some related works, we discuss the problem domain for our new results and our solutions for particular problems within that domain. Then, some highlights and speedup in our implementations are discussed, after which we present comprehensive experimental results to evaluate our new algorithm.

Related work
The study of network community detection dates back to the well-known Kernighan-Lin algorithm from the early 1970s [5]. At that time, the network community detection problem was often formulated as a graph-partitioning problem, which aims at "dividing the vertices" into a predefined number of nonoverlapping "groups of predefined size, such that the number of edges lying between the groups is minimal" [6]. Many methods that produce good-quality solutions were proposed, but they were based on combinatorial optimization algorithms and were not scalable. Later, when it was discovered that graph partitioning is an important problem for balanced distribution of work loads in parallel computing, computer scientists developed many algorithms, such as METIS [7], SCOTCH [8], and Chaco [9], for graph partitioning of parallel communication problems. These algorithms usually follow a multilevel strategy, where a large graph is first coarsened to a smaller graph by recursively contracting multiple vertices into one vertex, and then the small graph is partitioned applying a method like the Kernighan-Lin algorithm, and finally the partition is mapped back to the original graph with some refinement. Most of these algorithms (e.g., all three we mentioned above) scale well to very large networks containing millions of nodes.
A spectral clustering method that minimizes normalized cut was proposed as an image-segmentation algorithm [10], and it soon became popular in the area of graph clustering. However, due to the time-consuming eigenvector/singular vector computation in this algorithm, it is not scalable to the case when the number of communities is large. The Graclus algorithm [11] by Dhillon, Guan, and Kulis solved this issue by utilizing the mathematical equivalence between general cut or association objectives (including normalized cut and ratio association) and the weighted kernel k-means objective [12] and applying a multilevel framework. Kuang, Ding, and Park discovered that the SymNMF (symmetric nonnegative matrix factorization) objective function is also equivalent to normalized cut and ratio association objective functions with a relaxation different from that in spectral clustering [13,14]. This algorithm has better interpretability like many other NMF-based methods.
Girvan and Newman [15] produced pioneering work developing graph-partitioning/ clustering methods from a community detection viewpoint, which finds "groups of vertices which probably share common properties and/or play similar roles within the graph" [6]. Around that time period, many new algorithms were invented. Later, Newman and Girvan [16] proposed the modularity measurement for community detection, on which the biggest family of community detection algorithms is based [17]. A scalable example in this family of algorithms is the Louvain algorithm [18]. Several algorithms such as Walktrap [19] and Infomap [20] are based on random walk on graphs, with the idea that in a random walk, the probability of staying inside a community is higher than going to another community. The paper [6] provides a comprehensive review of the algorithms that appeared up to 2010.
The early overlapping community detection algorithms [21,22] were not effective on large graphs. Lancichinetti et al. [23] proposed a scalable overlapping community detection algorithm-order statistics local optimization method (OSLOM), which was based on a measurement similar to modularity but was able to handle overlapping communities. Yang and Leskovec studied the properties of large-scale overlapping communities [4] and proposed the BigClam algorithm [1]. They provided some large-scale datasets with ground-truth communities available to researchers, which have become standard test datasets. The BigClam algorithm seeks to fit a probabilistic generative model that satisfies certain community properties discovered in their studies [1,24]. Whang et al. [2] proposed another overlapping community detection algorithm called NISE, based on seed set expansion, which starts with a seed set generated by Graclus or other methods and uses random walk to obtain overlapping communities. These algorithms that are dedicated to community detection have demonstrated better performance in terms of discovering ground-truth communities compared with the traditional graph-partitioning algorithms.
Recently, [17] proposed a new nonoverlapping community detection algorithm, scalable community detection (SCD). Network communities of good quality should have stronger intraconnections than interconnections. Previous algorithms measure such strength of connectivity only through the number of edges. The uniqueness of SCD is that it is based on a triangular structure of edges. The goal is to identify communities where each node forms more triangles with nodes inside the community than those formed with nodes outside the community.
On the other hand, nonnegative matrix factorization (NMF)-based methods exhibit superior interpretability in many application areas such as text mining [25,26] and image analysis [27,28]. In this paper, we will show that our NMF-based algorithm has competitive performance and scalability for community detection. Although our algorithm currently only handles nonoverlapping community detection, it has achieved comparable or even better-quality than the state-of-the-art overlapping community detection algorithms (such as BigClam) in our extensive tests. Our algorithm is inspired by SymNMF [13,14] for graph clustering and HierNMF2 [25] for fast document clustering.

Problem definition
As mentioned above, there is no universally accepted definition for network communities. Rather than defining community detection as optimizing some specific measurement criteria, we believe it would be more effective and flexible to understand the problem by clearly defining what makes a community detection result good.
In this paper, we focus on link-based community detection, and use functional communities (if known) as ground truth. In general, to evaluate link-based community detection results, we may ask several questions. The first is whether the result is coherent from the point of view of network links (Q1). There are many measurement scores defined based on the presumption that a community should have more intraconnections than interconnections. These measures include normalized cut, ratio cut, conductance, etc. The second question is whether the result agrees with prior knowledge, especially the ground truth if it is known (Q2). There have been largely two approaches. One is to manually analyze the results with some known meta-information (such as the entity each node represents), which is not scalable to large networks. Another approach is to compare the community detection result using some measures such as F1 score, which assumes the existence of ground truth. Finally, we would like to know whether the result reveals some new and useful information about the network (Q3). This is mostly relevant for the study of small networks [15]. For large networks, it is almost impossible to manually check all communities discovered. However, the answers to Q2 may guide our focus to more interesting parts of the network.

Ground-truth communities
Ground-truth communities of large networks were not available to researchers until Yang and Leskovec [1] defined ground-truth communities (as found in SNAP 1 ) for several real-world large networks, including several social networks, paper coauthorship networks, and product copurchase networks. In their work, the ground-truth communities are defined using functional communities already present in the data. For example, in social networks, user groups can be treated as communities; in paper coauthorship networks, two authors publishing in the same venue can be seen to be in the same community; in product copurchase networks, product category can be naturally used as communities. These functional communities are not necessarily directly related to network structures. For example, the product category is an inherent property of a product, which can never be affected by copurchasing activities. Therefore, it is not reasonable to expect that link-based community detection algorithms can fully recover functional communities. On the other hand, many studies show that there are close relations between link-based communities and functional communities. The paper [4] shows that many linkage-based measurements (such as normalized cut, conductance, etc.) also have good performance on functional communities. Also, the results of link-based community detection algorithms can sometimes recover functional communities.
Based on the above observations, we conclude that link-based community detection algorithms have the ability to partially recover functional communities, but such ability is inherently limited by the essential differences between functional communities and link-based communities. This should be kept in mind when comparing link-based community detection results against ground-truth communities.

Overlapping vs nonoverlapping communities
Real-world network communities are usually overlapping. For example, it is common that one user joins a variety of groups in a social network. However, our current focus is on nonoverlapping community detection, since nonoverlapping community detection is also very useful for revealing the network structure, and our algorithm is designed to detect nonoverlapping communities efficiently. The results of a good-quality nonoverlapping community detection algorithm can be used as an effective starting point for overlapping community detection [2,3].

Hierarchical rank-2 symmetric NMF
We present an algorithm called HierSymNMF2 for hierarchical community detection. HierSymNMF2 uses a fast SymNMF algorithm [14] with rank 2 (SymNMF2) for binary community detection and recursively apply SymNMF2 to further binary split one of the communities into two communities in each step. This process is repeated until a preset number of communities is discovered, or there are no more communities that are worthy of any further binary split. Our approach starts with a low rank approximation (LRA) of the data based on the nonnegative matrix factorization (NMF), which reduces the dimension of the data while keeping key information. In addition, the results of NMF-based methods directly provide information regarding the assignment of data to clusters/communities.
Given the vast amounts of nonnegative data available for extracting critical information, the NMF has found a wealth of applications in such domains as image, text, and chemical data processing. It can be shown that applying algorithms to such data without constraining the solution can result in uninterpretable results such as negative chemical concentrations and possibly false negative and/or false positive detections, which could lead to meaningless results [29]. For text analytics, a corpus of text documents can be represented by a nonnegative term-document matrix. Likewise, for graph analytics, the nonnegative adjacency matrix is used as an input to NMF algorithms. NMF seeks an approximation of such nonnegative matrices with a product of two nonnegative low rank matrices. With various constraints and regularization terms on the NMF objective function, there are many variants of NMF, which are appropriate for a large variety of problem domains. A common formulation of NMF is the following: where X ∈ R m×n + , W ∈ R m×k + , H ∈ R k×n + (R + is the set of all real nonnegative numbers), and k ≪ min(m, n). In this formulation, each data item is represented by a column of the matrix X, and each column in the matrix H can be seen as a low rank representation of the data item. Nonnegativity constraints allow such a low rank representation to be more interpretable than other low rank approximations such as SVD. This formulation can be applied to areas such as document clustering [26] and can be solved efficiently for very large m and n [30]. However, when k reaches a value on the order of thousands, NMF algorithms become slow. To solve this issue, [25] developed a divide-and-conquer method that relies on rank-2 NMF, where k = 2, which exhibits significant speedups. The framework of this divide-and-conquer method is shown in Algorithm 1. In this divide-and-conquer framework, the task of splitting one cluster into two clusters is performed by rank-2 NMF, which reduces the superlinear time complexity with respect to k to linear [25].
A variant of NMF, SymNMF [13,14], which is the symmetric version of NMF, can be used for graph clustering. The formulation of SymNMF is where S ∈ R n×n is a symmetric similarity matrix of graph nodes: H ∈ R n×k + and k ≪ n . Some choices of the input matrix S for SymNMF are adjacency matrix S G and normalized adjacency matrix D −1/2 S G D −1/2 , where D = diag(d 1 , . . . , d n ), and d i = n j=1 S G ij is the degree of node i. When S is the adjacency matrix, (2) is a relaxation of maximizing the ratio association; when S is the normalized adjacency matrix, (2) is a relaxation of minimizing the normalized cut [13] (see Appendices A, B for a complete proof ). Sym-NMF is an effective algorithm for graph clustering, but for large k, improvements in computational efficiency are necessary.
The algorithm we introduce in this paper uses the framework shown in Algorithm 1, where a cluster is a community, and the task of splitting a community is performed by our rank-2 version of SymNMF. The decision to choose the next node to split is based on a criterion discussed in the next section. In the following sections, we denote S as the similarity matrix representing a graph G, and S c as the matrix representation of a community, i.e., a subgraph of G (the corresponding submatrix of S).

Splitting a community using rank-2 SymNMF
Splitting a community is achieved by rank-2 SymNMF of S c ≈ HH T where H ∈ R n×2 + . The result H naturally induces a binary split of the community: suppose H = (h ij ), then where c i is the community assignment of the ith graph node.
A formal formulation of rank-2 SymNMF is the following optimization problem: where H ∈ R n×2 + . This is a special case of SymNMF when k = 2, which can be solved by a general SymNMF algorithm [13,14]. However, by combining the alternating nonnegative least squares (ANLS) algorithm for SymNMF from [14] and the fast algorithm for rank-2 NMF from [25], we can obtain a fast algorithm for rank-2 SymNMF.
First, we rewrite (3) into asymmetric form plus a penalty term [31] as follows: where W and H ∈ R n×2 + , and α > 0 is a scalar parameter for the tradeoff between the approximation error and the difference between W and H. Formulation (4) can be solved using a two-block coordinate descent framework, alternating between the optimizations for W and H. When we solve for W, (4) can be reformulated as where I 2 is the 2 × 2 identity matrix. Similarly, when we solve for H, (4) can be reformulated as We note that both (5) and (6) are in the following form: where F ∈ R m×2 + , G ∈ R m×n + . This formulation can be efficiently solved by an improved active-set-type algorithm described in [25], which we call rank2nnls-fast. The idea behind rank2nnls-fast can be summarized as follows: the optimization problem (7) can be decomposed into n independent subproblems in the following form: To solve (8) efficiently, we note that when g � = 0, there will be only three possible cases for y = [y 1 , y 2 ], where only one of y 1 and y 2 is 0 or both are positive. These three cases can easily be solved by the usual least-squares algorithms, e.g., normal equations. Details can be found in Algorithm 2 in [25].

Choosing a node to split based on normalized cut
The "best" community to split further is chosen by computing and comparing splitting scores for all current communities corresponding to the leaf nodes in the hierarchy. The proposed splitting scores are based on normalized cut. We make this choice because (1) normalized cut determines whether a split is structurally effective since it measures the difference between intraconnections and interconnections among network nodes; and (2) for SymNMF, when S is the normalized adjacency matrix, the SymNMF objective function is equivalent to (a relaxation of ) minimizing the normalized cut, which is the preferred choice in graph clustering [14].
Suppose we have a graph G = (V , E), where the weight of an edge (u, v) is w(u, v). Note that for an unweighted graph, where which measures the number of edges inside the subgraph induced by A i (intraconnection); and measures the number of edges between A i and the remaining nodes in the graph (interconnection). Note that in the definition of within(A i ) (10), each edge within A i is counted twice. In the special case, k = 2, we have From Eq. (9), it is evident that when each community has many more intraconnections than interconnections, there is a small normalized cut.
For example, the graph shown in Fig. 1 originally has three communities A 1 , A 2 and A 3 , and the corresponding normalized cut is The community A 3 is now split into two smaller communities B 1 and B 2 and normalized cut can be used to measure the goodness of this split. We consider three possibilities: (1) isolate A 3 and compute normalized cut of the split as where the subscript A 3 means only consider the edges inside A 3 . We denote the above criterion by ncut_local. (2) A more global criterion is to also consider the edges that go across A 3 : This criterion is denoted by ncut_global. (3) Minimize the global normalized cut using a greedy strategy. Specifically, choose the split that results in the minimal increase in the global normalized cut: We denote this criterion by ncut_global_diff and will compare the performance of these three criteria in later sections.

Implementation
In the previous work on rank-2 NMF [32] that takes a term-document matrix as input in the context of text clustering, sparse-dense matrix multiplication (SpMM) was the main computational bottleneck for computing the solution. However, this is not the case with rank-2 SymNMF or HierSymNMF2 for community detection problems on typical large-scale networks. Suppose we have an n × n adjacency matrix with z nonzeros as an input to rank-2 SymNMF. In Algorithm 2, i.e., Nonnegative Least Squares (NLS) with two unknowns, SpMM costs 2z floating-point operations (flops), while searching for the optimal active set (abbreviated as opt-act) costs 12n flops. Of the 12n flops for optact, 8n flops are required for solving n linear systems each of size 2 × 2 corresponding to line 1 in Algorithm 2, and the remaining 4n flops are incurred by lines 4-5. Note that comparison operations and the memory I/O required by opt-act are ignored. Fig. 1 A graph for illustrating normalized cut and our splitting criteria. The structure of the graph is inspired by Figure 1 from [15] The above rough estimation of computational complexity reveals that if z ≤ 6n, or equivalently, if each row of the input adjacency matrix contains no more than 6 nonzeros on average, then SpMM will not be the major bottleneck of the rank-2 SymNMF algorithm. In other words, when the input adjacency matrix is extremely sparse, which is the typical case we have seen on various datasets (Table 1), then further acceleration of the algorithmic steps in opt-act will achieve higher efficiency. Figure 2 (upper) shows the proportions of runtime corresponding to SpMM, optact, and other algorithmic steps implemented in Matlab, which demonstrate that both SpMM and opt-act are the targets for performance optimization.

Multithreaded SpMM
SpMM is a required routine in lines 1, 4, and 5 of Algorithm 2. The problem can be written as where A ∈ R n×n is a sparse matrix and X, Y ∈ R n×k are dense matrices. 2 Most open-source and commercial software packages for sparse matrix manipulation have a single-threaded implementation for SpMM, for example, Matlab 3 , Eigen 4 , and Armadillo 5 (the same is also true for SpMV, sparse matrix-vector multiplication). For the Intel Math Kernel Library 6 , while we are not able to view the source, our simple tests have shown that it can exploit only one CPU core for computing SpMM. Part of the reason for the lack of parallel implementation of SpMM in generic software packages is that the best implementation for computing SpMM for a particular matrix A depends on the sparsity pattern of A.
In this paper, we present a simple yet effective implementation to compute SpMM for a matrix A that represents an undirected network. We exploit two important facts in order to reach high performance: • Since the nodes of the network are arranged in an arbitrary order, the matrix A is not assumed to have any special sparsity pattern. Thus, we can store the matrix A in the commonly used generic storage, the compressed sparse column (CSC) format, as is 2 The more general form of SpMM is Y ← Y + A · X. Our algorithm only requires the more simplistic form Y ← A · X, and thus, for this case, we wrote a specific routine saving n operations for addition. 3  practiced in the built-in sparse matrix type in Matlab. As a result, nonzeros of A are stored column-by-column.
• The matrix A is symmetric. This property enables us to build an SpMM routine for A T X to compute AX.
The second fact above is particularly important: When A is stored in the CSC format, computing AX with multiple threads would incur atomic operations or mutex locks to avoid race conditions between different threads. Implementing multithreaded A T X is much easier, since A T can be viewed as a matrix with nonzeros stored row-by-row, and we can divide the rows of A T into several chunks and compute the product of each row  chunk with X on one thread. Our customized SpMM implementation is described in Algorithm 3.
In addition, the original adjacency matrix often has the value "1" as every nonzero entry, that is, all the edges in the network carry the same weight. Thus, multiplication operations are no longer needed in SpMM with such a sparse matrix. Therefore, we have developed a specialized routine for the case where the original adjacency matrix is provided as input to HierSymNMF2.

C/Matlab hybrid implementation of opt-act
The search for the optimal active set, opt-act, is the most time-consuming step in the algorithm for NLS with two columns (Fig. 2 (lower)) when the input matrix is extremely sparse. Our overall program was written in Matlab, and the performance of the optact portion was optimized with native C code. The optimization exploits multiple CPU cores using OpenMP, and the software vectorization is enabled by calling AVX (advanced vector extensions) intrinsics.
It turns out that a C/Matlab hybrid implementation is the preferred choice for achieving high performance with native C code. Intel CPUs are equipped with AVX vector registers, since the Sandy Bridge architecture and these vector registers are essential for applying the same instructions to multiple data entries (known as instruction-level parallelism or SIMD). For example, a 256-bit AVX register can process four double-precision floating point numbers (64-bit each) in one CPU cycle, which amounts to four times speed-up over a sequential program. AVX intrinsics are external libraries for exploiting AVX vector registers in native C code. These libraries are not part of the ANSI C standard but retain mostly the same interface on various operating systems (Windows, Linux, etc). However, to obtain the best performance from vector registers, functions in the AVX libraries often require the operands having aligned memory addresses (32-byte aligned for double precision numbers). The function calls for aligned memory allocation, which is completely platform dependent for native C code, means that our software would not be easily portable across various platforms if aligned memory allocation were managed in the C code. Therefore, in order to strike the right balance between computational efficiency and software portability, our strategy is to allocate memory within Matlab for the vectors involved in opt-act, since Matlab arrays are memory aligned in a cross-platform fashion.
Finally, note that the opt-act step in lines 1, 4, and 5 of Algorithm 2 contains several division operations, which cost more than 20 CPU cycles each and are much more expensive than multiplication operations (1 CPU cycle). This large discrepancy in time cost would be substantial for vector-scalar operations. Therefore, we replace vectorscalar division, in the form of x/α where x is a vector and α is a scalar, by vector-scalar multiplication, in the form of x · (1/α).

Methods for comparison
We compare our algorithm with some recent algorithms mentioned in the "Related work" section. We use eight threads for all methods that support multithreading. For NISE, we are only able to use one thread because its parallel version exits with errors in our experiments. For all the algorithms, default parameters are used if not specified. To better communicate the results, below are the labels that denote each algorithm, which will be used in the following tables: • h2-n(g)(d)-a(x): These labels represent several versions of our algorithm. Here h2 stands for HierSymNMF2, n for the ncut_local criterion, ng for the ncut_ global criterion, and ngd for the ncut_global_diff criterion (see previous sections for the definitions of these criteria); 'a' means that we compute the real normalized cut using the original adjacency matrix; and 'x' indicates that an approximated normalized cut is computed using the normalized adjacency matrix, which usually results in faster computations. We stop our algorithm after k − 1 binary splits where k is the number of communities to find. Theoretically, this will generate k communities. However, we remove fully disconnected communities, as outliers since they are often far from being significant because of their unusually small size and they correspond to all-zero submatrices in the graph adjacency matrix, which does not have a meaningful rank-2 representation. Therefore, the final number of communities are usually slightly smaller than k, as will be shown in "Experiment results" section.
• NISE: An improved version of NISE that is published in 2016 [3].

Internal measures: average normalized cut/conductance
Normalized cut (9) is a measurement of the extent that communities have more intraconnections than interconnections and is shown to be an effective score [4]. Since our algorithm implicitly minimizes the normalized cut, it is natural to use normalized cut as an internal measure of community/clustering quality. One drawback of normalized cut is that it tends to increase when the number of communities increases. In Appendix B, we prove that the normalized cut strictly increases when one community is split into two. In practice, we observed that the normalized cut increases almost linearly with respect to the number of communities. Some community detection algorithms automatically determine the number of communities; hence, it is not fair to compare normalized cut for such algorithms against others that detect a preassigned number of communities. Therefore, it makes more sense to use the average normalized cut, i.e., the normalized cut divided by the number of communities. In addition, since the average normalized cut can be treated as a per-community property, it also applies to overlapping communities. Given k communities A 1 , . . . , A k (which may be overlapping), we define the average normalized cut as Conductance [33], which is shown to be an effective measure [4], is defined for a community as Conductance Hence the average normalized cut is actually equal to the average conductance (per community).

External measures: precision, recall, and F-score
Alternatively, we can measure the qualities of detected communities by comparing them with ground truth. Suppose k communities A 1 , . . . , A k were detected, and the ground truth has k ′ communities B 1 , . . . , B k ′. We compute the confusion matrix C = (c ij ) k×k ′, where c ij = |A i ∩ B j |. Then, pairwise scores can be defined as Although a global best match (i.e., finding a one-to-one mapping) between detected communities and ground-truth communities would be ideal, finding such a match is time consuming. We used per-community best match as a heuristic alternative. Specifically, we define the average F 1 score [1] as The average precision and average recall can be defined in a similar way. When comparing detected communities against ground truth, we remove nodes without ground-truth labels from the detected communities to achieve meaningful comparisons.

Datasets
The data used for the experimental results of this paper are mostly from SNAP datasets [4,34]. In our study, we found that the ground-truth information in SNAP is incomplete; for example, a large percentage of nodes do not belong to any ground-truth community. Table 1 shows some statistics regarding the number of communities to which each node belongs. Although all of these datasets can be conveniently accessed on the SNAP website as a graph with ground-truth communities, DBLP06 is the only dataset with a complete raw dataset openly available to the public. The other five datasets (Youtube, Amazon, Live-Journal, Friendster, and Orkut) were obtained by crawling the web, and they are far from being complete. Crawling large complex graphs is challenging by itself that may need extensive and specialized research efforts. We do not aim to solve this issue in this paper. The Orkut and Youtube datasets can be acquired from [35]. Detailed descriptions are available explaining the crawling procedure and analysis of the completeness. It has been concluded that the Orkut and Youtube datasets are not complete. Such incompleteness in crawled datasets is expected due to intrinsic restrictions of web crawling such as rate limit and privacy protection. The Friendster data were crawled by the ArchiveTeam, and the LiveJournal data come from [36]. The Amazon data were crawled by the SNAP group [37]. However, information on how the data were collected and processed and the information on analysis of data completeness are not available.
Possible reasons that many nodes in these datasets do not belong to any communities are (1) SNAP removed communities with less than three nodes, which caused some nodes to "lose" their memberships; (2) the well-known incompleteness of crawled datasets; (3) for social networks (Youtube, LiveJournal, Friendster, and Orkut), it is common that a user does not join any user groups; (4) SNAP used the dataset from [36] to generate the DBLP06 dataset, which was published in 2006. At that time, the DBLP database was not as mature and complete as it is today. Another issue of the above datasets is that all nodes are anonymized, which ensures protection of user privacy, but limits our ability to interpret community detection results.
The DBLP data are openly accessible, and are provided using a highly structured format-XML. We reconstructed the coauthorship network and ground-truth communities from a recent DBLP snapshot to obtain a more recent and complete DBLP dataset with all of the meta-information preserved (see the following subsection). Although the other datasets which we currently cannot improve are also valuable, our goal is to obtain new information from comparison of community detection results and ground-truth communities, rather than simply recovering the ground-truth communities.

Constructing the DBLP15 dataset
DBLP is an online reference for bibliographic information on major computer science publications [38]. As of June 17, 2015, DBLP has indexed 4316 conferences, 1417 journals, and 1,573,969 authors [39]. The whole DBLP dataset is provided in a well-formatted XML file. The snapshot/release version of the data we use can be accessed at http:// dblp.dagstuhl.de/xml/release/dblp-2015-06-02.xml.gz. The structure of this XML file is illustrated in Fig. 3. The root element is the dblp element. We call the children of the root elements Level 1 elements and the children of Level 1 elements Level 2 elements, and so on. Level 1 elements represent the individual data records [40], such as article and book, etc. Since publication-venue relation makes more sense for the journal and conference papers, and these two types of publications occupy most of DBLP, we consider only article and inproceedings elements when constructing our dataset. Level 2 elements contain the meta-information about the publications, such as title, authors, journal/proceeding names, etc.
Our goal is to obtain a coauthorship network and ground-truth information (venueauthor relation) from the XML file. Although the XML file is highly structured, such a task is still not straightforward due to the ambiguity of entities, such as conflicts or changes of author names, various abbreviations, or even journal name change. DBLP resolves the author ambiguity issue by means of a unique number for each author. However, the venue ambiguity is still an issue in DBLP: there are no unique identifiers for venues. Fortunately, each record in DBLP has a unique key, and most paper keys contain the venue information as follows: However, there are still a few exceptions. To examine the validity of venue identifiers efficiently, we manually examine the identifiers not listed in the journal and conference index provided by the DBLP website, since such indices seem to be maintained by humans and assumed to be reliable. Using this process, we found 5240 unique venues (journals or conferences). Now unique identifiers for both authors and venues make extracting the network and community information very reasonable. The next step is to create a node for each author, and create a link between two authors if they have ever coauthored in the same publication. For community information, each venue is a community, and an author belongs to a community if he/she has published in the corresponding venue.
A few authors do not have any coauthor in the DBLP database, and become isolated nodes in the generated network. Thus, we remove these authors. However, after removing those authors, some venues/communities become empty because all of their authors are removed. Hence, we remove those empty communities. After this cleaning, we obtained 1,509,944 authors in 5147 communities (venues).
This cleaned network has 51,328 (weakly) connected components, where the largest connected component contains 1,357,781 nodes, which makes 89.9% of all nodes. The remaining 51,327 connected components are all small, the largest of which has only 37 nodes. We take the largest connected component as the network to study. By  extracting the largest connected component, we obtain a network with 1,357,781 nodes,  6,369,212 edges and 5146 ground-truth communities. The ground-truth communities were divided into connected components, obtaining 93,824 communities. The divided ground-truth communities were used for comparison with detected communities.

Experiment results
We run our experiments on a server with two Intel E5-2620 processors, each having six cores, and 377-GB memory. The results are listed in Tables 2, 3, 4, 5, 6, 7, 8, and 9.
In the "internal measures" table, "coverage" measures the percentage of nodes which are assigned to at least one community; "algorithm time" and "total time" provide the runtime information. We list two measures of runtime since our algorithm (and also NISE) implemented in MATLAB directly uses a processed matrix in memory as its input. Other algorithms must first read the graph stored as an edge list or an adjacency list and convert the graph to the appropriate internal representation. Therefore, we use "algorithm time" to measure the algorithm runtime without the time needed for reading and converting the graph, which is reported by the algorithms themselves. The "total time" is the wall clock time for running the algorithm, including the time for reading  and converting the graph, which is measured with an external timer. BigClam reports its algorithm time as the sum of time used in each core, and therefore, the results are not comparable. For completeness, we added the data-loading and data-preprocessing times, which are measured separately, to obtain a "total time" for the MATLAB algorithms (our algorithm and NISE).   The "number of clusters" in the "internal measures" table are different across different methods due to the following reasons. The SCD algorithm does not provide an interface for specifying the number of communities to detect, and instead detects the number of communities automatically. For other algorithms, we specify the number of communities to detect as 5000. The actual number of communities generated by HierSymNMF2   is usually less than 5000, as discussed in "Methods for comparison" section. Also, the number of communities generated by NISE are usually a little larger than 5000, which is also an expected behavior [3]. In the "external measures" table, the "reverse precision" and "reverse recall" refer to the scores computed as if the ground-truth communities are treated as detected communities and the detected communities are treated as the ground truth, respectively. Note that the number of clusters in "external measures" is less than the one in "internal measures" due to the removal of nodes that do not appear in the ground truth.
We have the following observations from the experimental results: (1) Our HierSym-NMF2 algorithm has significant advantages over other methods in average normalized cut on most datasets except the Amazon dataset. On the Amazon dataset, HierSym-NMF2 achieves much lower average normalized cut than SCD and BigClam, and the variant h2-ngd-a obtained comparable average normalized cut (0.1491) versus Graclus (0.1450), which is not as good as NISE (0.1118). (2) HierSymNMF2 runs slower than most other algorithms on DBLP06 and DBLP15. On the Youtube dataset, HierSym-NMF2 runs faster than BigClam, Graclus and NISE. On the Amazon dataset, Hier-SymNMF2 runs faster than NISE, but slower than other methods. (3) HierSymNMF2 achieves better F1 score than BigClam and NISE on all the datasets we used. Graclus has better F1 score than HierSymNMF2 on DBLP06, Amazon, and Youtube datasets but obtained an unusually low F1 score on the DBLP15 dataset. SCD achieves higher F1 scores than HierSymNMF2. However, SCD often discovers a significantly larger number of (nonoverlapping) communities than expected and has very unbalanced precision and recall scores compared with other algorithms. The SCD algorithm finds the number of communities as it finds the communities and the number of communities cannot be given to SCD as an input. The SCD algorithm starts by assigning an initial partitioning of the graph heuristically. In short, in the initial partitioning, each node and all its neighbors form a community, and special care is taken to ensure that no node belongs to more than one community. As a result, this initial step often creates many more communities than the optimal number, although later refining procedures may reduce the number of communities. As can be seen from the experimental results, when compared with Big-Clam, Graclus, NISE, and our proposed algorithms that take the number of communities as an input, a much larger number of communities that the SCD generates does not necessarily translate to a better overall community detection result in terms of either normalized cut or F1 scores.

Conclusions and discussion
Overall, HierSymNMF2 is an effective community detection method that optimizes the average normalized cut very well, although it is not the fastest method.
To address the quality issue of the ground-truth networks, we constructed a more complete and recent version of the DBLP dataset, where most nodes have at least one community membership and also the size of the data is significantly larger.
We partially answered the three questions we raised in "Problem definition" section. For Q1 and Q2, we used measurements that are commonly used in the current literature. One has to be careful when choosing community detection methods based on external measures such as F1 score because of the incompleteness of ground-truth communities and the difference between linkage-based communities (as detected) and functional communities (as in the ground truth). Therefore, it is important to use an internal measure such as the average normalized cut to evaluate an algorithm. In addition, we believe the current evaluation methods for large-scale community detection have certain limitations. These methods are mainly based on some quality scores, e.g., the F1 score. Such quality scores can be used to compare various algorithms. However, they do not provide much more information regarding the quality of results than the average performance. We assert that to better understand a community detection result, it is necessary to develop more comprehensive evaluation methods.
For Q3, we think that the large scale of the network and the large number of communities make the community detection results hard to interpret. One way to understand the community structure better would be to develop better methods for community visualization.
As a result of the growing popularity and utility of social media communications and other channels of communications between people and groups of people, there are vast amounts of data that contain latent community information. The amount of information is overwhelming and very demanding of our current technological capabilities, which may adversley impact the ability of stakeholders to make critical and timely decisions that are important in many domains such as natural disasters, local conflicts, healthcare, and law enforcement, to name a few. These domains typically involve groups of individuals with often hidden links. Thus, it is incumbent on the research community to develop fast and effective methods to first discover the communities formed by these links and then formulate useful summaries of the information provided by the algorithms and measures in order for decision makers to initiate appropriate actions as required.
Authors' contributions DK, HP, and BD initiated the idea of the algorithm and problem formulation. RD proposed and implemented the splitting criteria, and constructed the DBLP15 dataset. DK designed and implemented the acceleration scheme and the algorithm framework. RD and DK conducted the experiments. RD wrote the initial draft of the paper and DK wrote the "Implementation" section. BD and HP rewrote and added revised and new content. All four authors iterated to finalize the manuscript. The work is supervised by HP. All authors read and approved the final manuscript.

Competing interests
The authors declare that they have no competing interests.

Availability of data and materials
The six SNAP datasets listed in Table 1 are available at https://snap.stanford.edu/data/#communities. The DBLP15 dataset is available at https://github.com/smallk/smallk_data/tree/master/dblp_ground_truth.

Funding
The work of the authors was supported in part by the National Science Foundation (NSF) Grant IIS-1348152 and the Defense Advanced Research Projects Agency (DARPA) XDATA program Grant FA8750-12-2-0309. Any opinions, findings, and conclusions or recommendations expressed in this article are those of the authors and do not necessarily reflect the views of the NSF or DARPA.
We mentioned that the SymNMF formulation is a relaxation of optimization problems related to graph partitioning. Specifically, let S G be the adjacency matrix of graph G. When S = S G , Eq. (2) is a relaxation of maximizing the ratio association. When S = D −1/2 S G D −1/2 , where D = diag(d 1 , . . . , d n ) and d i = n j=1 S G ij is the degree of node i, Eq. (2) is a relaxation of minimizing the normalized cut.
We defined normalized cut of a graph partition (A 1 , . . . , A k ) as ncut(A 1 , . . . , A k ) (Eq. (9)) using the concept of within(A i ) (Eq. (10)), and out(A i ) (Eq. (11)). With the same notations, the ratio association of a graph partition is defined as where |A i | is the number of nodes in partition A i .
In the following sections, we will use the convenient Iverson bracket [41]: where P is any statement. Also, G = (V , E) is the graph under study, where V = {v 1 , . . . , v n }. The matrix S = (w ij ) is the adjacency matrix of graph G, where w ij = w(v i , v j ) is the weight of edge (v i , v j ). The tuple (A 1 , . . . , A j ) is a partition of graph G.

Relation between SymNMF and ratio association
We rewrite Eq.
and which means H T H = I. Therefore [42], If the restriction (18) is relaxed using H ≥ 0 , i.e., nonnegative H, we will arrive at our SymNMF formulation.

Relation between SymNMF and normalized cut
Let D = diag(d 1 , . . . , d n ) , where d i = n j=1 w ij is the degree of node i and let H = (h ij ) n×k , where Then, we have tr(H T SH ) = 1 ≤ i ≤ k 1 ≤ j, l ≤ n h ji w jl h li and which means H T H = I. Similar to (19), we have When the restriction (20) is relaxed to H ≥ 0, our SymNMF formulation is obtained.

Appendix B: Normalized cut increases when a community is split into two communities
In Fig. 1, the community A 3 was split into B 1 and B 2 , and the associated increase of normalized cut is  Note that this proof does not say that more communities always correspond to larger normalized cut in general (i.e., when the communities are not obtained through recursive splitting).