Community detection in networks:
Modularity optimization and maximum likelihood are equivalent

Modularity maximization is perhaps the most widely used method for community detection for networks. It operates by defining a benefit function, called the modularity, that measures the quality of divisions of a network into communities. One optimizes this benefit function over possible divisions of the network of interest to find the one that gives the highest score, taking this to be the definitive division of the network. Since the number of possible divisions of a network is exponentially large, we normally cannot perform the optimization exhaustively, so we turn instead to approximate optimization methods, of which many have been tried, including greedy algorithms Newman04a ; CNM04 , extremal optimization DA05 , spectral relaxation Newman06b , genetic algorithms SYHF09 , simulated annealing GSA04 ; MAD05 , and belief propagation ZM14 . The popular Louvain algorithm for community detection BGLL08 , which is built into a number of network analysis software packages, uses a multiscale modularity optimization scheme and is one of the fastest community detection methods in practice.

The definition of the modularity function is straightforward NG04 . We desire a benefit function which, given a network and a candidate division of that network into groups, returns a score that is larger if the division is a “good” one and smaller if it is “bad.” The heuristic notion used to define the modularity is that a good division is one that places most of the edges of a network within groups and only a few of them between groups.

Let us represent our network by its adjacency matrix. For an undirected network of $n$ nodes the adjacency matrix $\mathbf{A}$ is the real symmetric $n\times n$ matrix with elements $A_{ij}=1$ if there is an edge between nodes $i$ and $j$ and 0 otherwise. Further, let us consider a division of the network into $q$ nonoverlapping groups, numbered (in any order) from 1 to $q$, and let us denote by $g_{i}$ the number of the group to which node $i$ is assigned. Thus the complete vector $\mathbf{g}$ of group assignments specifies the division of the network. Then the number of edges that fall within groups, for this particular division, is equal to $\tfrac{1}{2}\sum_{ij}A_{ij}\delta_{g_{i}g_{j}}$, where $\delta_{ij}$ is the Kronecker delta and the leading factor of a half prevents double counting of edges.

The number of in-group edges alone, however, is not a good measure of the quality of a division, since it can be trivially maximized by putting all the nodes in one of the $q$ groups and none in any of the others. This would put 100% of edges inside groups but clearly doesn’t constitute a useful division of the network. Instead, therefore, modularity measures not just the number of edges within groups but the difference between that number and the expected number of such edges, were edges placed at random within the network.

Suppose we take our observed network and randomize the positions of its edges. We keep the total number of edges the same but we reposition them between the nodes at random, in a manner to be determined shortly. And suppose that, following this randomization, the probability that nodes $i$ and $j$ are connected by an edge is $P_{ij}$. Then the expected number of edges within groups, post-randomization, is $\tfrac{1}{2}\sum_{ij}P_{ij}\delta_{g_{i}g_{j}}$. The modularity is then proportional to the actual minus expected number of edges thus:

where $m$ is the total number of edges in the network and is included here by convention only—it makes $Q$ equal to a fraction of edges rather than an absolute number, which makes modularity values more easily comparable between networks of different sizes. For the purposes of maximizing the modularity, which is our main concern here, the factor of $m$ makes no difference at all. The position of the maximum does not depend on overall constant factors.

Note that if we now put all nodes in the same group, then $\delta_{g_{i}g_{j}}=1$ for all $i,j$ and

since, as we have said, the number of edges in the network is held constant during randomization, and hence $\sum_{ij}P_{ij}=\sum_{ij}A_{ij}=2m$. Thus we no longer get a high modularity score for putting all nodes in a single group together. The maximum of modularity occurs for some other (nontrivial) division of the nodes, which we take to be the best division of the network. This is the method of modularity maximization.

It remains to determine what $P_{ij}$ is. The value depends on the particular scheme we use to randomize the positions of the edges. The simplest scheme would be just to reposition the edges uniformly at random, every position having the same probability as every other. For a network of $n$ nodes there are ${n\choose 2}$ places to put an edge, and hence the probability of filling any of them with one of the $m$ edges is

