Pseudo-likelihood methods for community detection in large sparse networks

Introduction

Analysis of network data is important in a range of disciplines and applications, appearing in such diverse areas as sociology, epidemiology, computer science, and national security, to name a few. Network data here refers to observed edges between nodes, possibly accompanied by additional information on the nodes and/or the edges, for example, edge weights. One of the fundamental questions in analysis of such data is detecting and modeling community structure within the network. A lot of algorithmic approaches to community detection have been proposed, particularly in the physics literature; see NewmanPNAS , Fortunato2010 for reviews. These include various greedy methods such as hierarchical clustering (see Newman2004Review for a review) and algorithms based on optimizing a global criterion over all possible partitions, such as normalized cuts Shi00 and modularity Newman&Girvan2004 . The statistics literature has been more focused on model-based methods, which postulate and fit a probabilistic model for a network with communities. These include the popular stochastic block model Holland83 , its extensions to include varying degree distributions within communities Karrer10 and overlapping communities Airoldi2008 , Hall&Karrer&Newman2011 , and various latent variable models Handcock2007 , Hoff2007 .

The stochastic block model is perhaps the most commonly used and best studied model for community detection. For a network with $n$ nodes defined by its $n\times n$ adjacency matrix $A$, this model postulates that the true node labels $c=(c_{1},\ldots,c_{n})\in\{1,\ldots,K\}^{n}$ are drawn independently from the multinomial distribution with parameter $\pi=(\pi_{1},\ldots,\pi_{K})$, where $\pi_{i}>0$ for all $i$, and $K$ is the number of communities, assumed known. Conditional on the labels, the edge variables $A_{ij}$ for $i<j$ are independent Bernoulli variables with

where $P=[P_{ab}]$ is a $K\times K$ symmetric matrix. The network is undirected, so $A_{ji}=A_{ij}$, and $A_{ii}=0$ (no self-loops). The problem of community detection is then to infer the node labels $c$ from $A$, which typically also involves estimating $\pi$ and $P$.

Fitting block models is nontrivial, especially for large networks, since in principle the problem of optimizing over all possible label assignments is NP-hard. In the Bayesian framework, Markov Chain Monte Carlo methods have been developed Snijders&Nowicki1997 , Nowicki2001 , but they only work for networks with a few hundred nodes. Variational methods have also been developed and studied (see, e.g., Airoldi2008 , Celisseetal2011 , Mariadassouetal2010 , Bickel&Choi&etal2012 ), and are generally substantially faster than the Gibbs sampling involved in MCMC, but still do not scale to the order of a million nodes. Another Bayesian approach based on a belief propagation algorithm was proposed recently by Decelle et al. Decelleetal2011 , and is comparable to ours in theoretical complexity, but slower in practice; see more on this in Section 4.

In the non-Bayesian framework, a profile likelihood approach was proposed in Bickel&Chen2009 : since for a given label assignment parameters can be estimated trivially by plug-in, they can be profiled out and the resulting criterion can be maximized over all label assignments by greedy search. The same method is used in Karrer10 to fit the degree-corrected block model. The speed of the profile likelihood algorithms depends on exactly what search method is used and the number of iterations it is run for, but again these generally work well for thousands but not millions of nodes. A method of moments approach was proposed in Bickel&Chen&Levina2011 , for a large class of network models that includes the block model as a special case. The generality of this method is an advantage, but it involves counting all occurrences of specific patterns in the graph, which is computationally challenging beyond simple special cases. Some faster approximations for block model fitting based on spectral representations are also available Newman2006 , Rohe2011 , but the properties of these approximations are only partially known.

Profile likelihood methods have been proven to give consistent estimates of the labels when the degree of the graph grows with the number of nodes, under both the stochastic block models Bickel&Chen2009 and the degree-corrected version Zhaoetal2012 . To obtain “strong consistency” of the labels, that is, the probability of the estimated label vector being equal to the truth converging to 1, the average graph degree $\lambda_{n}$ has to grow faster than $\log n$, where $n$ is the number of nodes. To obtain “weak consistency,” that is, the fraction of misclassified nodes converging to 0, one only needs $\lambda_{n}\rightarrow\infty$. Asymptotic behavior of variational methods is studied in Celisseetal2011 and Bickel&Choi&etal2012 , and in Decelleetal2011 this belief propagation method is analyzed for both the sparse [$\lambda_{n}=O(1)$] and the dense ($\lambda_{n}\rightarrow\infty$) regimes, by nonrigorous cavity methods from physics, and a phase transition threshold, below which the labels cannot be recovered, is established. In fact, it is easy to see that consistency is impossible to achieve unless $\lambda_{n}\rightarrow\infty$, since otherwise the expected fraction of isolated nodes does not go to 0. The results one can get for the sparse case, such as Decelleetal2011 , can only claim that the estimated labels are correlated with the truth better than random guessing, but not that they are consistent. In this paper, for the purposes of theory we focus on consistency and thus necessarily assume that the degree grows with $n$. However, in practice we find that our methods are very well suited for sparse networks and work well on graphs with quite small degrees.

