Online Tensor Methods for Learning Latent Variable Models

Mixed Membership Stochastic Blockmodel, topic modeling, tensor method, stochastic gradient descent, parallel implementation, large datasets.

Introduction

The spectral or moment-based approach involves decomposition of certain empirical moment tensors, estimated from observed data to obtain the parameters of the proposed probabilistic model. Unsupervised learning for a wide range of latent variable models can be carried out efficiently via tensor-based techniques with low sample and computational complexities (Anandkumar et al., 2012). In contrast, usual methods employed in practice such as expectation maximization (EM) and variational Bayes do not have such consistency guarantees. While the previous works (Anandkumar et al., 2013b) focused on theoretical guarantees, in this paper, we focus on the implementation of the tensor methods, study its performance on several datasets.

We consider two problems: (1) community detection (wherein we compute the decomposition of a tensor which relates to the count of $3$-stars in a graph) and (2) topic modeling (wherein we consider the tensor related to co-occurrence of triplets of words in documents); decomposition of the these tensors allows us to learn the hidden communities and topics from observed data.

We recover hidden communities in several real datasets with high accuracy. When ground-truth communities are available, we propose a new error score based on the hypothesis testing methodology involving $p$-values and false discovery rates (Strimmer, 2008) to validate our results. The use of $p$-values eliminates the need to carefully tune the number of communities output by our algorithm, and hence, we obtain a flexible trade-off between the fraction of communities recovered and their estimation accuracy. We find that our method has very good accuracy on a range of network datasets: Facebook, Yelp and DBLP. We summarize the datasets used in this paper in Table 6. To get an idea of our running times, let us consider the larger DBLP collaborative dataset for a moment. It consists of $16$ million edges, one million nodes and $250$ communities. We obtain an error of $10\%$ and the method runs in about two minutes, excluding the $80$ minutes taken to read the edge data from files stored on the hard disk and converting it to sparse matrix format.

Compared to the state-of-the-art method for learning MMSB models using the stochastic variational inference algorithm of (Gopalan et al., 2012), we obtain several orders of magnitude speed-up in the running time on multiple real datasets. This is because our method consists of efficient matrix operations which are embarrassingly parallel. Matrix operations are carried out in the sparse format which is efficient especially for social network settings involving large sparse graphs. Moreover, our code is flexible to run on a range of graphs such as directed, undirected and bipartite graphs, while the code of (Gopalan et al., 2012) is designed for homophilic networks, and cannot handle bipartite graphs in its present format. Note that bipartite networks occur in the recommendation setting such as the Yelp dataset. Additionally, the variational implementation in (Gopalan et al., 2012) assumes a homogeneous connectivity model, where any pair of communities connect with the same probability and the probability of intra-community connectivity is also fixed. Our framework does not suffer from this restriction. We also provide arguments to show that the Normalized Mutual Information (NMI) and other scores, previously used for evaluating the recovery of overlapping community, can underestimate the errors.

We also employ the tensor method for topic-modeling, and there are many similarities between the topic and community settings. For instance, each document has multiple topics, while in the network setting, each node has membership in multiple communities. The words in a document are generated based on the latent topics in the document, and similarly, edges are generated based on the community memberships of the node pairs. The tensor method is even faster for topic modeling, since the word vocabulary size is typically much smaller than the size of real-world networks. We learn interesting hidden topics in New York Times corpus from UCI bag-of-words datasethttps://archive.ics.uci.edu/ml/datasets/Bag+of+Words with around $100,000$ words and $300,000$ documents in about two minutes. We present the important words for recovered topics, as well as interpret “bridging” words, which occur in many topics.

We present two implementations, viz., a GPU-based implementation which exploits the parallelism of SIMD architectures and a CPU-based implementation for larger datasets, where the GPU memory does not suffice. We discuss various aspects involved such as implicit manipulation of tensors since explicitly forming tensors would be unwieldy for large networks, optimizing for communication bottlenecks in a parallel deployment, the need for sparse matrix and vector operations since real world networks tend to be sparse, and a careful statistical approach to validating the results, when ground truth is available.

2 Related work

This paper builds on the recent works of Anandkumar et al (Anandkumar et al., 2012, 2013b) which establishes the correctness of tensor-based approaches for learning MMSB (Airoldi et al., 2008) models and other latent variable models. While, the earlier works provided a theoretical analysis of the method, the current paper considers a careful implementation of the method. Moreover, there are a number of algorithmic improvements in this paper. For instance, while (Anandkumar et al., 2012, 2013b) consider tensor power iterations, based on batch data and deflations performed serially, here, we adopt a stochastic gradient descent approach for tensor decomposition, which provides the flexibility to trade-off sub-sampling with accuracy. Moreover, we use randomized methods for dimensionality reduction in the preprocessing stage of our method which enables us to scale our method to graphs with millions of nodes.

There are other known methods for learning the stochastic block model based on techniques such as spectral clustering (McSherry, 2001) and convex optimization (Chen et al., 2012). However, these methods are not applicable for learning overlapping communities. We note that learning the mixed membership model can be reduced to a matrix factorization problem (Zhang and Yeung, 2012). While collaborative filtering techniques such as (Mnih and Salakhutdinov, 2007; Salakhutdinov and Mnih, 2008) focus on matrix factorization and the prediction accuracy of recommendations on an unseen test set, we recover the underlying latent communities, which helps with the interpretability and the statistical model can be employed for other tasks.

Although there have been other fast implementations for community detection before (Soman and Narang, 2011; Lancichinetti and Fortunato, 2009), these methods are not statistical and do not yield descriptive statistics such as bridging nodes (Nepusz et al., 2008), and cannot perform predictive tasks such as link classification which are the main strengths of the MMSB model. With the implementation of our tensor-based approach, we record huge speed-ups compared to existing approaches for learning the MMSB model.

To the best of our knowledge, while stochastic methods for matrix decomposition have been considered earlier (Oja and Karhunen, 1985; Arora et al., 2012), this is the first work incorporating stochastic optimization for tensor decomposition, and paves the way for further investigation on many theoretical and practical issues. We also note that we never explicitly form or store the subgraph count tensor, of size $O(n^{3})$ where $n$ is the number of nodes, in our implementation, but directly manipulate the neighborhood vectors to obtain tensor decompositions through stochastic updates. This is a crucial departure from other works on tensor decompositions on GPUs (Ballard et al., 2011; Schatz et al., 2013), where the tensor needs to be stored and manipulated directly.

Tensor Forms for Topic and Community Models

In this section, we briefly recap the topic and community models, as well as the tensor forms for their exact moments, derived in (Anandkumar et al., 2012, 2013b).

