Merge-split Markov chain Monte Carlo for community detection

03/16/2020 ∙ by Tiago P. Peixoto, et al. ∙ Central European University 0

We present a Markov chain Monte Carlo scheme based on merges and splits of groups that is capable of efficiently sampling from the posterior distribution of network partitions, defined according to the stochastic block model (SBM). We demonstrate how schemes based on the move of single nodes between groups systematically fail at correctly sampling from the posterior distribution even on small networks, and how our merge-split approach behaves significantly better, and improves the mixing time of the chain by several orders of magnitudes in typical cases. We also show how the scheme can be straightforwardly extended to nested versions of the SBM, yielding asymptotically exact samples of hierarchical network partitions.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

I Introduction

Community detection Fortunato (2010) is an important tool for the analysis of network data, as it allows the breaking down of the network structure into basic building blocks, thus providing a summary of the large-scale structure of a network. Among the many methods proposed for this task, those based on the statistical inference of generative models offer the most principled and robust choice Peixoto (2019a). This is due in large part to their ability to take into account and convey the statistical evidence available in the data, allowing us to prevent both overfitting — occurring when purely random fluctuations are mistaken by structure — as well as underfitting — occurring when statistically significant structure is mistaken by random fluctuations. These problems manifest themselves in non-statistical methods in their tendency to identify spurious communities in fully random networks Guimerà et al. (2004) as well as those that are non-modular but sparse Bagrow (2012), and in the failure at recognizing the existence of clear but relatively small communities in large networks Fortunato and Barthélemy (2007)

. The manner in which statistical inference allows us to overcome these obstacles is by allowing us to express not only our belief about how the data are generated, but also about how much information is needed to specify the model parameters. The combination of these two elements defines a posterior distribution of all possible community structures that are compatible with the data, ranked according to their plausibility. The posterior distribution inherently ascribes low probabilities to overly-complicated models, and yields the most parsimonious but sufficiently rich explanation for the data. Nevertheless, in case the model being used cannot perfectly encapsulate the structure in the data, the posterior distribution may still give comparatively high probabilities to multiple diverging explanations for the data. Choosing only one fit of the model, for example by selecting a single modular division with the strictly largest posterior probability, is a form of bias in this situation. Instead, it is preferable to consider the entire set of plausible partitions, by sampling from the posterior distribution, rather than maximizing it.

Although it is relatively easy to write down a mathematical expression for the posterior distribution of network partitions, at least up to a normalization constant, using a simple model family called the stochastic block model (SBM) Holland et al. (1983); Karrer and Newman (2011), actually sampling from this distribution is a different matter altogether. In fact, fitting the SBM in this way subsumes solving certain instances of combinatorical problems such as graph coloring that are known to be NP-hard, and indeed there is strong statistical evidence pointing to the existence of parameter regimes where, although it may be possible to characterize the posterior distribution of the SBM, this can only be done in non-polynomial time Decelle et al. (2011). Because of this inherent intractability, it is unlikely that there will ever be an algorithm that samples directly from the posterior distribution in an efficient manner, at least not for networks that have more than a handful of nodes, as the number of possible partitions grows super-exponentially with network size. Nevertheless, we can still approach the problem in practice using approximate methods. One attractive choice is Markov chain Monte Carlo (MCMC), which consists in moving from one partition of the network to another, with a probability conditioned only on the current state, in a manner that guarantees exact asymptotic convergence to the desired posterior distribution, provided one is able to wait long enough for the chain to “mix.” In this way, one should be able to sample effectively from the posterior in feasible instances of the problem, which would correspond to a short mixing time.

In the case of the SBM, MCMC sampling has been implemented in a variety of ways Guimerà and Sales-Pardo (2009); Peixoto (2014a); Newman and Reinert (2016); Riolo et al. (2017), but all approaches thus far employed involve moves from one partition to another that differ only in the group membership of a single node. As we will show, these “single-node” approaches fail systematically at efficiently exploring the posterior distribution of partitions, as they cannot overcome low probability “barriers” that prevent the chain from mixing. This is particularly relevant in situations where the number of groups is unknown a priori, and the movement of a single node at a time forces the Markov chain to move through low-probability states composed of groups containing only one node, if the chain is to visit partitions with different numbers of groups. These stopgap transitions prevent the MCMC from mixing in a reasonable amount of time, and result in abysmal performance even in easy instances of the problem.