(Technically this is the expected number of edges not the probability, but normally $m\ll{n\choose 2}$ so that probability and expected number are essentially the same.)

In practice, however, this choice does not work very well because it fails to respect the degrees of the nodes in the network. The probabilities of connections between nodes depend strongly on the total number of connections nodes have—their degrees—with nodes of high degree being much more likely to be connected than nodes of low degree. For reasons that will become clear in this paper, it is important to include this effect in the definition of modularity if things are to work correctly.

Instead, therefore, we consider a constrained randomization of the edges in the network in which we preserve the degree of every node, but otherwise position the edges at random. This kind of randomization is well known in the study of networks: it gives rise to the random graph ensemble known as the configuration model MR95 ; NSW01 . After randomization, the probability of connection between two nodes is equal to

where $k_{i}=\sum_{j}A_{ij}$ is the degree of node $i$ and $m$ is once again the number of edges in the network. (Again, this is technically the expected number of edges, but the probability and expected number are essentially the same.)

This is the choice that is most commonly used in the definition of the modularity. With this choice the modularity is given by

which is the form in which it is most commonly written.

There is a further twist, however, because even this definition does not always work well. As shown by Fortunato and Barthélémy FB07 , community detection by modularity maximization using the definition of (5), while it works in many situations, has one specific shortcoming: it is unable to find community structure in large networks with many small communities. In particular, if the number of communities in a network is greater than about $\sqrt{2m}$, then the maximum modularity will not correspond to the correct division. The maximum will instead tend to combine communities into larger groups and fail to resolve the smallest divisions in the network.

To address this problem, Reichardt and Bornholdt RB06a proposed a generalized modularity function

When the parameter $\gamma=1$, this is the same as the traditional modularity of Eq. (5), but other choices allow us to vary the relative weight given to the observed and randomized edge terms. If one places more weight on the observed edge term (by making $\gamma$ smaller), the maximum modularity division favors, and the method therefore tends to find, larger communities. If one places more weight on the randomized edge term (larger $\gamma$), the method finds smaller communities. There has not previously been any fundamental theory dictating what value of $\gamma$ one should use, but this is one of the questions on which we will shed light in this paper.

The other method of community detection we consider is the method of statistical inference, as applied to the stochastic block model and its variants. With this method, one fits a generative model of a network to observed network data, and the parameters of the fit tell us about the structure of the network in much the same way that fitting a straight line through a set of data points tells us about their slope.

The model most commonly used in this context is the so-called stochastic block model, which is a random graph model of a network with community structure HLL83 ; NS01 ; BC09 . One takes some number $n$ of nodes, initially without any edges, and divides them into $q$ groups in some way, with $g_{i}$ being the group to which node $i$ is assigned, as previously. Then one places edges between nodes independently at random, with the probability, which we denote $\omega_{rs}$, of an edge between a particular pair of nodes depending on the groups $r$ and $s$ to which the nodes belong. Thus there is a symmetric $q\times q$ matrix of parameters $\omega_{rs}$ which determine the probabilities of edges within and between every pair of groups. If the diagonal elements $\omega_{rr}$ of this matrix are larger than the off-diagonal elements, then networks generated by the model have a higher probability of edges within groups than between them and hence have traditional community structure.

In fact, the stochastic block model is often studied in a slightly different formulation in which one places not just a single edge between any pair of nodes but a Poisson distributed number of edges with mean $\omega_{rs}$. Thus $\omega_{rs}$ is the expected number, rather than the probability, of edges between nodes in groups $r$ and $s$, and the networks generated by the model can in principle have multiedges, meaning there can be more than one edge between the same pair of nodes. Moreover, one typically also allows the network to contain self-edges, edges that connect a node to itself, which are also Poisson distributed in number, with mean $\tfrac{1}{2}\omega_{rr}$ for a node in group $r$. (The factor of half is included solely because it makes the algebra simpler.) The inclusion of multiedges and self-edges in the model is unrealistic in the sense that most network data encountered in the real world contain neither. However, most real-world networks are also very sparse, meaning that only a tiny fraction of all possible edges that could exist actually do, and hence the values of the edge probabilities $\omega_{rs}$ are very small. In this situation, the density of multiedges and self-edges in the network will itself be small and the Poisson version of the model is virtually indistinguishable from the first (Bernoulli) version defined above. At the same time, the Poisson version is technically easier to handle than the Bernoulli version. In this paper we use the Poisson version. Similar results can be derived for the Bernoulli version, but the formulas are more complicated and the end result is little different.