The Dirichlet distribution allows us to specify the extent of overlap among the topics by controlling for sparsity in topic density function. A larger $\alpha_{0}$ results in more overlapped (mixed) topics. A special case of $\alpha_{0}=0$ is the single topic model.

We consider the first three order empirical moments, given by

We recall Theorem 3.5 of (Anandkumar et al., 2012):

where $\mu=[\mu_{1},\ldots,\mu_{k}]$ and $\mu_{i}=\Pr\left(x_{t}|h=i\right)$, $\forall t\in[l]$. In other words, $\mu$ is the topic-word matrix.

From the Lemma 1, we observe that the first three moments of a LDA topic model have a simple form involving the topic-word matrix $\mu$ and Dirichlet parameters $\alpha_{i}$. In (Anandkumar et al., 2012), it is shown that these parameters can be recovered under a weak non-degeneracy assumption. We will employ tensor decomposition techniques to learn the parameters.

2 Mixed Membership Model

In the mixed membership stochastic block model (MMSB), introduced by (Airoldi et al., 2008), the edges in a social network are related to the hidden communities of the nodes. A batch tensor decomposition technique for learning MMSB was derived in (Anandkumar et al., 2013b).

The community connectivity matrix is denoted by $P\in^{k\times k}$ where $P(a,b)$ measures the connectivity between communities $a$ and $b$, $\forall a,b\in[k]$. We model the adjacency matrix entries as either of the two settings given below:

This models a network with unweighted edges. It is used for Facebook and DBLP datasets in Section 6 in our experiments.

This models a network with weighted edges. It is used for the Yelp dataset in Section 6 to incorporate the review ratings.

The tensor decomposition approach involves up to third order moments, computed from the observed network. In order to compute the moments, we partition the nodes randomly into sets $X,A,B,C$. Let $F_{A}:=\Pi_{A}^{\top}P^{\top}$, $F_{B}:=\Pi_{B}^{\top}P^{\top}$, $F_{C}:=\Pi_{C}^{\top}P^{\top}$ (where $P$ is the community connectivity matrix and $\Pi$ is the membership matrix) and $\hat{\alpha}:=\left(\frac{\alpha_{1}}{\alpha_{0}},\ldots,\frac{\alpha_{k}}{\alpha_{0}}\right)$ denote the normalized Dirichlet concentration parameter. We define pairs over $Y_{1}$ and $Y_{2}$ as $\operatorname{Pairs}(Y_{1},Y_{2}):=G_{X,Y_{1}}^{\top}\otimes G_{X,Y_{2}}^{\top}$. Define the following matrices

We consider the first three empirical moments, given by

We now recap Proposition 2.2 of (Anandkumar et al., 2013a) which provides the form of these moments under expectation.

where $\otimes$ denotes the Kronecker product and $(F_{A})_{i}$ corresponds to the $i^{th}$ column of $F_{A}$.

We observe that the moment forms above for the MMSB model have a similar form as the moments of the topic model in the previous section. Thus, we can employ a unified framework for both topic and community modeling involving decomposition of the third order moment tensors $M_{3}^{\operatorname{Top}}$ and $M_{3}^{\operatorname{Com}}$. Second order moments $M_{2}^{\operatorname{Top}}$ and $M_{2}^{\operatorname{Com}}$ are used for preprocessing of the data (i.e., whitening, which is introduced in detail in Section 3.1). For the sake of the simplicity of the notation, in the rest of the paper, we will use $M_{2}$ to denote empirical second order moments for both $M_{2}^{\operatorname{Top}}$ in topic modeling setting, and $M_{2}^{\operatorname{Com}}$ in the mixed membership model setting. Similarly, we will use $M_{3}$ to denote empirical third order moments for both $M_{3}^{\operatorname{Top}}$ and $M_{3}^{\operatorname{Com}}$.

Learning using Third Order Moment

Our learning algorithm uses up to the third-order moment to estimate the topic word matrix $\mu$ or the community membership matrix $\Pi$. First, we obtain co-occurrence of triplet words or subgraph counts (implicitly). Then, we perform preprocessing using second order moment $M_{2}$. Then we perform tensor decomposition efficiently using stochastic gradient descent (Kushner and Yin, 2003) on $M_{3}$. We note that, in our implementation of the algorithm on the Graphics Processing Unit (GPU), linear algebraic operations are extremely fast. We also implement our algorithm on the CPU for large datasets which exceed the memory capacity of GPU and use sparse matrix operations which results in large gains in terms of both the memory and the running time requirements. The overall approach is summarized in Algorithm 1.

The whitening matrix $W$ is computed via truncated $k-$svd of the second order moments.

where $U_{M_{2}}$ and $\Sigma_{M_{2}}=\operatorname{diag}(\sigma_{M_{2},1},\ldots,\sigma_{M_{2},k})$ are the top $k$ singular vectors and singular values of $M_{2}$ respectively. We then perform multilinear transformations on the triplet data using the whitening matrix. The whitened data is thus

2 Stochastic Tensor Gradient Descent

In (Anandkumar et al., 2013b) and (Anandkumar et al., 2012), the power method with deflation is used for tensor decomposition where the eigenvectors are recovered by iterating over multiple loops in a serial manner. Furthermore, batch data is used in their iterative power method which makes that algorithm slower than its stochastic counterpart. In addition to implementing a stochastic spectral optimization algorithm, we achieve further speed-up by efficiently parallelizing the stochastic updates.

where $t\in X$ and denotes the index of the online data and $\bar{y}_{A}$, $\bar{y}_{B}$, and $\bar{y}_{C}$ denote the mean of the whitened data. Our goal is to find a symmetric CP decomposition of the whitened tensor.

where $v_{i}$ are the unknown components to be estimated, and $\theta>0$ is some fixed parameter.

In order to encourage orthogonality between eigenvectors, we have the extra term as $\theta\lVert\sum_{i\in[k]}\otimes^{3}v_{i}\lVert_{F}^{2}$. Since $\lVert\sum_{t\in X}{\mathcal{T}^{t}}\lVert_{F}^{2}$ is a constant, the above minimization is the same as minimizing a loss function $L(\mathbf{v}):=\frac{1}{n_{\mbox{\tiny\!X}}}\sum_{t}L^{t}(\mathbf{v})$, where $L^{t}(\mathbf{v})$ is the loss function evaluated at node $t\in X$, and is given by

The loss function has two terms, viz., the term $\lVert\sum_{i\in[k]}\otimes^{3}v_{i}\lVert_{F}^{2}$, which can be interpreted as the orthogonality cost, which we need to minimize, and the second term $\langle\sum_{i\in[k]}{\otimes^{3}v_{i}},{\mathcal{T}^{t}}\rangle$, which can be viewed as the correlation reward to be maximized. The parameter $\theta$ provides additional flexibility for tuning between the two terms.