In this work we present an alternative MCMC scheme that involves the simultaneous movement of several nodes at each step, via the merging, splitting and re-arrangements of groups. As we show, these moves can easily overcome the low-probability barriers that block the single-move approaches, and cause the chain to mix several orders of magnitude faster. This improvement makes MCMC usable for significantly larger networks, and allows us to explore community structure in empirical networks with more confidence and increased amount of detail.

This paper is defined as follows. In Sec. II we briefly review the inference approach to community detection. In Sec. III.1 we describe the MCMC approach with single-node moves and in Sec. III.2 we describe the version with merges and splits. In Sec. IV we investigate the performance of both algorithms in empirical networks, and in Sec. V we extend the algorithm to hierarchical models. We end in Sec. VI with a conclusion.

Ii Bayesian inference of community structure

The inference approach to community detection involves first the stipulation of a generative model for the network structure given a network partition , where is the membership of node in one of groups. For example, when using the degree-corrected stochastic block model (DC-SBM) Karrer and Newman (2011), it is assumed that a network with adjacency matrix is generated with probability


with the additional parameters and specifying the affinity between groups and expected node degrees, respectively. The inference procedure consists in sampling from the posterior distribution



is a prior probability of node partitions, and


is the marginal likelihood integrated over the remaining model parameters, weighted according to their own prior probabilities. If we make a non-informative choice for them,


we can compute the integral exactly as Peixoto (2017)


where , and . The last remaining quantity


is called the evidence, and servers as a normalization constant for the posterior distribution. Although it has an important role in the context of model selection, its value cannot be computed in closed form in the general case. Luckily this is not needed for MCMC, as we will shortly see.

The scheme above can be modified in a variety of ways, many of which can significantly improve the quality of the results. For example, it is possible to replace the noninformative priors by nested sequences of priors and hyperpriors 

Peixoto (2017) that enhance our capacity to identify small groups in large networks, more adequately describe broad degree sequences, and uncover hierarchical modular structures Peixoto (2014b). We can also incorporate the existence of edge covariates Peixoto (2018a), and the existence of latent edges for network reconstruction Peixoto (2018b, 2019b, 2020). The method we will describe in the following is suitable for all of these scenarios, and does not depend on the details of the model specification. It is able to asymptotically sample from an arbitrary desired posterior distribution, for which we use a shortcut notation


defined up to an unknown normalization constant.

Iii Markov chain Monte Carlo (MCMC)

The MCMC protocol consists of starting from an arbitrary partition of the network, and then sample a new partition from a proposal distribution


conditioned on the current partition. The new proposal is accepted with a Metropolis-Hastings (MH) Metropolis et al. (1953); Hastings (1970) probability


otherwise the move is rejected, and the next state of the chain is the same as the previous one. Note that in order to compute the ratio we need to know only up to a normalization constant, which cancels out. The use of this choice guarantees that the final transition probabilities , after considering the acceptance or rejection, fulfills the detailed balance condition


When this condition is combined with ergodicity, i.e. the proposals allow the eventual exploration of every possible partition, and lack of periodicity, i.e. the chain can return to a previous state in a number of steps that is not necessarily confined to a multiple of an integer larger than one, then, after a sufficiently large number of steps the partitions are visited by the chain with a probability that converges to the desired distribution .

What determines whether the MCMC method works in practice is the quality of the proposals . Clearly, if the proposal happens to be identical to the target distribution, the chain mixes in a single step. But if it were possible to propose partitions directly from the target distribution, we would not need MCMC in the first place. In reality we need to rely on proposals that are far away from the target distribution, but nevertheless allow the chain to converge to it in a feasible amount of time. In the following, we describe the conventional approach of proposing the change of a single node at a time, and then we introduce moves involving more than one node.