Our main contribution here is a new fast pseudo-likelihood algorithm for fitting the block model, as well as its variation conditional on node degrees that allows for fitting networks with highly variable node degrees within communities. The idea of pseudo-likelihood dates back to Besag74 , and in general amounts to ignoring some of the dependency structure of the data in order to simplify the likelihood and make it more tractable. The main feature of the adjacency matrix we ignore here is its symmetry; we also apply block compression, that is, divide the nodes into blocks and only look at the likelihood of the row sums within blocks. This leads to an accurate and fast approximation to the block model likelihood, which allows us to easily fit block models to networks with tens of millions of nodes. Another major contribution of the paper is the consistency proof of one step of the algorithm. The proof requires new and somewhat delicate arguments not previously used in consistency proofs for networks; in particular, we use the device of assuming an initial value that has a certain overlap with the truth, and then show the amount of overlap can be arbitrarily close to purely random. Finally, we propose spectral clustering with perturbations, a new clustering method of independent interest which we use to initialize pseudo-likelihood in practice. For sparse networks, regular spectral clustering often performs very poorly, likely due to the presence of many disconnected components. We perturb the network by adding additional weak edges to connect these components, resulting in regularized spectral clustering which performs well under a wide range of settings.

The rest of the paper is organized as follows. We present the algorithms in Section 2, and prove asymptotic consistency of pseudo-likelihood in Section 3. The numerical performance of the methods is demonstrated on a range of simulated networks in Section 4 and on a network of political blogs in Section 5. Section 7 concludes with discussion, and the Appendix contains some additional technical results.

Algorithms

The joint likelihood of $A$ and $c$ could in principle be maximized via the expectation–maximization (EM) algorithm, but the E-step involves optimizing over all possible label assignments, which is NP-hard. Instead, we introduce an initial labeling vector $e=(e_{1},\ldots,e_{n})$, $e_{i}\in\{1,\ldots,K\}$, which partitions the nodes into $K$ groups. Note that for convenience we partition into the same number of groups as we assume to exist in the true model, but in principle the same idea can be applied with a different number of groups; in fact dividing the nodes into $n$ groups with a single node in each group instead gives an algorithm equivalent to that of Newman&Leicht2007 .

The main quantity we work with are the block sums along the columns,

for $i=1,\ldots,n$, $k=1,\ldots,K$. Let $\mathbf{b}_{i}=(b_{i1},\ldots,b_{iK})$. Further, let $R$ be the $K\times K$ matrix with entries $\{R_{ka}\}$ given by

Let $R_{k\bolds\cdot}$ be the $k$th row of $R$, and let $P_{\bolds\cdot l}$ be the $l$th column of $P$. Let $\lambda_{lk}=nR_{k\bolds\cdot}P_{\bolds\cdot l}$ and $\Lambda=\{\lambda_{lk}\}$.

Our approach is based on the following key observations: for each node $i$, conditional on labels $c=(c_{1},\ldots,c_{n})$ with $c_{i}=l$: {longlist}[(B)]

$\{b_{i1},\ldots,b_{iK}\}$ are mutually independent;

$b_{ik}$, a sum of independent Bernoulli variables, is approximately Poisson with mean $\lambda_{lk}$. With true labels $\{c_{i}\}$ unknown, each ${\mathbf{b}}_{i}$ can be viewed as a mixture of Poisson vectors, identifiable as long as $\Lambda$ has no identical rows.

By ignoring the dependence among $\{\mathbf{b}_{i},i=1,\ldots,n\}$, using the Poisson assumption, treating $\{c_{i}\}$ as latent variables, and setting $\lambda_{l}=\sum_{k}\lambda_{lk}$, we can write the pseudo log-likelihood as follows (up to a constant):

For any labeling $e$, let $n_{k}(e)=\sum_{i}1(e_{i}=k)$, $n_{kl}(e)=n_{k}(e)n_{l}(e)$ if $k\neq l$, $n_{kk}(e)=n_{k}(e)(n_{k}(e)-1)$ and $O_{kl}(e)=\sum_{i,j}A_{ij}1(e_{i}=k,e_{j}=l)$. We suppress the dependence on $e$ whenever there is no ambiguity. The details of the algorithmic steps can be summarized as follows.

The pseudo-likelihood algorithm. Initialize labels $e$, and let $\hat{\pi}_{l}=n_{l}/n$, $\hat{R}=\operatorname{diag}(\hat{\pi}_{1},\ldots,\hat{\pi}_{K})$, $\hat{P}_{lk}=O_{lk}/n_{lk}$, $\hat{\lambda}_{lk}=n\hat{R}_{k\bolds\cdot}\hat{P}_{\bolds\cdot l}$, $\hat{P}=\{\hat{P}_{lk}\}$ and $\hat{\Lambda}=\{\hat{\lambda}_{lk}\}$. Then repeat $T$ times: {longlist}[(6)]

Compute the block sums $\{b_{il}\}$ according to (2).

Using current parameter estimates $\hat{\pi}$ and $\hat{\Lambda}$, estimate probabilities for node labels by

Given label probabilities, update parameter values as follows:

Return to step 2 unless the parameter estimates have converged.

Update labels by $e_{i}=\arg\max_{l}\hat{\pi}_{il}$ and return to step 1.

Update $\hat{P}$ as follows: $\hat{P}_{lk}=(\sum_{i,j}A_{ij}\hat{\pi}_{il}\hat{\pi}_{jk})/n_{lk}(e)$. In practice, in step 6 we only include the terms corresponding to $\hat{\pi}_{il}$ greater than some small threshold. The EM method fits a valid mixture model as long as the identifiability condition holds, and is thus guaranteed to converge to a stationary point of the objective function Wu1983 . Another option is to update labels after every parameter update (i.e., skip step 4). We have found empirically that the algorithm above is more stable, and converges faster. In general, we only need a few label updates until convergence, and even using $T=1$ (one-step label update) gives reasonable results with a good initial value. The choice of the initial value of $e$, on the other hand, can be important; see more on this in Section 2.3.