where $\beta^{t}$ is the learning rate. Computing the derivative of the loss function and substituting the result leads to the following lemma.

The stochastic updates for the eigenvectors are given by

In Equation (4), all our tensor operations are in terms of efficient sample vector inner products, and no tensor is explicitly formed. The multilinear operations are shown in Figure 1. We choose $\theta=1$ in our experiments to ensure that there is sufficient penalty for non-orthogonality, which prevents us from obtaining degenerate solutions.

After learning the decomposition of the third order moment, we perform post-processing to estimate $\widehat{\Pi}$.

3 Post-processing

Eigenvalues $\Lambda:=[\lambda_{1},\lambda_{2},\ldots,\lambda_{k}]$ are estimated as the norm of the eigenvectors $\lambda_{i}={\lVert{\phi_{i}}\rVert}^{3}$.

After we obtain $\Lambda$ and $\Phi$, the estimate for the topic-word matrix is given by

and in the community setting, the community membership matrix is given by

where $A^{c}:=X\cup B\cup C$. Similarly, we estimate $\hat{\Pi}_{A}$ by exchanging the roles of $X$ and $A$. Next, we obtain the Dirichlet distribution parameters

where $\gamma^{2}$ is chosen such that we have normalization$\sum_{i\in[k]}\hat{\alpha}_{i}:=\sum_{i\in[k]}\frac{\alpha_{i}}{\alpha_{0}}=1.$

Thus, we perform STGD method to estimate the eigenvectors and eigenvalues of the whitened tensor, and then use these to estimate the topic word matrix $\mu$ and community membership matrix $\widehat{\Pi}$ by thresholding.

Implementation Details

Note that for the topic model, the second order moment $M_{2}$ can be computed easily from the word-frequency vector. On the other hand, for the community setting, computing $M_{2}$ requires additional linear algebraic operations. It requires computation of matrices $Z_{B}$ and $Z_{C}$ in equation (7). This requires computation of pseudo-inverses of “Pairs” matrices. Now, note that pseudo-inverse of $\left(\operatorname{Pairs}\left(B,C\right)\right)$ in Equation (7) can be computed using rank $k$-SVD:

We then orthogonalize the third order moments to reduce the dimension of its modes to $k$. We perform linear transformations on the data corresponding to the partitions $A$, $B$ and $C$ using the whitening matrix. The whitened data is thus $y^{t}_{A}:=\left<W,G^{\top}_{t,A}\right>$, $y^{t}_{B}:=\left<W,Z_{B}G^{\top}_{t,B}\right>$, and $y^{t}_{C}:=\left<W,Z_{C}G^{\top}_{t,C}\right>$, where $t\in X$ and denotes the index of the online data. Since $k\ll n$, the dimensionality reduction is crucial for our speedup.

2 Efficient Randomized SVD Computations

When we consider very large-scale data, the whitening matrix is a bottleneck to handle when we aim for fast running times. We obtain the low rank approximation of matrices using random projections. In the CPU implementation, we use tall-thin SVD (on a sparse matrix) via the Lanczos algorithm after the projection and in the GPU implementation, we use tall-thin QR. We give the overview of these methods below. Again, we use graph community membership model without loss of generality.

Tall-thin SVD: This is used in the CPU implementation. The whitening matrix can be obtained by

The pseudo code for computing the whitening matrix $W$ using tall-thin SVD is given in Algorithm 2.

where $\Psi=S^{\top}\left(\operatorname{Pairs}\left(B,C\right)\right)S$ and $J=\left(\operatorname{Pairs}\left(B,C\right)\right)S$. The pseudo code for computing pseudoinverses is given in Algorithm 3.

The sparse representation of the data allows for scalability on a single machine to datasets having millions of nodes. Although the GPU has SIMD architecture which makes parallelization efficient, it lacks advanced libraries with sparse SVD operations and out-of-GPU-core implementations. We therefore implement the sparse format on CPU for sparse datasets. We implement our algorithm using random projection for efficient dimensionality reduction (Clarkson and Woodruff, 2012) along with the sparse matrix operations available in the Eigen toolkithttp://eigen.tuxfamily.org/index.php?title=Main_Page, and we use the SVDLIBC (Berry et al., 2002) library to compute sparse SVD via the Lanczos algorithm. Theoretically, the Lanczos algorithm (Golub and Van Loan, 2013) on a $n\times n$ matrix takes around $(2d+8)n$ flops for a single step where $d$ is the average number of non-zero entries per row.

Tall-thin QR: This is used in the GPU implementation due to the lack of library to do sparse tall-thin SVD. The difference is that we instead implement a tall-thin QR on $O$, therefore the whitening matrix is obtained as

The main bottleneck for our GPU implementation is device storage, since GPU memory is highly limited and not expandable. Random projections help in reducing the dimensionality from $O(n\times n)$ to $O(n\times k)$ and hence, this fits the data in the GPU memory better. Consequently, after the whitening step, we project the data into $k$-dimensional space. Therefore, the STGD step is dependent only on $k$, and hence can be fit in the GPU memory. So, the main bottleneck is computation of large SVDs. In order to support larger datasets such as the DBLP dataset which exceed the GPU memory capacity, we extend our implementation with out-of-GPU-core matrix operations and the Nystrom method (Gittens and Mahoney, 2013) for the whitening matrix computation and the pseudoinverse computation in the pre-processing module.

3 Stochastic updates

STGD can potentially be the most computationally intensive task if carried out naively since the storage and manipulation of a $O(n^{3})$-sized tensor makes the method not scalable. However we overcome this problem since we never form the tensor explicitly; instead, we collapse the tensor modes implicitly as shown in Figure 1. We gain large speed up by optimizing the implementation of STGD.To implement the tensor operations efficiently we convert them into matrix and vector operations so that they are implemented using BLAS routines. We obtain whitened vectors $y_{A},y_{B}$ and $y_{C}$ and manipulate these vectors efficiently to obtain tensor eigenvector updates using the gradient scaled by a suitable learning rate.

In STGD, note that the storage needed for the iterative part does not depend on the number of nodes in the dataset, rather, it depends on the parameter $k$, i.e., the number of communities to be estimated, since whitening performed before STGD leads to dimensionality reduction. This makes it suitable for storing the required buffers in the GPU memory, and using the CULA device interface for the BLAS operations. In Figure 3, we illustrate the data transfer involved in the GPU standard and device interface codes. While the standard interface involves data transfer (including whitened neighborhood vectors and the eigenvectors) at each stochastic iteration between the CPU memory and the GPU memory, the device interface involves allocating and retaining the eigenvectors at each stochastic iteration which in turn speeds up the spectral estimation.