iii.1 Single-node moves

The simplest kind of move proposal that can be done is the one that involves the change of a single node from its current group to another group , according to a local move proposal , leading to


with being the probability of choosing node to perform such a move. Note that we do not need to compute this whole expression when calculating the MH acceptance criterion, and the forward and reverse proposal must involve the same node, which also means that the probability will cancel out, and hence we are allowed to choose nodes with arbitrary frequencies, as long as they are nonzero to guarantee ergodicity.

A simple choice for the local move proposal would be to select the target group uniformly at random from all possibilities, including a new previously unoccupied group, where is the number of occupied groups in partition , i.e.


Although this is straightforward to implement, it will be inefficient as soon as the typical number of groups starts to moderately increase. Suppose, for example, that the typical number of groups is around but any given node can belong with non-negligible probability to at most groups. In such a situation of all move proposals will be rejected, leading to a substantial waste of effort, which becomes worse as the number of groups increases. We can be more efficient by proposing smarter moves that are less likely to be rejected. Following Ref. Peixoto (2014a), we choose to move to existing groups with a probability given by


where is the fraction of neighbors of node that belong to group . The above means that we inspect the local neighborhood of a node , by sampling one of its neighbors at random. We then consider its group membership and sample our new group based on the frequency with which it connects to nodes of type , which is proportional to total number of edges between these two groups plus a constant, i.e. . The constant must be nonzero to guarantee ergodicity, such that every possible move can be eventually chosen (a reasonable choice is simply ). This kind of move tends to concentrate more strongly on more plausible moves, thus increasing the acceptance ratio, although the extent of this improvement will depend on the network, and how far the chain has progressed from its initial state. However, this move proposal still does not completely guarantee ergodicity because it prevents the occupation of a new group (which in turn forbids the vacancy of any group, as the MH acceptance probability becomes zero in this case). We can finally incorporate this kind of move by augmenting our proposals as


where is the probability with which the population of a new group is attempted.

These smarter move proposals can be performed efficiently by incorporating only a small amount of extra book-keeping. If at any given time we keep a list of edges incident on each group, we can sample a new group from the proposal of Eq. 15 in time , provided we can sample a random neighbor also in constant time. The computation of the MH acceptance probability requires us to visit every neighbor, and hence can be done in time , where is the degree of the node being moved. The ratio can also be computed in time since the changes in Eq. 6 after the move involve a number of terms that is at most proportional to the degree of the node . This means that an entire “sweep” of the algorithm, which corresponds to single-node proposals, can be performed in time where is the number of edges. We note, however, that this is not true for all possible parametrizations of the SBM. For instance, the one used by Riolo and Newman Riolo et al. (2017) requires a number of updated terms on the order of where is the number of currently occupied groups. This number can increase algebraically with , and hence that approach will carry a performance degradation as the size of the network increases.

In Fig. 1 we show the achieved number of proposed and accepted moves per second for a variety of empirical networks collected from the KONECT repository Kunegis (2013), using an Intel Xeon Gold 6126 CPU, and with an implementation of the above algorithm using the C++ programming language Peixoto (2014c). We compare the single-node moves above with and (which amounts to fully random moves, for which the extra book-keeping is disabled), as well as a Gibbs (or “heat bath”) sampling version that corresponds to a move proposal given by