2 Pseudo-likelihood conditional on node degrees

For networks with hub nodes or those with substantial degree variability within communities, the block model can provide a poor fit, essentially dividing the nodes into low-degree and high-degree groups. This has been both observed empirically Karrer10 and supported by theory Zhaoetal2012 . The extension of the block model designed to cope with this situation, the degree-corrected block model Karrer10 , has an extra degree parameter to be estimated for every node, and writing out a pseudo-likelihood that lends itself to an EM-type optimization is more complicated. However, there is a simple alternative: consider the pseudo-likelihood conditional on the observed node degrees. Whether these degrees are similar or not will not then matter, and the fitted parameters will reflect the underlying block structure rather than the similarities in degrees.

The conditional pseudo-likelihood is again based on a simple observation: {longlist}[(C)]

If random variables $X_{k}$ are independent Poisson with means $\mu_{k}$, their distribution conditional on $\sum_{k}X_{k}$ is multinomial. Applying this observation to the variables $(b_{i1},\ldots,b_{iK})$, we have that their distribution, conditional on labels $c$ with $c_{i}=l$ and the node degree $d_{i}=\sum_{k}b_{ik}$, is multinomial with parameters $(d_{i};\theta_{l1},\ldots,\theta_{lK})$, where $\theta_{lk}=\frac{\lambda_{lk}}{\lambda_{l}}$. The conditional log pseudo-likelihood (up to a constant) is then given by

and the parameters can be obtained by maximizing this function via the EM algorithm for mixture models, as before. We again repeat the EM for a fixed number of iterations, updating the initial partition vector after the EM has converged. The algorithm is then the same as that for unconditional pseudo-likelihood, with steps 2 and 3 replaced by: {longlist}[$(3)$]

Based on current estimates $\hat{\pi}$ and $\{\hat{\theta}_{lk}\}$, let

Given label probabilities, update parameter values as follows:

3 Initializing the partition vector

One of the simplest possible ways to group nodes in a network is to separate them by degree, say by one-dimensional $K$-means clustering applied to the degrees as in Channarondetal2011 . This only works for certain types of block models, identifiable from their degree distributions, and in general $K$-means does not deal well with data with many ties, which is the case with degrees. Instead, we consider two-dimensional $K$-means clustering on the pairs $(d_{i},d_{i}^{(2)})$, where $d_{i}^{(2)}$ is the number of paths of length 2 from node $i$, which can be obtained by summing the rows of $A^{2}$.

3.2 Spectral clustering with perturbations

A more sophisticated clustering scheme is based on spectral properties of the adjacency matrix $A=\{A_{ij}\}$ or its graph Laplacian. Let $D=\operatorname{diag}(d_{1},\ldots,d_{n})$ be diagonal matrix collecting node degrees. A common approach is to look at the eigenvectors of the normalized graph Laplacian $L=D^{-1/2}AD^{-1/2}$, choosing a small number, say $r=K-1$, corresponding to $r$ largest (in absolute value) eigenvalues, with the largest eigenvalue omitted; see, for example, Shi00 . These vectors provide an $r$-dimensional representation for nodes of the graph, on which we can apply $K$-means to find clusters; this is one of the versions of spectral clustering, which was analyzed in the context of the block model in Rohe2011 .

We found that this version of spectral clustering tends to do poorly at community detection when applied to sparse graphs, say, with expected degree $\lambda<5$. The $r$-dimensional representation seems to collapse to a few points, likely due to the presence of many disconnected components. We have found, however, that a simple modification performs surprisingly well, even for values of $\lambda$ close to 1. The idea is to connect all disconnected components which belong to the same community by adding artificial “weak” links. To be precise, we “regularize” the adjacency matrix $A$ by adding $\alpha/p\times\lambda/n$ multiplied by the adjacency matrix of an Erdos–Renyi graph on $n$ nodes with edge probability $p$, where $\alpha$ is a constant. We found that, empirically, $\alpha/p=0.25$ works well for the range of $n$ considered in our simulations, and that the results are essentially the same for all $p>0.1$ Thus we make the simplest and computationally cheapest choice of $p=1$, adding a constant matrix of small values, namely, $0.25(\lambda/n)1_{n}1_{n}^{T}$ where $1_{n}$ is the all-ones $n$-vector, to the original adjacency matrix. The rest of the steps, that is, forming the Laplacian, obtaining the spectral representation and applying $K$-means, are performed on this regularized version of $A$. We note that to obtain the spectral representation, one only needs to know how the matrix acts on a given vector; since $(A+0.25(\lambda/n)1_{n}1_{n}^{T})x=Ax+0.25(\lambda/n)(\sum_{i}x_{i})1_{n}$, the addition of the constant perturbation does not increase computational complexity. We will refer to this algorithm as spectral clustering with perturbations (SCP), since we perturb the network by adding new, low-weight “edges.”

Consistency results

By consistency we mean consistency of node labels (to be defined precisely below) under a block model as the size of the graph $n$ grows. For the theoretical analysis, we only consider the case of $K=2$ communities. We condition on the community labels $\{c_{i}\}$, that is, we treat them as deterministic unknown parameters. For simplicity, here we consider the case of balanced communities, each having $m=n/2$ nodes. An extension to the unbalanced case is provided in the supplementary material Amietal . The assumption of balanced communities naturally leads us to use the class prior estimates $\hat{\pi}_{1}=\hat{\pi}_{2}=1/2$ in (10). We call this assumption (E) (for equal class sizes): {longlist}[(E)]