We compare the running time of the CULA device code with the MATLAB code (using the tensor toolbox (Bader et al., 2012)), CULA standard code and Eigen sparse code in Figure 4. As expected, the GPU implementations of matrix operations are much faster and scale much better than the CPU implementations. Among the CPU codes, we notice that sparsity and optimization offered by the Eigen toolkit gives us huge gains. We obtain orders of magnitude of speed up for the GPU device code as we place the buffers in the GPU memory and transfer minimal amount of data involving the whitened vectors only once at the beginning of each iteration. The running time for the CULA standard code is more than the device code because of the CPU-GPU data transfer overhead. For the same reason, the sparse CPU implementation, by avoiding the data transfer overhead, performs better than the GPU standard code for very small number of communities. We note that there is no performance degradation due to the parallelization of the matrix operations. After whitening, the STGD requires the most code design and optimization effort, and so we convert that into BLAS-like routines.

4 Computational Complexity

We partition the execution of our algorithm into three main modules namely, pre-processing, STGD and post-processing, whose various matrix operation counts are listed above in Table 1.

The theoretical asymptotic complexity of our method is summarized in Table 2 and is best addressed by considering the parallel model of computation (JáJá, 1992), i.e., wherein a number of processors or compute cores are operating on the data simultaneously in parallel. This is justified considering that we implement our method on GPUs and matrix products are embarrassingly parallel. Note that this is different from serial computational complexity. We now break down the entries in Table 2. First, we recall a basic lemma regarding the lower bound on the time complexity for parallel addition along with the required number of cores to achieve a speed-up.

(JáJá, 1992) Addition of $s$ numbers in serial takes $O(s)$ time; with $\Omega(s/\log s)$ cores, this can be improved to $O(\log s)$ time in the best case.

Essentially, this speed-up is achieved by recursively adding pairs of numbers in parallel.

Lemma 7 follows by simply parallelizing the sparse inner products and applying Lemma 6 for the addition in the inner products. Note that, this can be generalized to the fact that given $c$ cores, the multiplication can be performed in $O(\max(psr/c,\log s))$ running time.

In preprocessing, given $c$ compute cores, we first do random projection using matrix multiplication. We multiply an $O(n)\times O(n)$ matrix $M_{2}$ with an $O(n)\times O(k)$ random matrix $S$. Therefore, this requires $O(nsk)$ serial operations, where $s$ is the number of non-zero elements per row/column of $M_{2}$. Using Lemma 7, given $c=\frac{nsk}{\log s}$ cores, we could achieve $O(\log s)$ computational complexity. However, the parallel computational complexity is not further reduced with more than $\frac{nsk}{\log s}$ cores.

After the multiplication, we use tall-thin SVD for CPU implementation, and tall-thin QR for GPU implementation.

We perform Lanczos SVD on the tall-thin sparse $O(n)\times O(k)$ matrix, which involves a tri-diagonalization followed with the QR on the tri-diagonal matrix. Given $c=\frac{nsk}{\log s}$ cores, the computational complexity of the tri-diagonalization is $O(\log s)$. We then do QR on the tridiagonal matrix which is as cheap as $O(k^{2})$ serially. Each orthogolization requires $O(k)$ inner products of constant entry vectors, and there are $O(k)$ such orthogolizations to be done. Therefore given $O(k)$ cores, the complexity is $O(k)$. More cores does not help since the degree of parallelism is $k$.

Alternatively, we perform QR in the GPU implementation which takes $O(sk^{2})$. To arrive at the complexity of obtaining $Q$, we analyze the Gram-Schmidt orthonormalization procedure under sparsity and parallelism conditions. Consider a serial Gram-Schmidt on $k$ columns (which are $s$-dense) of $O(n)\times O(k)$ matrix. For each of the columns $2$ to $k$, we perform projection on the previously computed components and subtract it. Both inner product and subtraction operations are on the $s$-dense columns and there are $O(s)$ operations which are done $O(k^{2})$ times serially. The last step is the normalization of $k$ $s$-dense vectors with is an $O(sk)$ operation. This leads to a serial complexity of $O(sk^{2}+sk)=O(sk^{2})$. Using this, we may obtain the parallel complexity in different regimes of the number of cores as follows.

Parallelism for inner products : For each component $i$, we need $i-1$ projections on previous components which can be parallel. Each projection involves scaling and inner product operations on a pair of $s$-dense vectors. Using Lemma 6, projection for component $i$ can be performed in $O(\max(\frac{sk}{c},\log s))$ time. $O(\log s)$ complexity is obtained using $O(sk/\log s)$ cores.

Parallelism for subtractions: For each component $i$, we need $i-1$ subtractions on a $s$-dense vector after the projection. Serially the subtraction requires $O(sk)$ operations, and this can be reduced to $O(\log k)$ with $O(sk/\log k)$ cores in the best case. The complexity is $O(\max(\frac{sk}{c},\log k))$.

Combing the inner products and subtractions, the complexity is $O\left(\max(\frac{sk}{c},\log s)+\max(\frac{sk}{c},\log k)\right)$ for component $i$. There are $k$ components in total, which can not be parallel. In total, the complexity for the parallel QR is $O\left(\max(\frac{sk^{2}}{c},\log s)+\max(\frac{sk^{2}}{c},\log k)\right)$.

The serial time complexity of SVD is $O(n^{2}k)$ but with randomized dimensionality reduction (Gittens and Mahoney, 2013) and parallelization (Constantine and Gleich, 2011), this is significantly reduced.

4.2 STGD

In STGD, we perform implicit stochastic updates, consisting of a constant number of matrix-matrix and matrix-vector products, on the set of eigenvectors and whitened samples which is of size $k\times k$. When $c\in[1,k^{3}/\log k]$, we obtain a running time of $O({k^{3}/c})$ for computing inner products in parallel with $c$ compute cores since each core can perform an inner product to compute an element in the resulting matrix independent of other cores in linear time. For $c\in(k^{3}/\log k,\infty]$, using Lemma 6, we obtain a running time of $O(\log k)$. Note that the STGD time complexity is calculated per iteration.

4.3 Post-processing

Finally, post-processing consists of sparse matrix products as well. Similar to pre-processing, this consists of multiplications involving the sparse matrices. Given $s$ number of non-zeros per column of an $O(n)\times O(k)$ matrix, the effective number of elements reduces to $O(sk)$. Hence, given $c\in[1,nks/\log s]$ cores, we need $O({nsk/c})$ time to perform the inner products for each entry of the resultant matrix. For $c\in(nks/\log s,\infty]$, using Lemma 6, we obtain a running time of $O(\log s)$.