Note that this move proposal always has a MH acceptance probability of , however to implement it we need to compute the probability for the move of node to every other group, resulting in a required time for every move. We also compare with the approach by Riolo and Newman Riolo et al. (2017) that proposes single-node moves according to group sizes, using the C implementation made available by the authors. As we see in the figure, our method described above achieves millions of proposals per second, which can be sustained for networks of very large size. The extra book-keeping required by the smart move proposals incurs a barely noticeable slow-down of the rate of proposals. In fact, in certain cases it seems to speed up the algorithm, which may seem counter-intuitive. This happens because the proposal in these cases often chooses to move a node to the same group it is currently in, in which case our implementation simply skips the move altogether, saving some computation time. Both the Gibbs sampling and the Riolo-Newman method show a considerable degradation in performance as the size of the network increases, which is due to their explicit linear time dependency on the number of groups, as discussed previously. In terms of the acceptance rate, the smart move proposals achieve a roughly constant number of accepted moves in the range of to per second, independent of the size of the network. When performing fully random moves with , the performance degrades noticeably, due to the problem with a large number of groups as we described previously. Despite the ideal acceptance rate with Gibbs sampling, the actual number of acceptances obtained with it is far lower in most cases even when compared to fully random moves, which is due entirely to the higher computational cost of the move proposals. In terms of the rate of move acceptances, the performance of the Riolo-Newman algorithm is often equivalent to Gibbs sampling, even though it is based on a MH scheme. This is due to a combination of the performance degradation that comes from random moves with the extra time required to compute the likelihood ratio when using their model parametrization.

Figure 1: Number of single-node move proposals (top panel) and acceptances (bottom panel) per second, for several empirical networks, with number of nodes , edges and occupied groups indicated in the bottom axis legend, and average degree on top. The different curves correspond to the single-node moves described in the text, as well as the algorithm by Riolo and Newman Riolo et al. (2017), as show in the legend. Missing points indicate a lack of recorded events in the allotted time.

Although this kind of performance may seem encouraging at first, it is not difficult to see that ultimately this kind of approach cannot work very well. Consider, for example, starting the Markov chain from a partition consisting of a single group containing all the nodes. The only possible way forward is for one of the nodes to placed on a different group on its own. However, for a majority of networks, such move will almost always be rejected. By means of a concrete example, consider the American football network of Ref. Girvan and Newman (2002), composed of 115 nodes and 613 edges, which has a clear community structure. For this network, starting from a single group partition, the easiest node to move to a new group needs to overcome a posterior probability ratio . This means that such a move will be accepted for the first time, on average, only after around attempts, and even when this happens, the next accepted move is overwhelmingly more likely to be the one that puts this node back into its single group, moving the chain back to its initial state. The Markov chain, even though it is in theory ergodic, effectively behaves as if it were not, because it is forced to go through such states with very low probability, before reaching those that are more likely. It is clear from this example that the same problem can happen also when the starting state has more than one group, for example when two clear communities are merged into one: in order for these nodes to be split into two groups, this process must begin by placing one of the nodes into a new group by itself, which likewise will be almost always rejected. The same problem occurs also in the reverse direction: in order for an arbitrary division of a community into two groups to be merged into one, even if this final state has a much higher probability, it can only be reached by making one of the groups progressively smaller, until it is left with a single node, which then is then finally vacated. Even though this last step is likely to be accepted, the ones preceding it are progressively less probable, and would almost never be observed in any practical amount of time.

Attempts to alleviate this lack of effective ergodicity include sophisticated heuristics for initialization of the Markov chain, such as the agglomerative approach of Ref. 

Peixoto (2014a), which tries to place the initial position of the chain closer to the most likely partitions, in the hope these low probability barriers will not be relevant there. However, this assumption can be easily violated in situations where posterior distribution ascribes high probabilities to partitions with different number of groups, and an efficient algorithm would need to visit them all. In the following, we address this limitation by allowing the chain to skip over these barriers by performing the merging and splitting of groups as a single step of the Markov chain.

iii.2 Merges and splits

In order to overcome the limitations of single-node moves, we may introduce moves where entire groups are merged together. We implement this by selecting an existing group with a uniform probability and moving all its nodes to a different group with probability , yielding a transition proposal


The simplest merge choice is one where the second group is chosen uniformly at random


However, we have already seen the suboptimal acceptance rate of such random choices when doing single-node moves, and we should expect the same problems to occur when doing merges as well. Following the same line as before, a smarter merge proposal is


where is given by Eq. 14. This amounts to a random selection of a node in group , and using it to perform the same smart move proposal of group as before, but ignoring the situation where .

Such merge proposals are straightforward to implement, but in order for this to incorporated into the MH scheme we need to be able to propose “split” moves that go in the reverse direction, i.e. we need to be able to select a subset of the nodes in group and move them to a new group