Assume each class contains $m=n/2$ nodes, and set $\hat{\pi}_{1}=\hat{\pi}_{2}=1/2$. Without loss of generality, we can take $c_{i}=1$ for $i\in\{1,2,\ldots,m\}$.

As an intermediate step in proving consistency for the block model introduced in Section 1, we first prove the result for a directed block model. Recall that for the (undirected) block model introduced earlier, one has

In the directed case, we assume that all the entries in the adjacency matrix are drawn independently, that is,

We will use different symbols for the adjacency and edge-probability matrices in the two cases. This is to avoid confusion when we need to introduce a coupling between the two models. In both cases, we have assumed that diagonal entries of the adjacency matrices are also drawn randomly (i.e., we allow for self-loops as valid within-community edges). This is convenient in the analysis with minor effect on the results.

The directed model is a natural extension of the block model when one considers the pseudo-likelihood approach; in particular, it is the model for which the pseudo-likelihood assumption of independence holds. It is also a useful model of independent interest in many practical situations, in which there is a natural direction to the link between nodes, for example, in email, web, routing and some social networks. The model can be traced back to the work of Holland and Leinhardt HolLei81 and Wang and Wong wang1987 in which it has been implicitly studied in the context of more general exponential families of distributions for directed random graphs.

Our approach is to prove a consistency result for the directed model, with an edge-probability matrix of the form

Note that the only additional restriction we are imposing is that $\widetilde{P}$ has the same diagonal entries. Both $a$ and $b$ depend on $n$ and can in principle change with $n$ at different rates. This is a slightly different parametrization from the more conventional $P_{n}=\rho_{n}S$ Bickel&Chen2009 , where $S$ (and $\pi$) do not depend on $n$, and $\lambda_{n}=\rho_{n}\pi^{T}S\pi$. We use this particular parametrization here because we only consider the case $K=2$, and it makes our results more directly comparable to those obtained in the physics literature, for example, Decelleetal2011 .

A coupling between the directed and the undirected model that we will introduce allows us to carry the consistency result over to the undirected model, with the edge-probability matrix

Asymptotically, the two edge-probability matrices have comparable (to first order) expected degree and out-in-ratio (as defined by Decelleetal2011 ), under mild assumptions. The average degrees for $\widetilde{P}$ and $P$ are $a+b$ and $2(a+b)-\frac{1}{m}(a^{2}+b^{2})$, respectively. The latter is $\sim 2(a+b)$ as long as $\frac{1}{2m}\frac{a^{2}+b^{2}}{a+b}\leq\frac{a+b}{n}\to 0$. The condition is satisfied as soon as the average degree of the directed model has sublinear growth: $a+b=o(n)$. The same holds for out-in-ratios.

For our analysis, we consider an E-step of the CPL algorithm. It starts from some initial estimates $\hat{a}$, $\hat{b}$ and $\hat{\pi}=(\hat{\pi}_{1},\hat{\pi}_{2})$ of parameters $a$, $b$ and $\pi$, together with an initial labeling $e$, and outputs the label estimates

The key assumption of our analysis is that the initial labeling has a certain overlap with the truth (we will show later that the amount of overlap is not important). One situation where this might naturally arise is survey data, when some small fraction of nodes has been surveyed about their community membership. Another possibility is to run some other crude algorithm first to obtain a preliminary result. More formally, we consider an initial labeling $e=(e_{i})\in\{1,2\}^{n}$, which is balanced (i.e., assigns equal number of nodes to each label) and matches exactly $\gamma m$ labels in community $1$, for some $\gamma\in(0,1)$. We do not assume that we know which labels are matched, or the value of $\gamma$. It is easy to see that this is equivalent to $e$ matching exactly $\gamma m$ labels in each of the two communities. Assuming $\gamma m$ to be an integer, let $\mathcal{E}^{\gamma}=\mathcal{E}^{\gamma}_{n}$ denote the collection of such labelings,

Let us consider the directed case first. As our measure of performance (i.e., the loss function), we take the following (directed-case) mismatch ratio

where $\hat{c}_{i}(e)$ are computed based on the directed adjacency matrix $\widetilde{A}$, and $\{(1\,2),(2\,1)\}$ is the set of permutations of $\{1,2\}$, with $\phi$ accounting for the fact that the labels assigned to the communities are only determined up to a permutation. The counterpart for the undirected case is denoted by $M_{n}(e)$. Note that the notion of consistency based on convergence of this quantity matches the “weak” consistency discussed in Zhaoetal2012 , rather than the “strong” consistency used by Bickel&Chen2009 . Define

and let $h(p)=-p\log p-(1-p)\log(1-p)$, $p\in$ be the binary entropy function. Let us also consider the collection of estimates $(\hat{a},\hat{b})$ which have the same ordering as true parameters $(a,b)$,

where $\kappa_{\gamma}(n):=\frac{1}{n}[\log(\frac{n}{4\pi\gamma(1-\gamma)})+\frac{1}{3n}]=o(1)$.

In particular, if $\tau_{n}^{2}\to\infty$, we have $u_{n}\to\infty$ and the CPL estimate is uniformly consistent.