Note that $nk^{2}$ is the complexity of computing the exact SVD and we reduce it to $O(k)$ when there are sufficient cores available. This is meant for the setting where $k$ is small. This $k^{3}$ complexity of SVD on $O(k\times k)$ matrix can be reduced to $O(k)$ using distributed SVD algorithms e.g. (Kannan et al., 2014; Feldman et al., 2013). We note that the variational inference algorithm complexity, by Gopalan and Blei (Gopalan and Blei, 2013), is $O(mk)$ for each iteration, where $m$ denotes the number of edges in the graph, and $n<m<n^{2}$. In the regime that $n\gg k$, our algorithm is more efficient. Moreover, a big difference is in the scaling with respect to the size of the network and ease of parallelization of our method compared to variational one.

Validation methods

The test statistic used for the $p$-value testing of the estimated communities is

The right $p$-value is obtained via the probability of obtaining a value (say $t_{ij}$) greater than the test statistic $T_{ij}$, and it is defined as

Note that $T_{ij}$ has Student’s $t$-distribution with degree of freedom $n-2$ (i.e. $T_{ij}\sim t_{n-2}$). Thus, we obtain the right $p$-valueThe right $p$-value accounts for the fact that when two communities are anti-correlated they are not paired up. Hence note that in the special case of block model in which the estimated communities are just permuted version of the ground truth communities, the pairing results in a perfect matching accurately..

In this way, we compute the $\mathsf{\mathbf{P}_{\text{val}}}$ matrix as

2 Evaluation metrics

Validating the results requires a matching of the true membership $\Pi$ with estimated membership $\widehat{\Pi}$. Let $\mathsf{P_{\text{val}}}(\Pi_{i},\widehat{\Pi}_{j})$ denote the right $p$-value under the null hypothesis that $\Pi_{i}$ and $\widehat{\Pi}_{j}$ are statistically independent. We use the $p$-value test to find out pairs $\Pi_{i},\widehat{\Pi}_{j}$ which pass a specified $p$-value threshold, and we denote such pairs using a bipartite graph $G_{\{\mathsf{P_{\text{val}}}\}}$. Thus, $G_{\{\mathsf{P_{\text{val}}}\}}$ is defined as

and the edges of $G_{\{\mathsf{P_{\text{val}}}\}}$ satisfy

A simple example is shown in Figure 5, in which $\Pi_{2}$ has statistically significant dependence with $\widehat{\Pi}_{1}$, i.e., the probability of not rejecting the null hypothesis is small (recall that null hypothesis is that they are independent). If no estimated membership vector has a significant overlap with $\Pi_{3}$, then $\Pi_{3}$ is not recovered. There can also be multiple pairings such as for $\Pi_{1}$ and $\{\widehat{\Pi}_{2},\widehat{\Pi}_{3},\widehat{\Pi}_{6}\}$. The $p$-value test between $\Pi_{1}$ and $\{\widehat{\Pi}_{2},\widehat{\Pi}_{3},\widehat{\Pi}_{6}\}$ indicates that probability of not rejecting the null hypothesis is small, i.e., they are independent. We use $0.01$ as the threshold. The same holds for $\Pi_{2}$ and $\{\widehat{\Pi}_{1}\}$ and for $\Pi_{4}$ and $\{\widehat{\Pi}_{4},\widehat{\Pi}_{5}\}$. There can be a perfect one to one matching like for $\Pi_{2}$ and $\widehat{\Pi}_{1}$ as well as a multiple matching such as for $\Pi_{1}$ and $\{\widehat{\Pi}_{2},\widehat{\Pi}_{3},\widehat{\Pi}_{6}\}$. Or another multiple matching such as for $\{\Pi_{1},\Pi_{2}\}$ and $\widehat{\Pi}_{3}$.

Let $\operatorname{Degree}_{i}$ denote the degree of ground truth community $i\in[k]$ in $G_{\{\mathsf{P_{\text{val}}}\}}$, we define the recovery ratio as follows.

The perfect case is that all the memberships have at least one significant overlapping estimated membership, giving a recovery ratio of $100\%$.

For performance analysis of our learning algorithm, we use an error function given as follows:

where $E_{\{\mathsf{P_{\text{val}}}\}}$ denotes the set of edges based on thresholding of the $p$-values.

The error function incorporates two aspects, namely the $l_{1}$ norm error between each estimated community and the corresponding paired ground truth community, and the error induced by false pairings between the estimated and ground-truth communities through $p$-value testing. For the former $l_{1}$ norm error, we normalize with $n$ which is reasonable and results in the range of the error in $$. For the latter, we define the average error function as the summation of all paired memberships errors divided by the true number of communities$k$. In this way we penalize falsely discovered pairings by summing them up. Our error function can be greater than 1 if there are too many falsely discovered pairings through$p$-value testing (which can be as large as$k\times\widehat{k}$).

Bridgeness in overlapping communities is an interesting measure to evaluate. A bridge is defined as a vertex that crosses structural holes between discrete groups of people and bridgeness analyzes the extent to which a given vertex is shared among different communities (Nepusz et al., 2008). Formally, the bridgeness of a vertex $i$ is defined as

Note that centrality measures should be used in conjunction with bridge score to distinguish outliers from genuine bridge nodes (Nepusz et al., 2008). The degree-corrected bridgeness is used to evaluate our results and is defined as

Experimental Results

The specifications of the machine on which we run our code are given in Table 3.

We perform experiments for both the stochastic block model ($\alpha_{0}=0$) and the mixed membership model. For the mixed membership model, we set the concentration parameter $\alpha_{0}=1$. We note that the error is around $8\%-14\%$ and the running times are under a minute, when $n\leq 10000$ and $n\gg k$The code is available at https://github.com/FurongHuang/Fast-Detection-of-Overlapping-Communities-via-Online-Tensor-Methods.

We observe that more samples result in a more accurate recovery of memberships which matches intuition and theory. Overall, our learning algorithm performs better in the stochastic block model case than in the mixed membership model case although we note that the accuracy is quite high for practical purposes. Theoretically, this is expected since smaller concentration parameter $\alpha_{0}$ is easier for our algorithm to learn (Anandkumar et al., 2013b). Also, our algorithm is scalable to an order of magnitude more in $n$ as illustrated by experiments on real-world large-scale datasets.