. We denote this subset via a binary vector

, where if node is selected to move from group to , otherwise . This leads us to a move proposal that can be written as


A simple, but ultimately naive way of moving forward is to sample in two stages: at first we select the number of nodes to be moved uniformly at random in the interval with probability


and then sample the nodes uniformly at random, with a probability


yielding thus


Although this approach is easy to implement, it does not work in practice, because such fully random split proposals are almost never accepted, since the number of good splits, even when they exist, is exponentially outnumbered by bad ones. Furthermore, the above puts a low probability for every possible split, which means that it will also cause good merges to be rejected, as they cannot be easily reversed under this scheme.

At this point, it is important to realize that finding group splits is just a smaller scale version of our original problem of finding good overall partitions of the whole network, so we should not expect it be accomplished in one fell swoop. This inherent difficulty in proposing splits makes a merge-split MCMC less straightforward than one might realize at first. Although we might consider employing one of the many possible heuristics to find such good splits, we stumble upon two difficulties: 1. Besides sampling from it, we need also to be able to compute the probability of the proposal, in order to decide whether it can be accepted; 2. Given a merge proposal, we need to compute the corresponding split probability in the reverse direction, without actually sampling from the split proposal. These requirements limit the kind of schemes we have at our disposal. Luckily, as it has been shown before by Jain and Neal Jain and Neal (2004) in the context of Dirichlet process mixture models, we can overcome these limitations and perform splits with a much higher acceptance rate by making use of a split “staging” step as an auxiliary variable, as we describe in the following.

iii.2.1 Auxiliary variables and split staging

Let us recall the detailed balance condition that needs to be fulfilled for the Markov chain to converge to the target distribution,


Now, let us consider an augmented version of the posterior distribution obtained by sampling an arbitrary auxiliary variable with probability

. If we use MCMC to sample from the joint distribution

, we can marginalize it to obtain the original distribution , therefore sampling from this augmented space subsumes sampling from the original one. The usefulness of introducing this auxiliary variable comes from the fact that we can use it to condition our move proposals, as we will see. The detailed balance condition in this case reads


We can choose the joint transition by first making a transition and then sampling the auxiliary variable from its conditional distribution, i.e.


Then the above detailed balance condition boils down to


which, importantly, is independent of the probability . This condition is the same as the detailed balance for the original system, with the difference that the transitions are conditioned on arbitrary values of the auxiliary variable. We can incorporate this into the Metropolis-Hastings framework by conditioning our proposals , and accepting them with probability


which will enforce the condition above. The key advantage here is that we do not need to know how to compute the probability ; we need only to be able to sample from this distribution. In other words, the transition probabilities themselves can be randomly sampled, without affecting detailed balance, and without requiring us to compute their probability. This gives us freedom to implement more elaborate schemes, and leads us to the notion of split staging, where as an auxiliary variable we use a tentative split , obtained from some arbitrary distribution (which we will discuss shortly), and we perform the final split via a single Gibbs sweep from that initial state, whose probability can be easily computed as




with the shortcut notation given by


This yields a total split proposal probability


which is conditioned on the stage split , sampled from an arbitrary distribution. Note that this choice guarantees the strict reversibility of every possible merge since any split has a strict nonzero probability, regardless of the staging split .

We need now only to determine how to sample the staging split , but we can proceed in a variety of ways, since we do not need to compute the resulting probabilities. After experimentation, we determined that the following scheme works well in a variety of cases. We begin by sampling a pre-staging split uniformly at random from one of the following three algorithms:

  1. Random split:

    1. We sample the split size uniformly at random from the interval .

    2. We choose uniformly at random from all splits of size .

  2. Sequential spreading:

    1. We move all nodes in group to an empty group .

    2. In random order, we chose the nodes sequentially and move them either to group or group , with probability proportional to for a node , with being the current working partition, except for the first two moves, which are to groups and , to prevent leaving either group empty. In the end, to a node with we set , and if , we set .

  3. Sequential coalescence:

    1. We move each node in group to its own previously empty group , which in the end all have a single node each.

    2. We proceed like 2(b) above.