We think of $\gamma$ as fixed, but it is possible to let $\gamma=\gamma_{n}\to\frac{1}{2}$, making the problem harder as $n$ grows. We still get consistency as long as $(1-2\gamma_{n})^{2}\tau_{n}^{2}\to\infty$.

In the balanced case, the CPL iteration has a simple intuitive interpretation, as will become clear during the proof of Theorem 1. One starts with an initial assignment of labels to nodes. Then, each node updates its label by taking a majority vote among its neighbors. In the case where $b=0$, it is intuitively clear that for $a$ large enough, this procedure

increases the number of correct labels relative to the initial assignment. Figure 1 illustrates these ideas. In the general case where $b\neq 0$, Theorem 1 states that $\tau_{n}^{2}$ is the key parameter that needs to grow for the procedure to succeed.

While the labels are of primary interest in community detection, one may also be interested in consistency of the estimated parameters. Under strong consistency in the sense of Bickel&Chen2009 , consistency of the natural plug-in estimates of the block model parameters follows easily, but here we only show weak consistency of the labels. However, in the directed model the pseudo-likelihood function we defined is in fact exactly the likelihood of $\mathbf{b}_{i}$’s. Parameter estimates (say $\hat{a}$ and $\hat{b}$) obtained by the EM algorithm converge to a local maximum of this function. As a consequence of Theorem 1, these estimates are also consistent (for $a$ and $b$). Since the likelihood is smooth with bounded derivatives, one may be able to use standard arguments to show that the estimated parameters are a unique local maximum in a neighborhood of the truth, and even derive their asymptotic normality along; see, for example, Theorem 6.2.1, page 384 of Bickel&Doksum . We do not pursue this direction here.

Assume (E), and let $\gamma\in(0,1)\setminus\{\frac{1}{2}\}$. Let the adjacency matrix $A$ be generated according to the undirected model (6) with edge-probability matrix (9), and assume $a\neq b$. In addition, assume

where $\kappa_{\gamma}(n)=o(1)$ is as defined in Theorem 1.

In particular, if $\tau_{n}^{2},a_{\gamma}\to\infty$, we have $u_{n},v_{n}\to\infty$, and the CPL estimate is uniformly consistent.

The proofs of both theorems can be found in Section 6.

Condition (17) can be met for a fixed $\varepsilon\in(0,1)$ by choosing $\gamma$ sufficiently small and an upper bound on $b/a$ in terms of $\gamma$. For example, for $\varepsilon=\frac{1}{2}$ and $\gamma<\frac{1}{8}$, we have (17) if

The parameter $\tau_{n}^{2}$ controlling consistency is the same as the one reported in Decelleetal2011 and Mosseletal2012 . There the concern is with recovering a labeling which is positively correlated with the truth, and the threshold of success is observed to be $\tau_{n}^{2}\geq 2$. A similar lower bound was given in Chaudhuri&Chung&Tsiatas2012 for spectral clustering. Here, we are concerned with moving from a positively correlated labeling to one with an asymptotically vanishing mismatch ratio [i.e., $\widetilde{M}_{n}(e)=o_{p}(1)$], which is why we need $\tau_{n}^{2}\to\infty$.

These results can be extended to the case of unbalanced communities. Such an extension is provided for the directed block model in the supplementary material Amietal . There we consider the model with two communities of sizes $n_{1}$ and $n_{2}$ (not necessarily equal) and an edge-probability matrix

which relaxes our earlier assumption $a_{1}=a_{2}$ in (8). The class of initial labelings is also enlarged to include those that have $\gamma_{k}$-overlap with community $k$, that is, $\mathcal{E}^{\gamma_{1},\gamma_{2}}:=\{e\dvtx\sum_{i}1_{\{e_{i}=k,c_{i}=k\}}=\gamma_{k}n_{k},k=1,2\}$, with $\gamma_{1}\neq\gamma_{2}$. In this situation, one needs more assumptions on the initial estimate $\hat{P}$ used in the CPL iteration than in the balanced case. Supplementary material Amietal gives the details. While we do not discuss the undirected case in this general setting, ideas used in the proof of Theorem 2 can be used to carry the results from the directed to the undirected case.

Numerical results

The matrix $P$ is constructed as follows. It is controlled by two parameters: the “out-in-ratio” $\beta$ Decelleetal2011 , which we will vary from 0 to 0.2, and the weight vector $w$, which determines the relative degrees within communities. We consider two values of $w$: $w=(1,1,1)$ (no information about communities is contained in node degrees) and $w=(1,5,10)$ (degrees themselves provide relevant information for clustering). If $\beta=0$, we set $P^{(0)}=\operatorname{diag}(w)$, a diagonal matrix. Otherwise, we set the diagonal of $P^{(0)}$ to $\beta^{-1}w$ and set all off-diagonal elements to $1$. We then fix the overall expected network degree $\lambda$, which is the natural parameter to control Bickel&Chen2009 and which we will vary from 1 to 15. Then we rescale $P^{(0)}$ to obtain this expected degree, giving the final $P$

To compare our results to the true labels, we will use normalized mutual information (NMI). One can think of the confusion matrix $R$ as a bivariate probability distribution, and of its row and column sums $R_{i+}$ and $R_{+j}$ as the corresponding marginals. Then the NMI is defined by Yao03 as $\operatorname{NMI}(c,e)=-\sum_{i,j}R_{ij}\log\frac{R_{ij}}{R_{i+}R_{+j}}(\sum_{i,j}R_{ij}\log R_{ij})^{-1}$, and is always a number between 0 and 1 (perfect match). It is useful to have a few benchmark values of NMI for reference: for example, for large $n$, matching $50\%$, $70\%$ and $90\%$ of the labels correspond to values of NMI of approximately $0.12$, $0.26$ and $0.58$, respectively.