Note that we threshold the estimated memberships to clean the results. There is a tradeoff between match ratio and average error via different thresholds. In synthetic experiments, the tradeoff is not evident since a perfect matching is always present. However, we need to carefully handle this in experiments involving real data.

We perform experiments for the bag of words dataset (Bache and Lichman, 2013) for The New York Times. We set the concentration parameter to be $\alpha_{0}=1$ and observe top recovered words in numerous topics. The results are in Table 4. Many of the results are expected. For example, the top words in topic # 11 are all related to some bad personality.

We also present the words with most spread membership, i.e., words that belong to many topics as in Table 5. As expected, we see minutes, consumer, human, member and so on. These words can appear in a lot of topics, and we expect them to connect topics.

We describe the results on real datasets summarized in Table 6 in detail below. The simulations are summarized in Table 7.

The results are presented in Table 7. We note that our method, in both dense and sparse implementations, performs very well compared to the state-of-the-art variational method. For the Yelp dataset, we have a bipartite graph where the business nodes are on one side and user nodes on the other and use the review stars as the edge weights. In this bipartite setting, the variational code provided by Gopalan et al (Gopalan et al., 2012) does not work on since it is not applicable to non-homophilic models. Our approach does not have this restriction. Note that we use our dense implementation on the GPU to run experiments with large number of communities $k$ as the device implementation is much faster in terms of running time of the STGD step.On the other hand, the sparse implementation on CPU is fast and memory efficient in the case of sparse graphs with a small number of communities while the dense implementation on GPU is faster for denser graphs such as Facebook. Note that data reading time for DBLP is around 4700 seconds, which is not negligible as compared to other datasets (usually within a few seconds). Effectively, our algorithm, excluding the file I/O time, executes within two minutes for $k=10$ and within ten minutes for $k=100$.

The ground truth on business attributes such as location and type of business are available (but not provided to our algorithm) and we provide the distribution in Figure 6 on the left side. There is also a natural trade-off between recovery ratio and average error or between attempting to recover all the business communities and the accuracy of recovery. We can either recover top significant communities with high accuracy or recover more with lower accuracy. We demonstrate the trade-off in Figure 6 on the right side.

We select the top ten categories recovered with the lowest error and report the business with highest weights in $\widehat{\Pi}$. Among the matched communities, we find the business with the highest membership weight (Table 9). We can see that most of the “top” recovered businesses are rated high. Many of the categories in the top ten list are restaurants as they have a large number of reviewers. Our method can recover restaurant category with high accuracy, and the specific restaurant in the category is a popular result (with high number of stars). Also, our method can also recover many of the categories with low review counts accurately like hobby shops, yoga, churches, galleries and religious organizations which are the “niche” categories with a dedicated set of reviewers, who mostly do not review other categories.

The top bridging nodes recovered by our method for the Yelp dataset are given in the Table 8. The bridging nodes have multiple attributes typically, the type of business and its location. In addition, the categories may also be hierarchical: within restaurants, different cuisines such as Italian, American or Pizza are recovered by our method. Moreover, restaurants which also function as bars or lounges are also recovered as top bridging nodes in our method. Thus, our method can recover multiple attributes for the businesses efficiently.

There are $5$ attributes associated with all the $11537$ businesses, which are “open”, “Categories”, “Location”, “Review Counts” and “Stars”. We model ground truth communities as a combination of “Categories” and “Location”. We select business categories with more than $20$ members and remove all businesses which are closed. $10010$ businesses are remained. Only $28588$ users are involved in reviews towards the $10010$ businesses. There are $3$ attributes associated with all the $28588$ users, which are “Female”, “Male”, “Review Counts” and “Stars”. Although we do not directly know the gender information from the dataset, a name-gender guesser https://github.com/amacinho/Name-Gender-Guesser by Amac Herdagdelen. is used to estimate gender information using names.

We provide some sample visualization results in Figure 7 for both the ground truth and the estimates from our algorithm. We sub-sample the users and businesses, group the users into male and female categories, and consider nail salon and tire businesses. Analysis of ground truth reveals that nail salon and tire businesses are very discriminative of the user genders, and thus we employ them for visualization. We note that both the nail salon and tire businesses are categorized with high accuracy, while users are categorized with poorer accuracy.

Our algorithm can also recover the attributes of users. However, the ground truth available about users is far more limited than businesses, and we only have information on gender, average review counts and average stars (we infer the gender of the users through their names). Our algorithm can recover all these attributes. We observe that gender is the hardest to recover while review counts is the easiest. We see that the other user attributes recovered by our algorithm correspond to valuable user information such as their interests, location, age, lifestyle, etc. This is useful, for instance, for businesses studying the characteristics of their users, for delivering better personalized advertisements for users, and so on.

A snapshot of the Facebook network of UNC (Traud et al., 2010) is provided with user attributes. The ground truth communities are based on user attributes given in the dataset which are not exposed to the algorithm. There are $360$ top communities with sufficient (at least 20) users. Our algorithm can recover these attributes with high accuracy; see main paper for our method’s results compared with variational inference result (Gopalan et al., 2012).

We also obtain results for a range of values of $\alpha_{0}$ (Figure 8). We observe that the recovery ratio improves with larger $\alpha_{0}$ since a larger $\alpha_{0}$ can recover overlapping communities more efficiently while the error score remains relatively the same.

For the Facebook dataset, the top ten communities recovered with lowest error consist of certain high schools, second majors and dorms/houses. We observe that high school attributes are easiest to recover and second major and dorm/house are reasonably easy to recover by looking at the friendship relations in Facebook. This is reasonable: college students from the same high school have a high probability of being friends; so do colleges students from the same dorm.

The DBLP data contains bibliographic recordshttp://dblp.uni-trier.de/xml/Dblp.xml with various publication venues, such as journals and conferences, which we model as communities. We then consider authors who have published at least one paper in a community (publication venue) as a member of it. Co-authorship is thus modeled as link in the graph in which authors are represented as nodes. In this framework, we could recover the top authors in communities and bridging authors.

Conclusion

In this paper, we presented a fast and unified moment-based framework for learning overlapping communities as well as topics in a corpus. There are several key insights involved. Firstly, our approach follows from a systematic and guaranteed learning procedure in contrast to several heuristic approaches which may not have strong statistical recovery guarantees. Secondly, though using a moment-based formulation may seem computationally expensive at first sight, implementing implicit “tensor” operations leads to significant speed-ups of the algorithm. Thirdly, employing randomized methods for spectral methods is promising in the computational domain, since the running time can then be significantly reduced.