From the pre-stage , we obtain by performing sequential Gibbs sweeps for every node in group and , where they are allowed to move only between these two groups, forbidding moves that would leave either group empty. Note that although we could in principle compute the probability of proposing the pre-merge step, this becomes intractable after the Gibbs sweeps are performed, since we would need to compute the probability of reaching the final state via every possible intermediate trajectory. But this is precisely what is not needed in this scheme.

With the stage split in place, the final proposal is obtained by performing one more Gibbs sweep as described previously, and it is then accepted with probability


where the reverse transition corresponds to a merge proposal. Likewise a merge can be accepted with probability


which means we always need to generate a staging split for each merge candidate, using the algorithm above, to compute the reverse proposal probability.

As the number of Gibbs sweeps made during the staging step increases, the more likely it becomes that the split will be proposed with a probability proportional to the target distribution , which would be optimal. In practice, however, we do not want to choose this value too large, as it is not worth to spend too much time in a single Monte Carlo step. Although the optimal value is likely to vary for each network, we found that a value around offers a good trade-off between speed and proposal quality in most experiments we made.

iii.2.2 Joint merge-split moves

We also consider an additional kind of move that keeps the number of groups constant, and is composed of a merge of group into , which is then split again by moving some nodes from to back to , following the same scheme as before. This amounts to a proposal probability


and the final move is accepted with probability


This kind of move is able to re-distribute the nodes between two groups, without going through states with a smaller or lower number of groups, which we found to improve mixing in circumstances where the number of groups tends not to vary significantly in the posterior distribution.

iii.3 Algorithm overview and overall complexity

Even though the split and merge moves are strictly sufficient to guarantee ergodicity and detailed balance, the mixing time of the Markov chain is improved if we combine all types of moves considered above. We do so by introducing the relative move probabilities , , and , such that, e.g. the probability of a single-node move is given by


and likewise for the other kinds of moves. At each step we choose one of the moves above with the corresponding probability, and accept or reject according to the appropriate MH criterion. We note that for merge and splits we need to incorporate these move probabilities into the rejection criterion, e.g. for merge proposals we have


and likewise for splits. In our experiments, we found it more efficient to propose single nodes more often, and we have chosen and for our analysis.

If the model parametrization of Ref. Peixoto (2017) is used, the move of a single node can be done in time , where is its degree, and hence the time taken to perform a merge from group to is , and likewise, a split of a group takes a time . If we make a choice , then an entire sweep of move proposals can be done in linear time . A C++ implementation of the algorithm is freely available as part of the graph-tool Python library Peixoto (2014c).

Single-node MCMC
Merge-split MCMC
Figure 2: Comparison of runs with the single-node and merge-split MCMC algorithms, for a network of friendships between high school students Moody (2001), starting for and states (except for the single-node algorithm, where the initialization is replaced with an agglomerative heuristic) showing the values of the effective number of groups as a function of the number of sweeps (where one sweep equals move proposals), the resulting distribution of values and unnormalized posterior probability for each run, drawn with different colors.

Iv Performance for empirical networks

In Fig. 2 we show a comparison between the single-node and merge-split MCMC algorithms for a network of friendships between high school students Moody (2001). We compare six runs of each algorithm, half of which start from the single-node groups state, and the other half from the single-group

state, except for the single-node MCMC runs which would get trivially trapped in this state, in which case we initialize from an estimate of the ground state using the agglomerative heuristic described in Ref. 

Peixoto (2014a). To analyse the convergence of the chain, we compute the effective number of groups , where


is the entropy of the group membership distribution. As the results show, the single-node algorithm fails to equilibrate the chain even after an inordinate number of sweeps, as the different starting states do not converge to the same values. The initial states with a large number of groups remain at a relatively larger value, due the “merge” barriers we have described previously. More worryingly, even the improved starting states fail to make the algorithm work, as the chains clearly remain trapped in different sets of partitions. It is important to note that individual chains do not show any sign of a lack of equilibration, as they are seemingly trapped in metastable states, but they completely fail to accurately represent the target distribution. It is only when we compare different runs of the algorithm that we see a problem.