All figures show the performance of the following methods: $K$-means clustering on 1- and 2-degrees (DC), spectral clustering (SC), spectral clustering with perturbations (SCP), unconditional pseudo-likelihood (UPL) initialized with either DC or SCP, and conditional pseudo-likelihood (CPL), with the same two initial values for labelings. The number of outer iterations for UPL and CPL is set to $T=20$; $n$, $\lambda$, $\rho$ and the number of replications $N$ are specified in the figures.

Figures 2 and 3 show results on estimating the node labels with varying $\beta$ and $\lambda$, respectively. Generally, smaller $\beta$ and larger $\lambda$ make the problem easier, as we expect. In principle, degree-based clustering gives no information about the labels with uniform weights $w$, and only a moderate amount of information with nonuniform weights, so it serves as an example of a poor starting value for pseudo-likelihood. Regular spectral clustering performs well with uniform weights, but very poorly with nonuniform weights; we conjecture that this is due to a limitation of $K$-means. Spectral clustering with perturbation, on the other hand, performs very well in all scenarios. Apart from being a useful general method on its own, it also serves as an example of a good starting value for pseudo-likelihood.

Figures 2 and 3 show that pseudo-likelihood achieves large gains over a poor starting value, giving surprisingly good results even when starting from the uninformative degree clustering in the case of $w=(1,1,1)$. One exception is unconditional pseudo-likelihood with $\rho=0.9$ and $w=(1,1,1)$, which shows that conditioning is necessary to accommodate variation in degrees when the starting value is not very good. When spectral clustering with perturbation is used as a starting value, which is already very good, UPL and CPL do not have much room to do better, although UPL still provides a noticeable improvement, being overall the best method when initialized with SCP. It appears that a good starting value overcomes the limitations of the regular block model for networks with hubs, effectively ruling out the competing solution which divides nodes by degree.

Finally, Figure 4 shows run times for all the methods for the case of the regular block model ($\rho=0$) with different community weights [$w=(1,1,1)$ and $w=(1,5,10)$]. The times shown for UPL and CPL do not include the time to compute the initial value, which is shown separately. For the case $w=(1,1,1)$, all methods take roughly the same amount of time. For the case $w=(1,5,10)$, spectral clustering (SC) takes considerably more time than the rest. On the other hand, SCP takes nearly the same time as it takes for $w=(1,1,1)$, and it slightly outperforms DC for larger values of $n$. This might be explained, in part, by the sparse matrix multiplication required for DC, which is both time and memory-consuming for large $n$. Generally, SCP provides an excellent starting value, with low computational complexity in a variety of situations.

We have also done some brief comparisons with the belief propagation (BP) method of Decelleetal2011 . Direct fair comparison is difficult because of the different platform for the belief propagation code and the different way in which it handles initial values; generally, we found that while the computing time of belief propagation scales with $n$ at the same rate as ours, BP is slower by a constant factor of about 10. In terms of accuracy of community detection, in the examples we tried BP was either similar to or a little worse than pseudo-likelihood.

Example: A political blogs network

This dataset on political blogs was compiled by Adamic and Glance Adamic05 soon after the 2004 U.S. presidential election. The nodes are blogs focused on US politics, and the edges are hyperlinks between these blogs. Each blog was manually labeled as liberal or conservative in Adamic05 , and we treat these as true community labels. Following Karrer10 , we ignore directions of the hyperlinks and analyze the largest connected component of this network, which has 1222 nodes and the average degree of 27. The distribution of degrees is highly skewed to the right (the median degree is 13, and the maximum is 351).

The results in Figure 5 show that the conditional pseudo-likelihood produces a result closest to the truth, as one would expect in view of highly variable degrees. Its result is also very close to those obtained by profile maximum likelihood for the degree-corrected block model and by two different modularities Karrer10 , Zhaoetal2012 . Unconditional pseudo-likelihood, on the other hand, puts high-degree nodes in one group and low-degree nodes in the other. This is very close to the block model solution Karrer10 . This example confirms that the unconditional and conditional pseudo-likelihood methods are correctly fitting the block model and the degree-corrected block model, respectively.

Proofs of consistency results

Due to symmetry, we can assume without loss of generality that $\gamma\in(0,\frac{1}{2})$. Similarly, we can assume $a>b$. Then, for any $(\hat{a},\hat{b})\in\mathcal{P}_{a,b}$ we have $\hat{a}>\hat{b}$. These will be our standing assumptions throughout the proofs. To see that the assumptions are not restrictive, one can check that the proof goes through, without change, if $\gamma\in(\frac{1}{2},1)$ and $b>a$. For the other two cases, namely, $\gamma\in(0,\frac{1}{2})$ and $b>a$, or $\gamma\in(\frac{1}{2},1)$ and $a>b$, the proof goes through by switching the estimated labels when matching them with the true labels. That is, we compare estimated community $1$ to true community $2$ and vice versa. These can seen by examining (21) and the discussion that follows.

Under the equal priors assumption (E), the CPL estimate (10) simplifies to

where $\{\widetilde{b}_{im}\}$ are obtained by block compression of the directed adjacency matrix $\widetilde{A}$.

Let us focus on $i\in\mathcal{C}_{1}$ from now on. Then $\hat{c}_{i}(e)=1$ if