The definition of the model above is in terms of its use to generate networks. As applied to community detection, however, the model is used in the “reverse” direction to infer structure by fitting it to data. In this context, one hypothesizes that an observed network, with adjacency matrix $\mathbf{A}$, was generated from the stochastic block model, and attempts to work out what values of the model parameters must have been used in the generation. The parameters in this case are the edge probabilities $\omega_{rs}$ and the group memberships $g_{i}$.

Given particular values of the parameters we can write down the probability, or likelihood, that the observed network was generated from the block model thus:

where $\mathbf{\Omega}$ denotes the complete matrix of values $\omega_{rs}$ and we have adopted the standard convention that a self-edge is represented by a diagonal adjacency matrix element $A_{ii}=2$ (and not 1 as one might at first imagine).

The position of the maximum of this quantity with respect to $\mathbf{\Omega}$ and $\mathbf{g}$ tells us the values of the parameters most likely to have generated the observed network. Typically we are interested only in the group assignments $\mathbf{g}$, which tell us how the network divides into groups. The probabilities $\mathbf{\Omega}$ are usually uninteresting and can be discarded once the likelihood has been maximized.

Alternatively (and usually more conveniently), we can maximize the logarithm of the likelihood:

The terms $\tfrac{1}{2}A_{ii}\log\tfrac{1}{2}$, $\log(\tfrac{1}{2}A_{ij}!)$, and $\log A_{ij}!$ are all independent of the parameters and do not affect the position of the maximum, so they can be ignored, and the log-likelihood then simplifies to

The optimal division of the network into communities is then given by maximizing this quantity with respect to both $\mathbf{g}$ and $\mathbf{\Omega}$.

As with the modularity, however, this is not the whole story. This approach fares poorly when applied to most real-world networks because it doesn’t respect the node degrees in the network. The stochastic block model as described here (in either Bernoulli or Poisson versions) generates networks that have a Poisson degree distribution, which is very different from the broad distributions seen in empirical networks. This means that, typically, the model does not fit observed networks well for any choice of parameter values. It’s as if one were trying to fit a straight line through an inherently curved set of data points. Even the best fit of such a line will not be a good fit. There are no good fits when the model you are fitting is simply wrong.

The conventional solution to this problem, in the present situation, is to use a slightly different model, the degree-corrected block model KN11a , which can fit networks with any degree distribution. In this model the nodes are again assigned to groups $g_{i}$ and edges placed independently at random between them, but now the expected number of edges between nodes $i$ and $j$ is $(k_{i}k_{j}/2m)\omega_{rs}$ (where as before $r,s$ are the groups to which the nodes belong and $k_{i},k_{j}$ are the degrees) or a half that number for self-edges. The factor of $2m$ in the denominator is optional, but convenient since, as we have said, $k_{i}k_{j}/2m$ is the probability of an edge in the configuration model and hence, with this definition, $\omega_{rs}$ quantifies the probability of edges relative to the configuration model.

Following the same line of reasoning as before, and again neglecting constants that have no effect on the position of the likelihood maximum, the log-likelihood for this model is

Community detection now involves the maximization of this quantity with respect to the parameters $\mathbf{\Omega},\mathbf{g}$ to find the best fit of the model to the observed network. This maximization can be achieved in a number of ways. As with the modularity, there are too many possible group assignments $\mathbf{g}$ to maximize exhaustively on any but the smallest of networks, but researchers have successfully applied a variety of approximate methods, including label switching algorithms BC09 , Kernighan–Lin style greedy algorithms KL70 ; KN11a , spectral methods Newman13b , Monte Carlo NS01 ; Peixoto14b , and belief propagation DKMZ11a ; Yan14 .