This paper paves the way for several interesting directions for further research. While our current deployment incorporates community detection in a single graph, extensions to multi-graphs and hypergraphs are possible in principle. A careful and efficient implementation for such settings will be useful in a number of applications. It is natural to extend the deployment to even larger datasets by having cloud-based systems. The issue of efficient partitioning of data and reducing communication between the machines becomes significant there. Combining our approach with other simple community detection approaches to gain even more speedups can be explored.

Acknowledgement

The first author is supported by NSF BIGDATA IIS-1251267, the second author is supported in part by UCI graduate fellowship and NSF Award CCF-1219234, and the last author is supported in part by Microsoft Faculty Fellowship, NSF Career award CCF-1254106, NSF Award CCF-1219234, and ARO YIP Award W911NF-13-1-0084. The authors acknowledge insightful discussions with Prem Gopalan, David Mimno, David Blei, Qirong Ho, Eric Xing, Carter Butts, Blake Foster, Rui Wang, Sridhar Mahadevan, and the CULA team. Special thanks to Prem Gopalan and David Mimno for providing the variational code and answering all our questions. The authors also thank Daniel Hsu and Sham Kakade for initial discussions regarding the implementation of the tensor method. We also thank Dan Melzer for helping us with the system-related issues.

References

Appendix

A Stochastic Updates

where $x\in X$ and $t$ denotes the index of the online data.

The stochastic gradient descent algorithm is obtained by taking the derivative of the loss function $\frac{\partial L^{t}(\mathbf{v})}{\partial v_{i}}$:

for $i\in[k]$, where $y_{A}^{t}$, $y_{B}^{t}$ and $y_{C}^{t}$ are the online whitened data points as discussed in the whitening step and $\theta$ is a constant factor that we can set.

The iterative updating equation for the stochastic gradient update is given by

for $i\in[k]$, where $\beta^{t}$ is the learning rate, $\phi^{t}_{i}$ is the last iteration eigenvector and $\phi^{t}_{i}$ is the updated eigenvector. We update eigenvectors through

Now we shift the updating steps so that they correspond to the centered Dirichlet moment forms, i.e.,

B Proof of correctness of our algorithm:

We now prove the correctness of our algorithm.

Thus, our whitening matrix is computed. Now, our whitened tensor is $\mathcal{T}$ is given by

C GPU Architecture

The algorithm we propose is very amenable to parallelization and is scalable which makes it suitable to implement on processors with multiple cores in it. Our method consists of simple linear algebraic operations, thus enabling us to utilize Basic Linear Algebra Subprograms (BLAS) routines such as BLAS I (vector operations), BLAS II (matrix-vector operations), BLAS III (matrix-matrix operations), Singular Value Decomposition (SVD), and iterative operations such as stochastic gradient descent for tensor decomposition that can easily take advantage of Single Instruction Multiple Data (SIMD) hardware units present in the GPUs. As such, our method is amenable to parallelization and is ideal for GPU-based implementation.

From a higher level point of view, a typical GPU based computation is a three step process involving data transfer from CPU memory to GPU global memory, operations on the data now present in GPU memory and finally, the result transfer from the GPU memory back to the CPU memory. We use the CULA library for implementing the linear algebraic operations.

The GPUs achieve massive parallelism by having hundreds of homogeneous processing cores integrated on-chip. Massive replication of these cores provides the parallelism needed by the applications that run on the GPUs. These cores, for the Nvidia GPUs, are known as CUDA cores, where each core has fully pipelined floating-point and integer arithmetic logic units. In Nvidia’s Kepler architecture based GPUs, these CUDA cores are bunched together to form a Streaming Multiprocessor (SMX). These SMX units act as the basic building block for Nvidia Kepler GPUs. Each GPU contains multiple SMX units where each SMX unit has 192 single-precision CUDA cores, 64 double-precision units, 32 special function units, and 32 load/store units for data movement between cores and memory.

Each SMX has L$1$, shared memory and a read-only data cache that are common to all the CUDA cores in that SMX unit. Moreover, the programmer can choose between different configurations of the shared memory and L$1$ cache. Kepler GPUs also have an L$2$ cache memory of about $1.5$MB that is common to all the on-chip SMXs. Apart from the above mentioned memories, Kepler based GPU cards come with a large DRAM memory, also known as the global memory, whose size is usually in gigabytes. This global memory is also visible to all the cores. The GPU cards usually do not exist as standalone devices. Rather they are part of a CPU based system, where the CPU and GPU interact with each other via PCI (or PCI Express) bus.

In order to program these massively parallel GPUs, Nvidia provides a framework known as CUDA that enables the developers to write programs in languages like C, C++, and Fortran etc. A CUDA program constitutes of functions called CUDA kernels that execute across many parallel software threads, where each thread runs on a CUDA core. Thus the GPU’s performance and scalability is exploited by the simple partitioning of the algorithm into fixed sized blocks of parallel threads that run on hundreds of CUDA cores. The threads running on an SMX can synchronize and cooperate with each other via the shared memory of that SMX unit and can access the Global memory. Note that the CUDA kernels are launched by the CPU but they get executed on the GPU. Thus compute architecture of the GPU requires CPU to initiate the CUDA kernels.

CUDA enables the programming of Nvidia GPUs by exposing low level API. Apart from CUDA framework, Nvidia provides a wide variety of other tools and also supports third party libraries that can be used to program Nvidia GPUs. Since a major chunk of the scientific computing algorithms is linear algebra based, it is not surprising that the standard linear algebraic solver libraries like BLAS and Linear Algebra PACKage (LAPACK) also have their equivalents for Nvidia GPUs in one form or another. Unlike CUDA APIs, such libraries expose APIs at a much higher-level and mask the architectural details of the underlying GPU hardware to some extent thus enabling relatively faster development time.

Considering the tradeoffs between the algorithm’s computational requirements, design flexibility, execution speed and development time, we choose CULA-Dense as our main implementation library. CULA-Dense provides GPU based implementations of the LAPACK and BLAS libraries for dense linear algebra and contains routines for systems solvers, singular value decompositions, and eigen-problems. Along with the rich set of functions that it offers, CULA provides the flexibility needed by the programmer to rapidly implement the algorithm while maintaining the performance. It hides most of the GPU architecture dependent programming details thus making it possible for rapid prototyping of GPU intensive routines.

The data transfers between the CPU memory and the GPU memory are usually explicitly initiated by CPU and are carried out via the PCI (or PCI Express) bus interconnecting the CPU and the GPU. The movement of data buffers between CPU and GPU is the most taxing in terms of time. The buffer transaction time is shown in the plot in Figure 9. Newer GPUs, like Kepler based GPUs, also support useful features like GPU-GPU direct data transfers without CPU intervention. Our system and software specifications are given in Table 3.