where $R(e)$ is defined in (3). It is then not hard to see that after row normalization of $\hat{\Lambda}=[nR(e)\hat{P}]^{T}$, we obtain $\hat{\theta}_{11}(e)=\hat{\theta}_{22}(e)=\gamma\frac{\hat{a}}{\hat{a}+\hat{b}}+(1-\gamma)\frac{\hat{b}}{\hat{a}+\hat{b}}$, and $\hat{\theta}_{12}(e)=\hat{\theta}_{21}(e)=\gamma\frac{\hat{b}}{\hat{a}+\hat{b}}+(1-\gamma)\frac{\hat{a}}{\hat{a}+\hat{b}}$.

Since by assumption $\hat{a}>\hat{b}$ and $\gamma\in(0,\frac{1}{2})$, it follows that $\hat{\theta}_{11}<\hat{\theta}_{21}$. Then, (21) is equivalent to $\widetilde{b}_{i1}(e)-\widetilde{b}_{i2}(e)<0$. Recalling that $\widetilde{b}_{ik}(e)=\sum_{j=1}^{m}\widetilde{A}_{ij}1\{e_{i}=k\}=\sum_{j\in\mathcal{S}_{k}}\widetilde{A}_{ij}$, we can write the condition as

and $\sigma(e)=(\sigma_{1}(e),\ldots,\sigma_{n}(e))$. Let $\Sigma^{\gamma}=\Sigma^{\gamma}_{n}$ be the set of all $\sigma(e)$ with $e\in\mathcal{E}^{\gamma}$, that is,

Since we are focusing on $i\in\mathcal{C}_{1}$, we are concerned with $\widetilde{M}_{n,1}(e)$. In a slight abuse of notation, $\widetilde{M}_{n}(e)$ in (22) is in fact an upper bound on the mismatch ratio as defined in (12), since here we are using a particular permutation—the identity.

Let us define, for $\sigma\in\{-1,+1\}^{n}$ and $r\geq 0$,

where the inequality is due to treating the ambiguous case $\widetilde{\xi}_{i}(\sigma)=0$ as error. We now set out to bound this in probability. Let us start with a tail bound on $\widetilde{\xi}_{i}(\sigma)$ for fixed $\sigma$ and $i$.

For any $\sigma\in\Sigma^{\gamma}$ and $t\in(0,3(a+b)]$, we have

Noting that for $t/3\leq(a+b)$, we have $2(v+t/3)\leq 4(a+b)$ completes the proof.

We also need a tail bound on $\widetilde{N}_{n,1}(\sigma;r)$. Let us define

Note that these probabilities do not depend on the particular value of $\sigma\in\Sigma^{\gamma}$, due to symmetry. We have the following lemma.

Follows from Lemma 5 in the Appendix, by noting that$\{1\{\widetilde{\xi}_{i}(\sigma)\geq-r\}\}_{i=1}^{m}$ are independent Bernoulli random variables.

Now we apply Lemma 1 with $t=(1-2\gamma)(a-b)\leq 3(a+b)$. Note that $\frac{a-b}{a+b}\leq 1\leq\frac{3}{1-2\gamma}$, for $\gamma\in(0,\frac{1}{2})$. Noting that the RHS of (23) does not depend on $i$, and using (24), we get

The cardinality of the set $\Sigma^{\gamma}$ is ${m\choose\gamma m}^{2}\leq(e^{m[h(\gamma)+\kappa_{\gamma}(2m)]})^{2}$ where $h(\cdot)$ is the binary entropy function, and $\kappa_{\gamma}(2m)=\kappa_{\gamma}(n)$ is as defined in the statement of the theorem. (See Lemma 6 in the supplementary material Amietal for a proof.) Applying Lemma 2 with $u=u_{n}$ and the union bound, we obtain

By symmetry the same bound holds for $\sup_{\sigma}\frac{1}{m}\widetilde{N}_{n,2}(\sigma;0)$. It follows from (22) that the same holds for $\sup_{e}M_{n}(e)$. This completes the proof of Theorem 1.

2 Proof of Theorem 2 (undirected case)

Our approach is to introduce a deterministic coupling between $A$ and $\widetilde{A}$, which allows us to carry over the results of the directed case. Let

In other words, the graph of $A$ is obtained from that of $\widetilde{A}$ by removing directions. Note that

which matches the relation between (8) and (9). From (26), we also note that

Let us now upper-bound $\xi_{i}(\sigma)$ in terms of $\widetilde{\xi}_{i}(\sigma)$. Based on (27), only those $\sigma_{j}$ that are equal to $1$ contribute to the upper bound. More precisely, let $D_{ij}=A_{ij}-\widetilde{A}_{ij}\geq 0$, and take $i\in\mathcal{C}_{1}$ from now on. Then

We further notice that $D_{ij}\leq\widetilde{A}_{ij}+\widetilde{A}_{ji}$. To simplify notation, let us define

where the dependence on $\sigma$ is due to $\mathcal{S}_{1}$ being derived from $\sigma$ [recall that $\mathcal{S}_{1}=\mathcal{S}_{1}(\sigma)=\{j\dvtx\sigma_{j}=1\}$]. Thus we have shown

Recall from definition (16) that $a_{\gamma}=\gamma a+(1-\gamma)b.$

Fix $\varepsilon>0$. For $i\in\mathcal{C}_{1}$, we have