A drawback of the equivalence we have demonstrated is that it applies only when we use the correct value of the resolution parameter $\gamma$, which normally we do not know. One can, however, make an empirical estimate of the value of $\gamma$ using an iterative scheme as follows.

First, one makes an initial guess about the value of $\gamma$. This guess need not be particularly accurate: $\gamma=1$ usually works fine. Then one finds the communities in the network by modularity maximization using this choice. This gives us some set of assignments $g_{i}$ of nodes to groups—likely not optimal—from which we can then make an estimate of the parameters $\omega_{\textrm{in}}$ and $\omega_{\textrm{out}}$ by noting that the expected total number of in-group edges $m_{\textrm{in}}$ in the (degree-corrected) planted partition model is

where $\kappa_{r}=\sum_{i}k_{i}\delta_{g_{i},r}$ is the sum of the degrees of the nodes in group $r$. Hence we can estimate $\omega_{\textrm{in}}$ from

using the observed value of $m_{\textrm{in}}$ as an estimate of the expected value. Similarly for $\omega_{\textrm{out}}$ we have

where $m_{\textrm{out}}$ is the number of edges running between distinct groups.

Using these estimates of $\omega_{\textrm{in}}$ and $\omega_{\textrm{out}}$ we can now calculate a new value of $\gamma$ from Eq. (15). Then we repeat the process, maximizing the modularity and recalculating $\omega_{\textrm{in}}$, $\omega_{\textrm{out}}$, and $\gamma$ until we achieve convergence. The consistency of modularity maximization, mentioned earlier, implies that this procedure should converge to the correct value of $\gamma$ (and the correct community structure) for sufficiently dense networks that actually are generated from the planted partition model. For all other networks (meaning, in practice, for all real-world applications of the method) we have no formal guarantees of correctness or convergence, though the same is also true of all other methods of community detection, including, but not limited to, community detection by statistical inference.

One might imagine that this would not be a very efficient method for calculating $\gamma$: it requires repeated maximization of different modularity functions until the correct value of $\gamma$ is reached. In practice, however, we have found that it converges very quickly. In most cases we have examined, $\gamma$ is calculated to within a few percent after just one iteration, and in no case have we found a need for more than ten iterations, so the method may in fact be quite serviceable. Figure 1 shows an example application to a set of artificially generated (“synthetic”) networks for which the true value of $\gamma$ is known and, as the figure shows, the algorithm is able to determine that value accurately in every case. Table 1 gives values of $\gamma$ calculated using the algorithm for a number of real-world networks that have been used as test cases in previous community detection studies.

The values of $\gamma$ vary in size, but there is an overall trend towards larger values in networks with larger numbers of communities, both among the synthetic networks and the real ones. This is perhaps not unexpected given that the resolution parameter $\gamma$ was originally introduced precisely in order to deal with networks with larger numbers of communities. Recall that larger values of $\gamma$, and specifically values larger than the traditional value of 1, are needed in networks where the number of communities exceeds the resolution limit at $\sqrt{2m}$. None of the networks studied here approach this limit, but nonetheless we should not find it surprising that the larger values of $q$ in both Fig. 1 and Table 1 are best treated using values $\gamma>1$.

Whether the algorithm given here is in fact a useful one in practice is a debatable point. As we have shown, it does no more than the likelihood maximization method, and the latter in principle gives better results, since it does not assume that all groups are statistically similar. Modularity maximization does have the advantage of being less nonlinear than maximum-likelihood methods, which allows for some faster algorithms such as spectral Newman06b ; Newman13b and multiscale BGLL08 algorithms. Still, the results derived here are primarily of interest not because of the algorithms they suggest, but because of the light they shed on the strengths and weaknesses of modularity maximization.