CULA exposes two important interfaces for GPU programming namely, standard and device. Using the standard interface, the developer can program without worrying about the underlying architectural details of the GPU as the standard interface takes care of all the data movements, memory allocations in the GPU and synchronization issues. This however comes at a cost. For every standard interface function call the data is moved in and out of the GPU even if the output result of one operation is directly required by the subsequent operation. This unnecessary movement of intermediate data can dramatically impact the performance of the program. In order to avoid this, CULA provides the device interface. We use the device interface for STGD in which the programmer is responsible for data buffer allocations in the GPU memory, the required data movements between the CPU and GPU, and operates only on the data in the GPU. Thus the subroutines of the program that are iterative in nature are good candidates for device implementation.

The pre-processing involves matrices whose leading dimension is of the order of number of nodes. These are implemented using the CULA standard interface BLAS II and BLAS III routines.

Pre-processing requires SVD computations for the Moore-Penrose pseudoinverse calculations. We use CULA SVD routines since these SVD operations are carried out on matrices of moderate size. We further replaced the CULA SVD routines with more scalable SVD and pseudo inverse routines using random projections (Gittens and Mahoney, 2013) to handle larger datasets such as DBLP dataset in our experiment.

After STGD, the community membership matrix estimates are obtained using BLAS III routines provided by the CULA standard interface. The matrices are then used for hypothesis testing to evaluate the algorithm against the ground truth.

D Results on Synthetic Datasets

Homophily is an important factor in social interactions (McPherson et al., 2001); the term homophily refers to the tendency that actors in the same community interact more than across different communities. Therefore, we assume diagonal dominated community connectivity matrix $P$ with diagonal elements equal to $0.9$ and off-diagonal elements equal to $0.1$. Note that $P$ need neither be stochastic nor symmetric. Our algorithm allows for randomly generated community connectivity matrix $P$ with support $$. In this way, we look at general directed social ties among communities.

We perform experiments for both the stochastic block model ($\alpha_{0}=0$) and the mixed membership model. For the mixed membership model, we set the concentration parameter $\alpha_{0}=1$. We note that the error is around $8\%-14\%$ and the running times are under a minute, when $n\leq 10000$ and $n\gg k$.

The results are given in Table 12. We observe that more samples result in a more accurate recovery of memberships which matches intuition and theory. Overall, our learning algorithm performs better in the stochastic block model case than in the mixed membership model case although we note that the accuracy is quite high for practical purposes. Theoretically, this is expected since smaller concentration parameter $\alpha_{0}$ is easier for our algorithm to learn (Anandkumar et al., 2013b). Also, our algorithm is scalable to an order of magnitude more in $n$ as illustrated by experiments on real-world large-scale datasets.

E Comparison of Error Scores

Normalized Mutual Information (NMI) score (Lancichinetti et al., 2009) is another popular score which is defined differently for overlapping and non-overlapping community models. For non-overlapping block model, ground truth membership for node $i$ is a discrete $k$-state categorical variable $\Pi_{\text{block}}\in[k]$ and the estimated membership is a discrete $\widehat{k}$-state categorical variable $\widehat{\Pi}_{\text{block}}\in[\widehat{k}]$. The empirical distribution of ground truth membership categorical variable $\Pi_{\text{block}}$ is easy to obtain. Similarly is the empirical distribution of the estimated membership categorical variable $\widehat{\Pi}_{\text{block}}$. NMI for block model is defined as

where $n_{m}$ is the number of nodes in community $m$. The same holds for $\widehat{\Pi}_{\text{mix}_{m}}$. The normalized conditional entropy between $\mathbf{\Pi}_{{\text{mix}}}$ and $\mathbf{\widehat{\Pi}}_{{\text{mix}}}$ is defined as

where $\Pi_{\text{mix}_{j}}$ denotes the $j^{th}$ entry of $\mathbf{\Pi}_{\text{mix}}$ and similarly for $\widehat{\Pi}_{\text{mix}_{i}}$. The NMI for overlapping community is

There are two aspects in evaluating the error. The first aspect is the $l_{1}$ norm error. According to Equation (26), the error function used in NMI score is $\frac{H\left(\widehat{\Pi}_{{\text{mix}}_{i}}\lvert\Pi_{{\text{mix}}_{j}}\right)}{H(\Pi_{{\text{mix}}_{j}})}$. NMI is not suitable for evaluating recovery of different sized communities. In the special case of a pair of extremely sparse and dense membership vectors, depicted in Figure 10, ${H(\Pi_{\text{mix}_{j}})}$ is the same for both the dense and the sparse vectors since they are flipped versions of each other (0s flipped to 1s and vice versa). However, the smaller sized community (i.e. the sparser community vector), shown in red in Figure 10, is significantly more difficult to recover than the larger sized community shown in blue in Figure 10. Although this example is an extreme scenario that is not seen in practice, it justifies the drawbacks of the NMI. Thus, NMI is not suitable for evaluating recovery of different sized communities.

In contrast, our error function employs a normalized $l_{1}$ norm error which penalizes more for larger sized communities than smaller ones.

The second aspect is the error induced by false pairings of estimated and ground-truth communities. NMI score selects only the closest estimated community through normalized conditional entropy minimization and it does not account for statistically significant dependence between an estimated community and multiple ground truth communities and vice-versa, and therefore it underestimates error. However, our error score does not limit to a matching between the estimated and ground truth communities: if an estimated community is found to have statistically significant correlation with multiple ground truth communities (as evaluated by the $p$-value), we penalize for the error over all such ground truth communities. Thus, our error score is a harsher measure of evaluation than NMI. This notion of “soft-matching” between ground-truth and estimated communities also enables validation of recovery of a combinatorial union of communities instead of single ones.

A number of other scores such as “separability”, “density”, “cohesiveness” and “clustering coefficient” (Yang and Leskovec, 2012) are non-statistical measures of faithful community recovery. The scores of (Yang and Leskovec, 2012) intrinsically aim to evaluate the level of clustering within a community. However our goal is to measure the accuracy of recovery of the communities and not how well-clustered the communities are.

Banerjee and Langford (Banerjee and Langford, 2004) proposed an objective evaluation criterion for clustering which use classification performance as the evaluation measure. In contrast, we look at how well the method performs in recovering the hidden communities, and we are not evaluating predictive performance. Therefore, this measure is not used in our evaluation.

Finally, we note that cophenetic correlation is another statistical score used for evaluating clustering methods, but note that it is only valid for hierarchical clustering and it is a measure of how faithfully a dendrogram preserves the pairwise distances between the original unmodeled data points (Sokal and Rohlf, 1962). Hence, it is not employed in this paper.