The merge-split MCMC, on the other hand, yields a substantially improved performance, where the diverging starting states converge only after a few sweeps to the same set of sampled states. The posterior distribution of effective number of groups and unnormalized posterior probability are identical between the different runs, serving as strong evidence for their successful equilibration. The pronounced multiple peaks in the distribution represent the multimodality of the posterior distribution that the single-node version of the algorithm cannot overcome.

The single-node algorithm can in fact work for networks which are sufficiently small. For example, if we apply it to the co-appearance network of characters of the novel Les Misérables by Victor Hugo Knuth (1993), we obtain essentially the same performance with both algorithms, as shown in Fig. 3. We can measure in more detail the mixing time of the Markov chain by computing the autocorrelation of a surrogate statistics, for example for we have


where is the value of after steps of the Markov chain, and is the mean value in the whole range. The value of for which approaches zero gives the decorrelation time after which we expect samples to be independent (provided the chain is not trapped in a metastable state). As seen in Fig. 3, the merge-split algorithm de-correlates faster, but the difference is very small in this example.

         Single-node          Merge-split
Figure 3: Posterior distribution of the effective number of groups obtained with the single-node and merge-split algorithms (top panel), as well as the respective autocorrelation functions (bottom panel), for the co-appearance network of characters of the novel Les Misérables by Victor Hugo Knuth (1993).
Single-node Merge-split
Network ESS ESS
Les Misérables 5,216,662 1,026,851 2,984,113 2,494,560
Football 13,036,736 2,532 5,490,028 190,267
High school 9,907,093 7 2,457,300 9,945
Table 1: Number of MCMC samples and effective sample size (ESS) obtained with the single-node and merge-split MCMC algorithms, for three empirical networks, as described in the text. Values marked with correspond to runs that failed to equilibrate even after a long time, and hence the ESS values should not be trusted.
         Single-node          Merge-split
Figure 4: Posterior distribution of the effective number of groups obtained with the single-node and merge-split algorithms (top panel), as well as the respective autocorrelation functions (bottom panel), for the network of American football games between teams Girvan and Newman (2002).
C. elegans
Network scientists
Malaria genes
Figure 5: Distribution of effective number of groups

at the lowest hierarchical level, obtained for several merge-split MCMC runs, drawn with different colors (left) and corresponding marginal hierarchical partition (right), for the neural network the

C. elegans organism White et al. (1986), the co-authorship network of network scientists Newman (2006), and the network of recombinant antigen genes from the human malaria parasite P. falciparum Larremore et al. (2013).

A related quantity that is useful to determine the actual performance of the MCMC is the effective sample size (ESS) defined as


where is the total number of samples taken from the MCMC run, and is the autocorrelation computed over only those samples. The ESS is interpreted as the effective number of uncorrelated samples taken during the run of the algorithm (provided it has equilibrated). If the chain is ideal we should have , but this is often not the case. In Table 1 we show the number of samples and the ESS values for the three networks considered previously, where the samples were gathered after several MCMC sweeps. For the Les Misérables network, the ESS with the merge-split algorithm is close to ideal, although the single-node version is still quite usable. For both the high school and football networks we see a large difference in the usefulness of both algorithms. In fact, we had already seen that for the high school network the single-node MCMC simply fails to converge, meaning that the actual ESS for that case should be in fact zero. For the football network we see in Fig 4 the autocorrelation obtained with both algorithms, as well as the posterior distribution. Not only is the decorrelation four orders of magnitude slower with the single-node algorithm, but also it in fact fails to correctly sample from all modes of the distribution, omitting values in the range , even though they do not have a negligible posterior probability. Therefore, even though this network is not very large, a detailed sampling from the posterior distribution using the single-node algorithm is already quite demanding. For network sizes slightly larger, such as the high school network, it breaks down completely, leaving the merges-split version as the only feasible alternative.

V Hierarchical partitions