The equality of the two probabilities follows by symmetry. Let us prove the bound for $\widetilde{A}_{i*}(\sigma)$. We apply Bernstein inequality. Note that

Since $\sum_{j\in\mathcal{S}_{1}}\operatorname{var}(\widetilde{A}_{ij})\leq\mu$, we obtain

Setting $t=\varepsilon\mu$ completes the proof.

which $\vee$ is the logical OR. This can be seen (as usual) by noting that if the RHS does not hold, then $\widetilde{\xi}_{i}(\sigma)+\widetilde{A}_{i*}(\sigma)+\widetilde{A}_{*i}(\sigma)<0$, implying $\xi_{i}(\sigma)<0$. Translating to indicator functions,

Averaging over $i\in\mathcal{C}_{1}$ (i.e., applying $m^{-1}\sum_{i=1}^{m}$), we get

where $\widetilde{Q}_{n,1*}(\sigma;t)=\sum_{i=1}^{m}1\{\widetilde{A}_{i*}(\sigma)\geq t\}$, and similarly for $\widetilde{Q}_{n,*1}(\sigma;t)$. Note that $\widetilde{Q}_{n,1*}(\sigma;t)$ and $\widetilde{Q}_{n,*1}(\sigma;t)$, while not independent, have the same distribution by symmetry, so we can focus on bounding one of them. The key is that each one is a sum of i.i.d. terms, for example, $\{\widetilde{A}_{i*}\}_{i=1}^{m}$.

We have a bound on $m^{-1}\widetilde{N}_{n,1}(\sigma;r)$ from Lemma 2. We can get similar bounds on the $\widetilde{Q}$-terms. To start, let

similar to (24), and note that these quantities too are independent of the particular choice of $\sigma\in\Sigma^{\gamma}$.

Follows from Lemma 5 in the Appendix, by noting that$\{1\{\widetilde{A}_{i*}(\sigma)\geq r/2\}\}_{i=1}^{m}$ is an independent sequence of Bernoulli variables.

The same bound holds for $\frac{1}{m}\widetilde{Q}_{n,*1}(\sigma;r/2)$. Recall the definition of $\bar{p}_{1}(r)$ from (24). Using (31) and Lemmas 2 and 4, we get

as long as $u_{n},v_{n}>1/e$. Now, take $r/2=(1+\varepsilon)a_{\gamma}$, so that Lemma 3 implies

Now, in Lemma 1, take $t=(1-2\gamma)(a-b)-2(1+\varepsilon)a_{\gamma}$. Note that the assumption

implies $t\geq(1-\varepsilon)(1-2\gamma)(a-b)>0$. In addition $t\leq(1-2\gamma)(a-b)\leq 3(a+b)$ as before. Thus, the chosen $t$ is valid for Lemma 1. Furthermore, $-(1-2\gamma)(a-b)+t=-r$. Hence, the lemma implies

The rest of the argument follows as in the directed case. This completes the proof of Theorem 2.

Discussion

The proposed pseudo-likelihood algorithms provide fast and accurate community detection for a range of settings, including large and sparse networks, contributing to the long history of empirical success of pseudo-likelihood approximations in statistics. For the theoretical analysis, we did not focus on the convergence properties of the algorithms, since standard EM theory guarantees convergence to a local maximum as long as the underlying Poisson or multinomial mixture is identifiable. The consistency of a single iteration of the algorithm was established for an initial value that is better than purely arbitrary, as long as, roughly speaking, the graph degree grows, and there are two balanced communities with equal expected degrees. The theory shows that this local maximum is consistent, and unique in a neighborhood of the truth, so in fact there is no need to assume that EM has converged to the global maximum, an assumption which is usually made in analyzing EM-based estimates. The theoretical analysis can be extended to the general two-community model with possibly unbalanced communities, as detailed in the supplementary material Amietal . Extending our argument to more than two communities also seems possible, but that would require extremely meticulous tracking of a large number of terms which we did not pursue.

We conjecture that additional results may be obtained under weaker assumptions if one focuses simply on estimating the parameters of the block model rather than consistency of the labels, just like one can obtain results for a labeling correlated with the truth (instead of consistent) under weaker assumptions discussed in Remark 5. For example, in a very recent paper Chatterjee2012 , results are obtained under very weak assumptions for the mean squared error of estimating the block model parameter matrix $P$ (which in itself does not guarantee consistency of the labels). While the primary interest in community detection is estimating the labels rather than the parameters, we plan to investigate this further to see if and how our conditions can be relaxed.

While in theory any “reasonable” initial value guarantees convergence, in practice the choice of initial value is still important, and we have investigated a number of options empirically. Spectral clustering with perturbations, which we introduced primarily as a method to initialize pseudo-likelihood, deserves more study, both empirically (e.g., investigating the optimal choice of the tuning parameter), and theoretically. This is also a topic for future work.

Appendix: Poisson-type tail bound

Here is a lemma which we used quite often in proving consistency results in Section 6.

where we have used $(1+x)^{m}\leq\exp(mx)$. The RHS is the Chernoff bound for a Poisson random variable with mean $\mu=\sum_{i}p_{i}$, and can be optimized to yield

Acknowledgment

We would like to thank Roman Vershynin (Mathematics, University of Michigan) for highly illuminating discussions.

Extension to unbalanced communities \slink[doi]10.1214/13-AOS1138SUPP \sdatatype.pdf \sfilenameaos1138_supp.pdf \sdescriptionThis supplement contains an extension of Theorem 1 to the case of unbalanced communities.

References