Political blogs
Global airports
Figure 6: Effective number of groups at the lowest level during several long runs of the merge-split MCMC algorithm using the hierarchical SBM, drawn with different colors, and corresponding distribution, for a network of political blogs Adamic and Glance (2005) and flights between airports from

The nested SBM Peixoto (2014b, 2017) consists of a hierarchical Bayesian model, where the parameters of the SBM are interpreted as a multigraph, which are then sampled from another SBM, and so on recursively, forming a nested hierarchy. The nodes at level are the (nonempty) group labels at level , and we have a set of partitions , sampled from the posterior distribution




is the marginal likelihood (we refer to Ref. Peixoto (2017) for more details).

For this model, we can employ the same MCMC described above, by sampling from the conditional posterior distribution at each hierarchy level ,


This means we can proceed at each step in the MCMC by choosing a level uniformly at random, and setting as our target distribution, and making move proposals the exact same way as before, by performing either a single-node, merge, split, or merge-split move. The only special consideration needed is when a new group is added at a given level, a node needs to be added to the level above, together with its group membership, and possibly so on recursively. This kind of new group addition needs to be able to reverse the disappearance of a group, together with its node at the next level, and so on. The way we proceed is to keep track of unoccupied groups and their own nodes in the hierarchy, even though they do not contribute to the posterior probability, with the only purpose of being able to reverse group vacancies and to propose placement of new groups. This book-keeping needs only to be virtual, since at any given time an unoccupied group can belong to either one of the occupied groups at the level above or an unoccupied one, with equal probability. Therefore, whenever a new group is occupied, its own group membership at level is chosen uniformly at random with probability


and so on recursively for the higher levels if a new group is chosen.

In Fig. 5 we see the above algorithm employed to sample from the hierarchical partitions of a few empirical networks, in particular the neural network the C. elegans organism White et al. (1986), the co-authorship network of network scientists Newman (2006), and the network of recombinant antigen genes from the human malaria parasite P. falciparum Larremore et al. (2013). For these results we considered several chains, where the lowest level of the hierarchy was initialized both at and , and in all cases they converged to the same overall distribution of hierarchical partitions.

The convincing equilibration results we observe for the examples above are not universal, and for some data the merge-split algorithm is seemingly not sufficient to fully equilibrate the chain, at least not in a short amount of time. Some examples of this are shown in Fig. 6 for a network of political blogs Adamic and Glance (2005) and flights between airports from Although the scheme succeeds in quickly converging to the same range of partitions, with no obvious memory of the starting state, the values sampled still show noticeable discrepancy between the runs, even after a long time. This long equilibration time is not related to the hierarchical model, and also occurs when using the “flat” model for these examples. The behavior in this case seems different from the metastable trapping observed with the single-node scheme, and is probably due to a different kind of multimodality that the merge and splits are not able to efficiently evade.

Vi Conclusion

We have showed how merge and split moves can be introduced into MCMC schemes that sample network partitions from a posterior distribution, thereby significantly improving the mixing time of the Markov chain in a variety of empirically relevant scenarios. We have also demonstrated how the often employed approach of moving the group membership of a single node at a time is systematically insufficient to characterize the posterior distribution of group assignments, due to its tendency of getting trapped in metastable states.

We have extended the algorithm also to the case of hierarchical partitions, allowing us to sample nested hierarchical clusterings with the same degree of efficiency. The method developed does not depend on details of the underlying model, and can be easily extended to other variations as well.

When investigating the performance of the merge-split scheme in a variety of empirical networks we also encountered situations where, although it is far less susceptible to getting trapped as the single node scheme, it may still need very long runs to sufficiently explore the posterior landscape of partitions. A better understanding of the reasons behind the sampling difficulty in these cases may prove important to a full characterization of the modular network structure in empirical settings, and in particular to determine the quality of fit of network models used for community detection. The latter is because this kind of roughness of the posterior distribution is not expected if the network is truly sampled from the assumed generative model, and may point to a unsuitability of the model for that particular data. We leave this line of inquiry for future work.