# Optimizing Edge Sets in Networks to Produce Ground Truth Communities Based
on Modularity

Daniel Kosmas  Department of Mathematics, Rensselaer Polytechnic Institute
John E. Mitchell  Department of Mathematics, Rensselaer Polytechnic Institute
Thomas C. Sharkey  Department of Industrial Engineering, Clemson University
Boleslaw K. Szymanski  Department of Computer Science, Rensselaer Polytechnic
Institute

###### Abstract

We consider two new problems regarding the impact of edge addition or removal
on the modularity of partitions (or community structures) in a network. The
first problem seeks to add edges to enforce that a desired partition is the
partition that maximizes modularity. The second problem seeks to find the
sparsest representation of a network that has the same partition with maximum
modularity as the original network. We present integer programming
formulations, a row generation algorithm, and heuristic algorithms to solve
these problems. Further, we demonstrate a counter-intuitive behavior of
modularity that make the development of heuristics for general networks
difficult. We then present results on a selection of social and illicit
networks from the literature.

##  1 Introduction

A well studied problem in network science is community detection. Given a
network $G=(V,E)$, where $m$ is the number of edges and $d_{i}$ is the degree
of node $i$, the community detection problem seeks to partition the node set
into disjoint clusters (or communities), $C_{1},\ldots,C_{k}$, satisfying
$\bigcup_{i=1}^{k}C_{i}=V$ and $C_{i}\cap C_{j}=\emptyset$ for all $i\neq j$,
in order to maximize some measure of this partition. Typically, this measure
would give preference to partitions where there are a large number of edges
within each cluster, but a small number of edges between the clusters. This
problem has applications in many domains, including biological networks [1],
social networks [1], and illicit network analysis [2].

A popular method to perform community detection is through modularity
maximization. Introduced by Newman and Girvan [3], modularity maximization
seeks to find the partition $P$ that maximizes the Newman-Girvan modularity
metric, which compares, for each cluster, its number of actual edges and
expected numbers of edges. Brandes et. al. [4] provides an equivalent
definition of modularity as

$Q=\sum_{C\in P}\left(\frac{|E(C)|}{m}-\left(\frac{|E(C)|+\sum_{C^{\prime}\in
P}|E(C,C^{\prime})|}{2m}\right)^{2}\right),$

where $E(C)$ is the number of edges within cluster $C$ and $E(C,C^{\prime})$
is the number of edges between clusters $C$ and $C^{\prime}$. The problem of
modularity maximization is known to be NP-hard [4].

In some scenarios, we may know the community structure of a network before all
of its edges are identified. However, community detection methods will find
the optimal partition in this partial network, and that optimal partition may
or may not agree with the true community structure since not all edges are
currently identified. Edge recommendation in social networks is a well-studied
problem [5], but it has not been studied when trying to enforce that a fixed
partition should maximize modularity across all potential partitions. In this
paper, we introduce the problem of finding the minimal number of edges to add
in order to enforce that the ground truth community structure has the maximum
modularity across all partitions.

Our analysis of this problem is motivated by problems in understanding and
disrupting illicit networks. In disrupting these networks, it is important to
predict how the network will react to disruptions since interdicting all
participants may not be possible. In such cases, the ‘ground truth’
communities of the participants would be known prior to the disruptions and
could have been obtained through intelligence-gathering activities on these
illicit networks. We then have the ‘remaining’ network after disruptions have
removed participants from the illicit network. The illicit network will then
react by adding edges into this network to restore the ground truth
communities and the problem we consider helps to understand the fewest number
of edges required to accomplish such restoration.

Additionally, it may be of interest to identify the so-called ‘backbone’ of
the network that underlies its current partition. This problem can be framed
as identifying a sparse sub-network which maintains the partition of the
original network. In other words, we seek to find the smallest number of edges
(amongst the current ones) such that the partition that maximizes the
modularity in the sub-network is the same as the partition that maximizes
modularity in the original network. This problem is important when edges are
suspected to exist in the network due to partial information and resources
must be allocated to confirm their presence. This again finds applications in
disrupting illicit supply networks as initial information may cause
‘suspected’ edges to appear and disruptions may not be implemented until the
underlying network structure is confirmed. Understanding the minimum number of
edges to produce the ground truth communities will allow scarce resources to
be allocated to best investigate the edges of the network.

We present approaches to formulate these problems as mixed integer nonlinear
programs (MINLPs). For the case where edges are unweighted, we are able to
formulate these problems as binary nonlinear programs. We then show how to
reformulate the binary nonlinear program as an equivalent binary linear
program. For the case where edges are weighted, we show how the MINLP can be
relaxed to a mixed integer linear program (MILP) and also discuss exact
approaches to solve this problem. We additionally develop heuristic approaches
to solve these problems, and demonstrate properties about modularity that make
the development of more accurate heuristics difficult. We apply our approaches
to several case studies as a proof of concept and show our heuristics perform
quite well.

###  1.1 Literature Review

Community detection has been a problem of interest in network analysis for
many years, with many different methods [1]. One such method, modularity
maximization, has been particularly well studied [6]. The goal of modularity
maximization is to find the partition that maximizes the Newman-Girvan
modularity metric, a global metric of community structure. Many heuristic
methods have been developed to quickly find partitions with relatively large
modularity scores. Similar metrics have also been proposed to capture the same
intuitive principles guiding modularity, while accounting for some of its
shortcomings.

Modularity maximization can be formulated as an integer programming (IP)
problem. Brandes et. al. [4] introduces this formulation, and Agarwal and
Kempe [7] propose a rounding technique to improve the solve time of these
problems. Let $A$ be the adjacency matrix of $G$, and for every pair of nodes
$i,j\in V$, let $x_{ij}=0$ if $i$ and $j$ are in the same cluster, and
$x_{ij}=1$ otherwise. The integer programming formulation of modularity
maximization is as follows:

$\displaystyle\max_{x\in\\{0,1\\}^{|V|^{2}}}~{}~{}~{}$
$\displaystyle\frac{1}{2m}\sum_{i,j\in
V}\left(A_{ij}-\frac{d_{i}d_{j}}{2m}\right)(1-x_{ij})$ s.t. $\displaystyle
x_{ij}\leq x_{ik}+x_{kj}$ $\displaystyle\text{ for all }i,j,k\in V$ (1)

There have since been other formulations of this problem to improve solve time
of the integer program (see [8] for such methods), but discussing them is
outside the scope of our paper. Such methods can easily be adapted to fit into
our framework.

One issue that arises in our problems, compared to this past work on using IP
to maximize modularity, is the fact that we are enforcing that the modularity
of the given partition is maximum amongst all partitions. This leads to a
number of constraints equal to the number of partitions, i.e., an exponential
number of constraints based on the number of nodes in the network. To avoid
this difficulty, our framework operates on ‘row generation,’ i.e., only
maintaining a subset of these constraints at a time and then generating new
constraints (rows) when needed. In each iteration of our framework, a master
problem is solved for a potential optimal solution based on a subset of the
partition constraints. A subproblem is then solved to verify optimality, which
fails when there is a partition not currently included that has a modularity
larger than the ground truth community structure. In such a case, a new
constraint is added to the master problem. Fischetti et. al. [9] demonstrates
this is an effective method to solve bilevel integer optimization problems.
More recently, community detection through the degree-corrected stochastic
block model, which is known to be equivalent to modularity maximization [10],
has been modeled as a MINLP in [11]. The authors linearized the model and
solved the resulting MILP using a row generation framework. However, our
problem has the added complication of needing to compare the ground truth to
all potential partitions in a graph where the problem is deciding the
existence of edges.

###  1.2 Paper Organization

The rest of the paper is organized as follows: in Section 2, our two problems
of interest are formally introduced. In Section 3, we describe a row
generation algorithm to solve the problem exactly, and reformulate the problem
as a MILP. In Section 4, we demonstrate properties of modularity that prove
challenging when designing heuristic solution methods. We also present our
heuristic approaches. In Section 5, we present our computational results of
our row generation algorithm and heuristics. We discuss our conclusions in
Section 6.

##  2 Description of Problems

###  2.1 Enforcing Ground Truth Communities as an Optimal Partition

Consider an unweighted, but not necessarily simple, network $G=(V,E)$ and
partition $T=(C_{1}^{T},C_{2}^{T},\ldots,C_{l}^{T})$, and let $K_{|V|}$ be the
set of all edges on $|V|$ nodes and $\mathcal{P}$ be the set of all possible
partitions. We seek to find the minimum set of edges $E^{*}\subset
K_{|V|}\setminus E$ to add to the network such that $T$ is the partition that
maximizes modularity in $G^{\prime}=(V,E\cup E^{*})$. Let $Q_{P}(E^{\prime})$
be the modularity of partition $P$ in the network $(V,E\cup E^{\prime})$. We
can define this problem as a mathematical program as follows:

$\displaystyle\min_{E^{\prime}\subset K_{|V|}\setminus E}~{}~{}~{}$
$\displaystyle|E^{\prime}|$ s.t. $\displaystyle
T\in\text{arg}\max_{P\in\mathcal{P}}\\{Q_{P}(E^{\prime})\\}$ (2)

By enumerating over all possible partitions, we can easily convert the
constraint in model (2.1) to the following:

$\displaystyle\min_{E^{\prime}\subset K_{|V|}\setminus E}~{}~{}~{}$
$\displaystyle|E^{\prime}|$ s.t. $\displaystyle Q_{T}(E^{\prime})\geq
Q_{P}(E^{\prime})$ $\displaystyle\text{ for all }P\in\mathcal{P}$ (3)

In this program, the constraint enforces that the modularity associated with
partition $T$ after adding the set of edges $E^{\prime}$ is at least as large
as the modularity of every other partition $P$ after adding the set of edges
$E^{\prime}$. Thus, $T$ is in the set of optimal partitions for the modularity
maximization problem. To enforce that $T$ is the unique optimal solution, the
constraint could be replaced with $Q_{\mathcal{T}}(E^{\prime})\geq
Q_{\mathcal{P}}(E^{\prime})+\epsilon$ for some small $\epsilon>0$.

We can further extend this problem for weighted networks, where we seek to
find the minimum set of weights to add such that the partition $T$ is the
partition that maximizes the modularity in the weighted networks. Consider a
weighted network $G=(V,E)$ with weights $w\in\mathbb{R}_{\geq 0}^{|E|}$. Let
$u\in\mathbb{R}_{\geq 0}^{K_{|V|}}$ be the bound on weights on the edges, and
$a\in\mathbb{R}_{\geq 0}^{K_{|V|}}$ be the amount of weight added to the
edges. Unlike before, $Q_{P}(E^{\prime})$ should be the modularity of
partition $P$ in the network $(V,E\cup E^{\prime})$ with edge weights $w+a$.

$\displaystyle\min_{E^{\prime}\subset K_{|V|}}~{}~{}~{}$
$\displaystyle\sum_{e\in K_{|V|}}a_{e}$ s.t. $\displaystyle
Q_{T}(E^{\prime})\geq Q_{P}(E^{\prime})$ $\displaystyle\text{ for all
}P\in\mathcal{P}$ (4) $\displaystyle w_{e}+a_{e}\leq u_{e}$
$\displaystyle\text{ for all }e\in K_{|V|}$ $\displaystyle a\geq 0$

###  2.2 Sparsest Representation of a Partition

We now introduce the problem of finding the minimum set of edges in a network
such that the sub-network composed of only these edges has the same optimal
partition as the original network. Consider a network $G=(V,E)$ with optimal
partition $T=P^{*}$. We can express this as the following mathematical
program:

$\displaystyle\min_{E^{\prime}\subseteq E}~{}~{}~{}$
$\displaystyle|E^{\prime}|$ s.t. $\displaystyle Q_{P^{*}}(E^{\prime})\geq
Q_{\mathcal{P}}(E^{\prime})$ $\displaystyle\text{ for all }P\in\mathcal{P}$
(5)

Unlike model (2.1), where we optimize over $K_{|V|}\setminus E$, here we
optimize over $E$. Because the constraints in both models are the same, we can
use a similar solution method to solve both problems.

##  3 Integer Programming Framework

###  3.1 Initial Framework

To solve model (2.1) by enumeration, we would need to verify the modularity of
every possible partition, the number of which is exponentially large in the
number of nodes. However, we expect that most of those partitions can be
rejected based on modularity properties. For example, if every node is in the
same cluster, or every node is in their own cluster, the modularity will
always be zero. As long as our given partition has a non-negative modularity,
we do not need to compare against these two partitions.

We propose a row generation method to solve these problems. In this framework,
we iteratively identify partitions of higher modularity than our given
modularity. We then use these partitions to determine what edges should be
added to the network so the given partition has a higher modularity than these
partitions.

Let $\mathcal{P}^{k}$ be the set of partitions identified by the $k^{th}$
iteration of our iterative procedure. Instead of solving model (2.1), we solve
the following:

$\displaystyle\min_{E^{\prime}\subset K_{|V|}\setminus E}~{}~{}~{}$
$\displaystyle|E^{\prime}|$ s.t. $\displaystyle Q_{T}(E^{\prime})\geq
Q_{\mathcal{P}}(E^{\prime})$ $\displaystyle\text{ for all
}P\in\mathcal{P}^{k}$ (6)

In model (3.1), we enforce that $T$ has a modularity at least as large as
every partition in $\mathcal{P}^{k}$. However, this does not enforce that the
given partition $T$ is in the set of partitions that maximize modularity. Let
$E^{k}$ be the optimal solution found in model (3.1). We now find the
partition that maximizes modularity in the new network with edges added. Let
$\bar{A}$ be the adjacency matrix of $G^{k}=(V,E\cup E^{k})$, $\bar{m}=|E\cup
E^{k}|$, and $\bar{d}_{i}$ be the degree of node $i$ in $G^{k}$. We find the
partition that maximizes modularity with the following:

$\displaystyle\max_{x\in\\{0,1\\}^{|V|^{2}}}~{}~{}~{}$
$\displaystyle\frac{1}{2\bar{m}}\sum_{i,j}\left(\bar{A}_{ij}-\frac{\bar{d}_{i}\bar{d}_{j}}{2\bar{m}}\right)(1-x_{ij})$
s.t. $\displaystyle x_{ij}\leq x_{ik}+x_{kj}$ $\displaystyle\text{ for all
}i,j,k\in V$ (7)

By solving model (3.1), we find the partition $\bar{P}$ that maximized
modularity in the new network, as well as its modularity. If
$Q_{T}(E^{k})=Q_{\bar{P}}(E^{k})$, then we have that $T$ is in the set of
optimal partitions, and we are done. If not, then we have found a partition,
$\bar{P}$, which has greater modularity than $T$ which must be considered when
adding edges. We set $\mathcal{P}^{k+1}=\mathcal{P}^{k}\cup\bar{P}$, and
repeat this process. Formalizing this procedure yields Algorithm 1.

Algorithm 1 Row Generation for Edge Addition

Initialize: network $G=(V,E)$, ground truth partition $T$, initial set of
partitions $\mathcal{P}^{1}$, iteration counter $k=1$.

while not converged do

Step 1. Solve model (3.1) for $E^{k}$.

if (3.1) is infeasible then

terminate; model (2.1) is infeasible.

end if

Step 2. Create test network $G^{k}=(V,E\cup E^{k})$. With test network
$G^{k}$, solve model (3.1) for optimal partition $\bar{P}$.

if $Q_{T}(E^{k})=Q_{\bar{P}}(E^{k})$ then

terminate; $E^{k}$ is the optimal solution to model (2.1).

else

Set $\mathcal{P}^{k+1}=\mathcal{P}^{k}\cup\bar{P}$, $k=k+1$, and return to
Step 1.

end if

end while

To analyze, each iteration of this procedure can result in one of three
outcomes. First, if model (3.1) is infeasible, then model (2.1) is also
infeasible. If model (3.1) is feasible, we have two cases. If
$Q_{T}(A)=Q_{\bar{P}}(A)$ in the $k^{th}$ iteration, we have found the optimal
set of edges to add to the network. Otherwise, we have
$Q_{T}(A)<Q_{\bar{P}}(A)$, and we continue the process. Since $|\mathcal{P}|$
is finite, our framework has a finite number of iterations, and will thus
solve model (2.1) in finite time.

While we choose to use modularity as our measure of community structure, note
that our framework can work for any optimization based clustering measure. In
model (2.1), and likewise model (3.1), the constraint enforcing that the given
partition $T$ has modularity as least as large as all other partitions would
be replaced with the desired measure of $T$ being at least as large as all
other partitions, and model (3.1) would be replaced with finding the partition
that optimizes that measure. In this paper, though, we are able to express
model (2.1) as an integer program, rather than a nonlinear integer program,
which we demonstrate over the next two subsections.

###  3.2 Nonlinear Integer Programming Representation

We now express model (2.1) as an integer program. Let $z_{ij}$ be the
indicator of whether or not the edge $(i,j)$ is included. We can rewrite the
modularity of a given partition with added edges as:

$Q_{P}(A)=\sum_{(i,j)\in
V}\left(\frac{A_{ij}+z_{ij}}{2(m+\sum_{(i,j)\in\bar{E}}z_{ij})}-\frac{(d_{i}+\sum_{l\in
V}z_{il})(d_{j}+\sum_{h\in
V}z_{jh})}{(2(m+\sum_{(s,t)\in\bar{E}}z_{st}))^{2}}\right)(1-x^{P}_{ij})\text{.}$
(8)

We can thus express the nonlinear integer program as:

$\displaystyle\min_{z\in\\{0,1\\}^{K_{|V|}\setminus E}}$
$\displaystyle\sum_{(i,j)\in K_{|V|}\setminus E}z_{ij}$ s.t.
$\displaystyle\sum_{(i,j)\in
V}\left(\frac{A_{ij}+z_{ij}}{2(m+\sum_{(i,j)\in\bar{E}}z_{ij})}-\frac{(d_{i}+\sum_{l\in
V}z_{il})(d_{j}+\sum_{h\in
V}z_{jh})}{(2(m+\sum_{(s,t)\in\bar{E}}z_{st}))^{2}}\right)(1-x^{T}_{ij})$ (9)
$\displaystyle\geq\sum_{(i,j)\in
V}\left(\frac{A_{ij}+z_{ij}}{2(m+\sum_{(i,j)\in\bar{E}}z_{ij})}-\frac{(d_{i}+\sum_{l\in
V}z_{il})(d_{j}+\sum_{h\in
V}z_{jh})}{(2(m+\sum_{(s,t)\in\bar{E}}z_{st}))^{2}}\right)(1-x^{P}_{ij})$
$\displaystyle\text{ for all }P\in\mathcal{P}$

We note that model (2.2) can be expressed similarly, where $z_{ij}$ is no
longer a binary vector indicating inclusion of edges, but a vector of
integers, expressing the amount of weight added to edges.

###  3.3 Linearizing the Problem

The constraints in model (3.2) is nonlinear, but can be linearized quite
nicely using standard techniques. Below, we describe the process of
linearizing these constraints. The equivalent constraint in the weighted
version of the problem can be linearized similarly, but we will highlight some
key differences after the derivation is complete.

The first point to note is that the $2(m+\sum_{(i,j)\in\bar{E}}z_{ij})$ in the
denominator is the same on both sides of the equality. Furthermore, it will
always be strictly greater than $0$, assuming the network has a nonempty edge
set. Thus, multiplying both sides of the inequality by this denominator does
not change the inequality, resulting in:

$\displaystyle\sum_{(i,j)\in V}\left(A_{ij}+z_{ij}-\frac{(d_{i}+\sum_{l\in
V}z_{il})(d_{j}+\sum_{h\in
V}z_{jh})}{2(m+\sum_{(s,t)\in\bar{E}}z_{st})}\right)(1-x^{T}_{ij})$ (10)
$\displaystyle\geq\sum_{(i,j)\in V}\left(A_{ij}+z_{ij}-\frac{(d_{i}+\sum_{l\in
V}z_{il})(d_{j}+\sum_{h\in
V}z_{jh})}{2(m+\sum_{(s,t)\in\bar{E}}z_{st})}\right)(1-x^{P}_{ij})\text{.}$

Again, using that $2(m+\sum_{(i,j)\in\bar{E}}z_{ij})>0$, we multiply both
sides of the inequality by $2(m+\sum_{(i,j)\in\bar{E}}z_{ij})$ to remove
fractions from our inequality. Additionally, the factors $1-x_{ij}^{T}$ and
$1-x_{ij}^{P}$ are multiplied by are the same terms. Thus, we can reduce the
inequality to the following:

$\displaystyle\sum_{i,j\in
V}\left((A_{ij}+z_{ij})(2(m+\sum_{(s,t)\in\bar{E}}z_{st}))-(d_{i}+\sum_{l\in
V}z_{il})(d_{j}+\sum_{h\in V}z_{jh})\right)x^{T}_{ij}$ (11)
$\displaystyle\leq\sum_{i,j\in
V}\left((A_{ij}+z_{ij})(2(m+\sum_{(s,t)\in\bar{E}}z_{st}))-(d_{i}+\sum_{l\in
V}z_{il})(d_{j}+\sum_{h\in V}z_{jh})\right)x^{P}_{ij}$

We now have a quadratic constraint equivalent to the original constraint. To
linearize this constraint, we use the McCormick inequalities [12]. For a pair
of edges $(i,j)$ and $(s,t)$, let $w_{ijst}=z_{ij}z_{st}$. The following
inequalities enforce that $w_{ijuv}$ takes the correct value:

$\displaystyle w_{ijst}$ $\displaystyle\leq z_{ij}\text{,}$ $\displaystyle
w_{ijst}$ $\displaystyle\leq z_{st}\text{,}$ $\displaystyle w_{ijst}$
$\displaystyle\geq 0\text{,}$ $\displaystyle w_{ijst}$ $\displaystyle\geq
z_{ij}+z_{st}-1\text{.}$

By multiplying out the quadratic terms, and replacing any bilinear terms with
the appropriate variable, we can express the term multiplied by $x_{ij}^{T}$
(and $x_{ij}^{P}$) as

$\displaystyle f_{ij}(z,w)$
$\displaystyle=2mA_{ij}+2A_{ij}\sum_{(s,t)\in\bar{E}}z_{st}+2mz_{ij}+2\sum_{(s,t)\in\bar{E}}w_{ijst}$
$\displaystyle-\left(d_{i}d_{j}+d_{i}\sum_{h\in V}z_{jh}+d_{j}\sum_{l\in
V}z_{il}+\sum_{h,l\in V}w_{iljh}\right)\text{.}$

We can now express the linearized constraint as:

$\displaystyle\sum_{i,j\in V}(x_{ij}^{P}-x_{ij}^{T})f_{ij}(z,w)\geq 0$ (12)

Thus, our integer programming formulation is:

$\displaystyle\min_{z\in\\{0,1\\}^{K_{|V|}\setminus E}}~{}~{}~{}$
$\displaystyle\sum_{(i,j)\in K_{|V|}\setminus E}z_{ij}$ s.t.
$\displaystyle\sum_{i,j\in V}(x_{ij}^{P}-x_{ij}^{T})f_{ij}(z,w)\geq 0\text{
for all }P\in\mathcal{P}$ (13) $\displaystyle w_{ijst}\leq z_{ij}\text{ for
all }(i,j),(s,t)$ $\displaystyle w_{ijst}\leq z_{st}\text{ for all
}(i,j),(s,t)$ $\displaystyle w_{ijst}\geq 0\text{ for all }(i,j),(s,t)$
$\displaystyle w_{ijst}\geq z_{ij}+z_{st}-1\text{ for all }(i,j),(s,t)$

######  Theorem 1.

Model (3.3) is equivalent to model (3.2).

In the framework of Algorithm 1, we would solve model (3.3) to determine
$E^{k}$, then solve model (3.1) to determine the partition that maximizes
modularity in the network $G^{k}=(V,E\cup E^{k})$. We additionally add a
constraint regarding the previous iteration’s objective value to model (3.3).
Let $c^{k}$ be the objective value in the $k^{th}$ iteration. We add the
constraint $\sum_{(i,j)\in K_{|V|}\setminus E}z_{ij}\geq c^{k}$, with $c^{0}$
initialized as $0$. By enforcing that the objective value is non-decreasing,
we can reduce the time of the branch-and-bound procedure needed to solve
iterations where the objective value stays the same as in the previous
iteration.

We note that, when extending these manipulations to the weighted version of
the problem, we use the generalized McCormick envelope [12]. For integer
variables $z_{ij},z_{st}$ with $l_{ij}\leq z_{ij}\leq u_{ij}$ and $l_{st}\leq
z_{st}\leq u_{st}$, we can represent $w_{ijst}$ as:

$\displaystyle w_{ijst}\leq l_{st}z_{ij}+u_{ij}z_{st}-l_{st}u_{ij}$
$\displaystyle w_{ijst}\leq l_{ij}z_{st}+u_{st}z_{ij}-l_{ij}u_{st}$
$\displaystyle w_{ijst}\geq l_{st}z_{ij}+l_{ij}z_{st}-l_{ij}l_{st}$
$\displaystyle w_{ijst}\geq u_{st}z_{ij}+l_{ij}z_{st}-u_{ij}u_{st}$

However, for non-binary $z_{ij}$ and $z_{st}$, we no longer have the guarantee
that $w_{ijst}=z_{ij}z_{st}$, and these substitutions will not result in an
equivalent model. A potential solution to this is to use binary expansion,
replacing each integer variable with a set of binary variables to represent
any value that integer variable can take [13]. More formally, for integer
variable $z_{ij}$ with $z_{ij}\leq u_{ij}$, we would create binary variables
$b_{ij}^{0},\ldots,b_{ij}^{k}$, with $k=\lfloor\log_{2}(u_{ij})\rfloor$. We
would then include $z_{ij}=\sum_{l=0}^{k}2^{l}b_{ij}^{l}$, and replace
$z_{ij}$ with $\sum_{l=0}^{k}2^{l}b_{ij}^{l}$ in the model, converting the
bilinear terms with integer variables to be bilinear terms with binary
variables. While this substitution would result in an equivalent model, it
also results in a large number of binary variables that does not scale well
with problem size. Preliminary computational results show this is not a
feasible approach.

###  3.4 Disjunctive Cuts and Additional Partitions

In order to better enforce that the ground truth partition is optimal, we
augment model (3.1) by defining valid disjunctive cuts, and augment Algorithm
1 by generating additional partitions. We first consider how to enforce that a
single node $v$ is clustered correctly. For ease of notation, consider a pair
of arbitrary partitions $P=(C_{1}\cup\\{v\\},C_{2},C_{3},\ldots,C_{l})$ and
$P^{\prime}=(C_{1},C_{2}\cup\\{v\\},C_{3},\ldots,C_{l})$, where $P$ is the
desired partition. Let $e_{i}(v)$ be the set of edges between $v$ and cluster
$C_{i}$ for $i=1,2$ and $e_{-}(v)$ be set of edges between $v$ and
$\bigcup_{j=3}^{l}C_{j}$. Let $E(C_{i},C_{-})$ be the set of edges between
$C_{i}$ and $\bigcup_{j=3}^{l}C_{j}$, and let $Q_{-}$ be the contribution to
modularity by $C_{3},\ldots,C_{l}$. We want to identify conditions on when
$Q_{P}<Q_{P^{\prime}}$ to derive valid inequalities that can be added to the
model (3.3) to prevent it from occurring.

$\displaystyle(4m^{2})\left(Q_{P}-Q_{P^{\prime}}\right)$
$\displaystyle=4m(|E(C_{1})|+|e_{1}(v)|)-\left(2|E(C_{1})|+2|e_{1}(v)|+|E(C_{1},C_{2})|+|e_{2}(v)|+|E(C_{1},C_{-})|+|e_{-}(v)|\right)^{2}$
$\displaystyle+4m|E(C_{2})|-\left(|2E(C_{2})|+|E(C_{1},C_{2})|+|e_{2}(v)|+|E(C_{2},C_{-})|\right)^{2}+4mQ_{-}$
$\displaystyle-4m|E(C_{1})|+\left(2|E(C_{1})|+|E(C_{1},C_{2})|+|e_{1}(v)|+|E(C_{1},C_{-})|\right)^{2}-4mQ_{-}$
$\displaystyle-4m(|E(C_{2})|+|e_{2}(v)|)+\left(2E(C_{2})+2|e_{2}(v)|+|E(C_{1},C_{2})|+|e_{1}(v)|+|E(C_{2},C_{-})|+|e_{-}(v)|\right)^{2}$
$\displaystyle=4m(|e_{1}(v)|-|e_{2}(v)|)+\left(2|E(C_{1})|+|e_{1}(v)|+|E(C_{1},C_{2})|+|E(C_{1},C_{-})|\right)^{2}$
$\displaystyle-[(2|E(C_{1})|+|e_{1}(v)|+|E(C_{1},C_{2})|+|E(C_{1},C_{-})|)+(|e_{1}(v)|+|e_{2}(v)|+|e_{-}(v)|)]^{2}$
$\displaystyle-\left(|2E(C_{2})|+|e_{2}(v)|+|E(C_{1},C_{2})|+|E(C_{2},C_{-})|\right)^{2}$
$\displaystyle+[(2|E(C_{2})|+|e_{2}(v)|+|E(C_{1},C_{2})|+|E(C_{2},C_{-})|)+(|e_{1}(v)|+|e_{2}(v)|+|e_{-}(v)|)]^{2}$
$\displaystyle=4m(|e_{1}(v)|-|e_{2}(v)|)-(|e_{1}(v)|+|e_{2}(v)|+|e_{-}(v)|)^{2}$
$\displaystyle-2\left(2|E(C_{1})|+|e_{1}(v)|+|E(C_{1},C_{2})|+|E(C_{1},C_{-}|\right)\left(|e_{1}(v)|+|e_{2}(v)|+|e_{-}(v)|\right)$
$\displaystyle+2\left(2|E(C_{2})|+|e_{2}(v)|+|E(C_{1},C_{2})|+|E(C_{2},C_{-}|\right)\left(|e_{1}(v)|+|e_{2}(v)|+|e_{-}(v)|\right)$
$\displaystyle+(|e_{1}(v)|+|e_{2}(v)|+|e_{-}(v)|)^{2}$
$\displaystyle=(4m-2|e_{1}(v)|-2|e_{2}(v)|-2|e_{-}(v)|)(|e_{1}(v)|-|e_{2}(v)|)$
$\displaystyle-2(2|E(C_{1})|+|E(C_{1},C_{-})|-2|E(C_{2})|-|E(C_{2},C_{-})|)(|e_{1}(v)|+|e_{2}(v)|+|e_{-}(v)|)$

Note that $0\leq|e_{1}(v)|+|e_{2}(v)|+|e_{-}(v)|\leq m$, so if
$|e_{1}(v)|\leq|e_{2}(v)|$ and
$2|E(C_{1})|+|E(C_{1},C_{-})|>2|E(C_{2})|+|E(C_{2},C_{-})|$, then
$Q_{P}<Q_{P^{\prime}}$. Thus, we must have that

$|e_{1}(v)|\geq|e_{2}(v)|+1\text{ or }2|E(C_{2})|+|(E(C_{2},C_{-})|\geq
2|E(C_{1})|+|(E(C_{1},C_{-})|\text{.}$

Likewise, if $|e_{1}(v)|<|e_{2}(v)|$ and $2|E(C_{1})|+|E(C_{1},C_{-})|\geq
2|E(C_{2})|+|E(C_{2},C_{-})|$, then $Q_{P}<Q_{P^{\prime}}$. This results in
needing to enforce

$|e_{1}(v)|\geq|e_{2}(v)|\text{ or }2|E(C_{2})|+|(E(C_{2},C_{-})|\geq
2|E(C_{1})|+|(E(C_{1},C_{-})|+1\text{.}$

To define these cuts regarding $T$, we must first discuss how to determine
which nodes are “classified incorrectly,” especially when our ground truth
community structure and the current optimal partition have a different number
of clusters. Let $\bar{P}$ be the partition that maximizes modularity, as
determined by Algorithm 1. We construct a complete bipartite network, where
one set of nodes represents the clusters in $T$, and the other set represents
the clusters in $\bar{P}$. For a pair of clusters $C_{i}^{T}$ and
$C_{j}^{\bar{P}}$, we set the edge weight for $(i,j)$ to be $|C^{T}_{i}\cap
C_{j}^{\bar{P}}|$, the number of nodes shared by those clusters. We then solve
the maximum weight matching problem on this bipartite network, indicating how
the clusters should be ordered. From this, we can easily determine which nodes
are incorrectly classified. Note that if the optimal partition has more
clusters than the ground truth community structure, then at least one cluster
in the optimal partition will not be matched with any cluster in the ground
truth. In this scenario, all nodes in this cluster are considered incorrectly
classified. Similarly, if the ground truth community structure has more
clusters than the optimal partition, at least one cluster in the ground truth
communities will not be matched with any cluster in the optimal partition. All
nodes in this unmatched cluster will also be considered incorrectly
classified. An example of this distinction being important is when two
clusters in the ground truth community structure are merged into a single
cluster in the optimal partition. While these nodes will all be clustered with
all of the nodes they should be clustered with, whichever ground truth cluster
is smaller will have its nodes considered incorrectly classified. This
indicates that edges (or weights) must be added to distinguish this smaller
cluster as its own individual cluster.

Once we have determined the set of incorrectly classified nodes $M$, we can
now define the disjunctive cuts. Let $v\in M$ be an incorrectly classified
node, and, without loss of generality, assume $v\in C_{1}^{T}$. Let $\bar{C}$
be the cluster $v$ belongs to in $\bar{P}$. Suppose $\bar{C}\setminus
C_{1}^{T}$ is nonempty, and consider $w\in\bar{C}\setminus C_{1}^{T}$. Without
loss of generality, assume $w\in C_{2}^{T}$. We define the partition
$P^{\prime}=(C_{1}^{T}\setminus\\{v\\},C_{2}^{T}\cup\\{v\\},C_{3}^{T},\ldots,C_{l}^{T})$.
By utilizing $\bar{P}$, we ensure that the partition generated provides a
meaningful contribution by clustering $v$ with a node with which it was
previously clustered. Besides defining the disjunctive cuts, we also include
the partition $P^{\prime}$ in $\mathcal{P}^{k}$. Computationally, this
inclusion improved solve time.

In the case when $\bar{C}\setminus C_{1}^{T}=\emptyset$, the above partition
is no longer properly defined, as $\bar{C}\subseteq C_{1}^{T}$. One option
would be to have $v$ in a cluster on its own in $P^{\prime}$. This partition,
however, is known to not be optimal as long as $d_{v}>0$ [4]. Instead, we
consider a partition that separates all of $\bar{C}$, not just a single node,
since $\bar{C}\subset C_{1}^{T}$. We can define a new partition
$P^{\prime}=(C_{1}^{T}\setminus\bar{C},\bar{C},C_{2}^{T},\ldots,C_{l}^{T})$.
This partition is well defined, since if $C_{1}^{T}\setminus\bar{C}=\emptyset$
and $\bar{C}\setminus C_{1}^{T}=\emptyset$, then $C_{1}^{T}=\bar{C}$ and thus
$v\notin M$. By including the constraint $Q_{T}\geq Q_{P^{\prime}}$, we
enforce that edges are added to merge $C_{1}^{T}\setminus\bar{C}$ and
$\bar{C}$.

At each iteration of Algorithm 1, we generate the set of incorrectly nodes
$M$, and remove nodes whose cluster in $\bar{P}$ is a subset of a cluster in
$T$. For each of these nodes $v$, we pick an arbitrary node $w$, where $v$ and
$w$ are in the same cluster in $\bar{P}$, but different clusters in $T$, to
define the disjunctive cuts and additional partition. Then, for every cluster
$C_{i}^{\bar{P}}\in\bar{P}$, we check if it is a subset of a cluster in
$C_{j}^{T}\in T$, and generate additional partitions for it when
$C_{i}^{\bar{P}}\subset C_{j}^{T}$.

##  4 Heuristic Approaches

The potential downsides of using model (3.3) include (i) that there could be a
large number of binary variables, particularly if the network is sparse, and
(ii) we must repeatedly solve model (3.3) as we generate new partitions to be
added to it. With this, even small networks can be computationally expensive
to solve, especially when trying to solve model (2.2). Thus, good heuristics
are important for obtaining high quality solutions in a reasonable amount of
time. Before we present our heuristic algorithms, we discuss some wonky
behaviours of modularity as a measure of community strength.

###  4.1 Counter-Intuitive Behavior of Modularity

As modularity is a measure of the strength of communities, there are two nice
properties we would hope to have. Intuitively, we would expect that adding an
edge between two nodes in the same cluster would increase the modularity score
of that partition. Likewise, we would expect that adding an edge between two
nodes in different clusters would decrease the modularity score of that
partition. Neither of these intuitive properties are true. Figure 1
demonstrates that adding an edge within a cluster can decrease modularity.

$0.6753$ Modularity

Add red edge: $0.6724$ Modularity

Figure 1: Example of adding an edge within a cluster decreasing modularity of
the partition

We construct this example by taking Zachary’s karate club network [14] and
removing all edges between distinct communities in the optimal partition of
the network. The modularity of the partition where each disjoint component is
its own community is $0.6753$. By adding the red edge, which is between two
nodes within the red cluster, the modularity decreases to $0.6724$.

We can likewise create an example where adding an edge across clusters
increase the modularity. Figure 2 demonstrates this behavior.

$0.1424$ Modularity

Add red edge: $0.1531$ Modularity

Figure 2: Example of adding an edge across clusters increasing modularity of
the partition

We create the example by starting with a complete network on $10$ nodes. We
then add six more nodes, and create two disjoint paths of edge-length two on
these nodes. We take one endpoint of each path, and connect them to one node
each in the complete sub-network. With the partition where the complete sub-
network and two paths are in their own clusters, we have a modularity of
$0.1424$. By adding the red edge, which connects the degree $1$ nodes in the
blue and green clusters, the modularity increases to $0.1531$.

A property of modularity being a global measure is that only the numbers of
edges within clusters and between clusters are relevant to computing
modularity, but not what specific nodes each edge connects [3]. We would then
expect that the change in modularity of a particular partition from adding an
edge is dependent only on the clusters the incident nodes are in, not the
nodes themselves. Exact computations for the change in modularity are
presented in Appendix A.

While this is a nice property for a global measure, this property makes
developing heuristics based on local properties difficult; we cannot
distinguish between edges to add to improve the modularity of our ground truth
communities. Further, the change in modularity _across different partitions_
will be different and thus understanding how adding an edge changes the
current best partition remains challenging without enumerating all partitions.
However, we still use our intuition in developing our heuristics.

We also wish to understand what edges may be necessary in creating a distinct
cluster. From Brandes et. al. [4], we know the following two properties.

######  Lemma 1.

A partition $P$ that maximizes modularity will not have disconnected clusters.

######  Lemma 2.

A partition $P$ that maximizes modularity will not have any clusters that
consist of a single node with degree $1$.

We can combine these two results into a more powerful result.

######  Theorem 2.

Given a partition $P=\\{C_{1},\ldots,C_{k}\\}$, a sparsest network to have $P$
as a partition that maximizes modularity has each cluster $C_{i}$ consisting
of a connected component whose edges form a star network.

This theorem helps us understand what edges would be useful to add to the
network to get the best gain for the smallest number of edges. In particular,
we will seek to ‘reinforce’ star-like sub-networks for each distinct cluster
in the ground truth. This also provides a lower bound on the objective value
of model (2.2).

###  4.2 Our Algorithms

Since we know exactly how adding edges to the network will impact modularity,
we can easily describe the ideal edges to add: ones that are within a single
cluster in ground truth partition $T$, and across two clusters in the
currently optimal partition $P^{*}$. Our procedure follows a simple scheme
adding edges to the network. First, we look for a “most incorrectly
classified” node to be one endpoint of the edge, which we will call $u$. We
then consider all nodes in the cluster in the ground truth partition $T$
containing $u$ that are not adjacent to $u$. We then choose the node $v$ that
minimizes the local increase in modularity. From the computations on change in
modularity in Appendix A, we expect that choosing such an edge will decrease
the gap between $Q_{P^{*}}$ and $Q_{T}$, which works towards enforcing that
$T$ is the partition that maximizes modularity.

First, we must choose which incorrectly classified node $u$ will be incident
to the added edge. We do so by comparing the number of nodes that $u$ should
be clustered with in $T$ but is currently not in the partition $P^{*}$ that
currently maximizes modularity, to the number of nodes it should be clustered
with. Let $M$ be the set of incorrectly classified nodes, where $C^{P}(u)$ is
the set of nodes in the cluster of $P$ containing $u$. We order the nodes in
$M$ in decreasing order, determined by the fraction

$q_{u}=\frac{|C^{T}(u)\setminus C^{P^{*}}(u)|}{|C^{T}(u)-1|}\text{.}$ (14)

Unless $C^{T}(u)=\\{u\\}$, the denominator of this term will be strictly
greater than $0$, and thus the term is well defined. Brandes et. al. [4]
proves that, as long as there are no isolated nodes in the network, the
partition that maximizes modularity will never have a cluster that is just a
single node. Additionally, they show that an isolated node can be placed in
any cluster (including its own) without impacting modularity. Thus, if a
ground truth partition ever has a cluster that is a single node $u$, that
problem will either be infeasible if $d_{u}>0$, or that node can be placed in
its own cluster in the partition that maximizes modularity, so it will never
be incorrectly classified.

By looking at the proportion of two terms of $q_{u}$, we can differentiate
between different size clusters. If there are two nodes $i$ and $j$ such that
$|C^{T}(i)\setminus C^{P^{*}}(i)|=|C^{T}(j)\setminus C^{P^{*}}(j)|$, and
$|C^{T}(i)|<|C^{T}(j)|$, then $\frac{|C^{T}(i)\setminus
C^{P^{*}}(i)|}{|C^{T}(i)-1|}>\frac{|C^{T}(j)\setminus
C^{P^{*}}(j)|}{|C^{T}(j)-1|}$. By normalizing the clusters by their size, we
allow for smaller clusters to more easily distinguish themselves and be
selected earlier in the heuristic.

Once we have chosen $u$, we must then choose the other endpoint node $v$. We
consider candidate nodes from the set $N(u)=\\{j\in V:j\in C^{T}(u)\setminus
C^{P^{*}}(u)\land\ A_{uv}=0\\}$. This set considers all non-adjacent nodes
that $u$ should be clustered with but currently is not. If this set is empty,
we consider a different node instead of $u$. From this set of nodes, we look
at the change in modularity attributed to $v$. Instead of summing over all
pairs of nodes in the cluster $C$, which results in a constant value
regardless of choice of $u,v$, we just sum over pairs of nodes where $v$ is at
least one of those nodes, which we define to be $Q_{v,T}(A)$. Multiplying by
$(2m)^{2}(2m+2)^{2}$, we compute this as follows:

$\displaystyle Q_{v,T}(\\{e\\})-Q_{v,T}(\emptyset)$
$\displaystyle=(2m)^{2}(2m+2)A_{vv}-(2m)^{2}(d_{v}+1)^{2}-(2m)(2m+2)^{2}A_{vv}+(2m+2)^{2}d_{v}^{2}$
$\displaystyle+2\left((2m)^{2}(2m+2)A_{uv}+(2m)^{2}(2m+2)-(2m)^{2}(d_{u}+1)(d_{v}+1)\right)$
$\displaystyle-2\left((2m)(2m+2)^{2}A_{uv}+(2m+2)^{2}d_{u}d_{v}\right)$
$\displaystyle+2\left(\sum_{u,v\neq y\in
C}\left((2m)^{2}(2m+2)A_{vy}-(2m)^{2}(d_{v}+1)d_{y}-(2m)(2m+2)^{2}A_{vy}+(2m+2)^{2}d_{v}d_{y}\right)\right)$
$\displaystyle=(-8m^{2}-8m)A_{vv}-8m^{2}d_{v}-4m^{2}+8md^{2}_{v}+4d^{2}_{v}$
$\displaystyle+2\left((-8m^{2}-8m)A_{uv}+8m^{3}+8m^{2}-(4m^{2})(d_{u}d_{v}+d_{u}+d_{v}+1)+(4m^{2}+8m+4)d_{u}d_{v}\right)$
$\displaystyle+2\left(\sum_{u,v\neq y\in
C}\left((-8m^{2}-8m)A_{vy}-4m^{2}(d_{v}d_{y}+d_{y})+(4m^{2}+8m+4)d_{v}d_{y}\right)\right)$
$\displaystyle+(-8m^{2}-8m)A_{vv}+(8m+4)d^{2}_{v}-8m^{2}d_{v}-4m^{2}$
$\displaystyle+2\left((-8m^{2}-8m)A_{uv}+(8m+4)d_{u}d_{v}+8m^{3}+4m^{2}-4m^{2}d_{u}-4m^{2}d_{v}\right)$
$\displaystyle+2\left(\sum_{u,v\neq y\in
C}\left((-8m^{2}-8m)A_{vy}+(8m+4)d_{v}d_{y}-4m^{2}dy\right)\right)$
$\displaystyle=16m^{3}+4m^{2}+(-8m^{2}-8m)A_{vv}-16m^{2}d_{v}$
$\displaystyle+2\sum_{y\neq v,y\in
C}\left((-8m^{2}-8m)A_{vy}+(8m+4)d_{v}d_{y}-4m^{2}d_{y}\right)$

We drop the $16m^{3}+4m^{2}$ term, since that will be constant across all
choices of $v$, reducing our rule to

$\delta_{v}=(-8m^{2}-8m)A_{vv}-16m^{2}d_{v}+2\sum_{y\neq v,y\in
C}\left((-8m^{2}-8m)A_{vy}+(8m+4)d_{v}d_{y}-4m^{2}d_{y}\right)\text{.}$

We choose to take the node $v^{*}=\text{arg}\min_{v\in N(u)}\delta_{v}$, as a
smaller local change in modularity indicates that node $v$ is better connected
within its cluster.

We note that our choice is closely related to the degree of each candidate
node. Intuitively, this aligns with Theorem 2. By adding edges to the highest
degree node within a cluster, we work towards creating a star sub-network in
the cluster. For incorrectly classified nodes connected to many nodes in other
clusters, we need to add more edges than just establishing a base framework.

If, for every $u\in M$, $N(u)=\emptyset$, no edge can be added such that it is
within a cluster in $T$ and across clusters in $P^{*}$. This may indicate that
some cluster of $T$ is a subset of a cluster of $P^{*}$. In this case, we may
be able to distinguish that cluster by adding edges to it. Instead of
considering $N(u)$, we consider $\tilde{N}(u)=\\{j\in V:j\in C^{T}(u)\land\
A_{uv}=0\\}$. This allows us to add edges that are both within $C^{T}(u)$ and
$C^{P^{*}}$, relaxing our restriction to include edges we expect to also be
beneficial to the modularity of $P^{*}$. If, for every $u\in M$,
$\tilde{N}(u)=\emptyset$, then every incorrectly classified node is adjacent
to every node it should be clustered with. In this case, we cannot add any
edges to the network that are incident to an incorrectly classified node,
providing a termination criterion under which the problem may be infeasible.

As with defining the disjunctive cuts, we need to understand which nodes are
incorrectly classified. We again solve the maximum weight matching problem on
our auxiliary bipartite graph to determine the correspondence between clusters
in $T$ and clusters in $P^{*}$. We formalize the full procedure in Algorithm
2.

Algorithm 2 Heuristic for Edge Addition

Initialize: initial network $G^{1}=(V,E)$, set of edges added
$E^{\prime}=\emptyset$ ground truth partition $T$, iteration counter $k=1$.

while not converged do

Step 1. Find partition $P^{k}$ that maximizes modularity of $G^{k}$. Identify
set $M$ of incorrectly classified nodes.

if $M=\emptyset$ then

Terminate, $E^{\prime}$ is a solution.

end if

Step 2. Match communities in the ground truth partition to communities in the
optimal partition.

Step 3. Order $M=\\{u_{1},\ldots u_{|M|}\\}$ in descending order based on
$\frac{|C^{T}(u)\setminus C^{P^{k}}(u)|}{|C^{T}(u)-1|}$. Set $chosen=false$,
$l=1$.

while $chosen=false$ do

if $l\leq|M|$ then

Step 3a. Identify candidate nodes $N(u)=\\{j\in V:C^{T}(u_{l})\setminus
C^{P^{k}}(u_{l})\land A_{ju_{l}}=0\\}$.

if $N(u)=\emptyset$ then

Set $l=l+1$, return to Step 3a.

else

Step 3b. Set $u=u_{l}$

for $i\in N(u)$ do

Step 3c. Choose $v=\text{arg}\min_{i\in N(u)}\delta_{i}$.

end for

end if

else if $l>|M|$ and $l\leq 2|M|$ then

Step 3a. Identify candidate nodes $\tilde{N}(u)=\\{j\in V:C^{T}(u_{l})\land
A_{ju_{l}}=0\\}$.

if $\tilde{N}(u)=\emptyset$ then

Set $l=l+1$, return to Step 3a.

else

Step 3b. Set $u=u_{l}$

for $i\in N(u)$ do

Step 3c. Choose $v=\text{arg}\min_{i\in\tilde{N}(u)}\delta_{i}$.

end for

end if

else

Terminate, problem is infeasible.

end if

end while

Step 4. Set $G^{k+1}=(V(G^{k}),E(G^{k})\cup(u,v))$,
$E^{\prime}=E^{\prime}\cup(u,v)$, $k=k+1$. Return to Step 1.

end while

Extending this heuristic to weighted networks is quite simple. While
computations of $m$ and $d_{i}$ change, the above computation remains the
same. We additionally allow for weight to be added to edges already existing
in the network, removing the condition $A_{uv}=0$ from the definition of
$N(u)$ and $\tilde{N}(u)$. In a weighted network, we note that we could have
instead added weight to an edge until the partition changes. Even if we were
to do so, this does not guarantee that an incorrectly classified node will be
correctly classified in future iterations nor that an edge will not be re-
visited by the heuristic in future iterations, as demonstrated in the example
in Figure 3. In this example, each edge has weight $1$, and the weight limit
for each edge is $10$.

(a) Optimal Partition

(b) Ground Truth

(c) Ground Truth with Added Weights

Figure 3: Example network where an edge needs additional weight added after
first visit by the heuristic approach

In this network, weight is added to the colored edges in Figure 3(c). In the
first iteration, a unit of weight is added to the red edge, $(1,5)$, causing
nodes $1$ and $5$ to be clustered together. The second iteration adds the blue
edge, $(2,4)$. The third and fourth iterations add a unit of weight to the
green edge $(5,7)$. The addition of this weight causes nodes $5$ and $7$ to be
clustered together, but separates node $1$ from node $5$. Our algorithm next
adds another weight to the red edge. At the beginning, only one unit weight is
needed for node $1$ to be classified correctly. However, after more
iterations, node $1$ is no longer classified correctly, and an edge needs to
be revisited.

When considering model (2.2), instead of selecting which edge to include from
the edge set, we turn our attention to removing edges from the network. The
idea behind this procedure is quite intuitive as well; order the edges in the
network, then, one by one, check if removing the $i^{th}$ edge results in a
change in the optimal partition. If the optimal partition does not change,
then that edge can be removed from the network. Note that, if the $i^{th}$
edge was removed, that edge could have been critical for changing the optimal
partition when removing any of the first $i-1$ edges. Thus, if any single edge
is removed, we must check the full edge list again to determine if more edges
can be removed. We repeat this procedure until it determines that removing any
single edge in the network results in a change of the optimal partition. We
formalize this approach in Algorithm 3.

Algorithm 3 Heuristic for Edge Removal

Initialize: initial network $G^{1}=(V,E)$, set of edges removed
$E^{\prime}=\emptyset$, optimal partition $T$, $check=true$, iteration counter
$k=1$.

while $check=true$ do

Step 1. Set $check=false$, $I=\emptyset$.

Step 2. Order edges $E(G^{k})=\\{e_{1},\ldots,e_{|E(G^{k})|}\\}$.

for $i\in 1:|E(G^{k})|$ do

Step 2. Create network
$G^{\prime}=(V(G^{k}),E(G^{K})\setminus(I\cup\\{e_{i}\\}))$.

Step 3. Find the optimal partition $P^{*}$ of $G^{\prime}$.

if $P^{*}=T$ then

Step 3a. Set $I=I\cup\\{e_{i}\\}$, $check=true$.

end if

end for

Step 4. Create network $G^{k+1}=(V(G^{k}),E(G^{k})\setminus I)$.

Step 5. Set $E^{\prime}=E^{\prime}\cup I$.

end while

We test the following rules for ordering the edges:

  1. 1.

a prespecified order,

  2. 2.

a random permutation,

  3. 3.

a fixed order based on contribution to modularity,

  4. 4.

a dynamically updated order based on contribution to modularity.

A weakness of the first rule is that the prespecified order may choose edges
whose removal may cause fewer edges to be removed in subsequent iterations.
The second rule randomizes the order in which the edges are processed in each
iteration. This provides a workaround for the weakness of the first rule, but
the quality of the solution obtained using this rule is determined by the
orders randomly generated.

The third and fourth rules use an edge’s contribution to modularity to
determine order. We order the edges in decreasing order of their values in the
modularity matrix, $M_{ij}=A_{ij}-\frac{d_{i}d_{j}}{2m}$. Intuitively, edges
with a large value of $M_{ij}$ have small products of degrees, meaning this
edge has a low probability to be generated in the corresponding random
network. We suspect that these edges are unlikely to be impactful in
determining the optimal partition, making them good candidates for removal.
For our third rule, we choose to pre-compute this ordering and use it in each
iteration, as with the first rule. For the fourth rule, we update this
ordering after each iteration. Since it is possible that the removal of
certain edges can cause changes of importance of other edges to the optimal
partition, the updates may result in more edges removed.

Since Algorithm 3 ensures that each subsequent network has the same optimal
partition as the original network, we can improve the solution found by
Algorithm 2 by taking the output solution $G^{\prime}$ as the input for
Algorithm 3. This removes edges that were added, but later become redundant
due to the addition of other edges. For this addition as post-processing, we
simply choose to process the edges in the order in which they were added.
Extending this algorithm to weighted networks is also straightforward. Instead
of removing one edge at a time, we remove a unit of weight on one edge at a
time.

##  5 Computational Results

We test our methods on multiple networks, some of which are popular in
community detection literature and some popular in illicit network analysis.
We retrieved this data from the CASOS Public Datasets [15] and Pajek Datasets
[16]. For the Montreal Gangs data set, which has disconnected components, we
delete the nodes that are not in the largest connected component. We model the
integer programming problem in AMPL with CPLEX 12.10 [17] as the solver.
Experiments are conducted on a laptop with an Intel® CoreTM i5-8250 CPU @ 1.6
GHz - 1.8 GHz and 16 GB RAM running Windows 10. We set a time limit of $7200$
seconds. We develop the heuristics in R using the iGraph package [18]. Due to
integer programming limitations, we additionally restrict the set of edges
considered for addition in thee experiments. To ensure feasibility, we start
with the set of edges identified by the heuristics with post-processing and
without edge restrictions as the initial set of allowable edges. We then
generate the list of all pairs of nodes with shortest path length $2$ between
them, and randomly select from that list to add to the set of allowable edges.
We consider $100$ edges for each network.

###  5.1 Results on Unweighted Networks

All networks except the Ciel data set have a known community structure. For
networks with a known community structure, we use it as the ground truth. For
the Ciel data set, we use the Louvain method [19] to find a near-optimal
partition, and use it as the ground truth. For networks that are originally
weighted, we remove edge weights and experiment on their unweighted versions.
Table 1 reports the networks we use for our experiments organized by the
number of edges in the network. Figures 4 and 5 demonstrate the differences in
the optimal partition and the ground truth community structure for example
networks.

Network | $|V|$ | $|E|$
---|---|---
Ciel [20] | 25 | 35
Caviar6 [20] | 27 | 47
Sawmill [21] | 36 | 62
Rhodes [22] | 22 | 66
Dolphins small [23] | 40 | 70
Montreal Gangs [24] | 29 | 75
Karate [14] | 34 | 78
Table 1: Size of test networks

(a) Optimal Partition

(b) Ground Truth

Figure 4: Optimal partition and ground truth community structure of Rhodes
Network

(a) Optimal Partition

(b) Ground Truth

Figure 5: Optimal partition and ground truth community structure of Karate
Club Network

Table 2 reports the initial number of incorrectly classified nodes as well as
the number of edges added to each network to reach the desired communities.
Table 3 reports the run time of each method. We denote the integer program
without disjunctive cuts and extra partitions as IP, and the integer program
with disjunctive cuts and additional partitions as IP+. Entries marked with *
did not converge to the optimal solution within the time limit, and the lower
bound on the optimal solution is reported.

Network | Init # Misclassified | #Edges Added (IP) | #Edges Added (IP+) | #Edges Added (Alg2) | #Edges Added(Alg2&3)
---|---|---|---|---|---
Ciel | 3 | 1 | 1 | 1 | 1
Caviar6 | 3 | 11 | 11 | 12 | 12
Sawmill | 10 | 9 | 9 | 16 | 13
Rhodes | 10 | 18* | 18* | 34 | 34
Dolphins small | 17 | 24* | 29* | 34 | 30
Montreal Gangs | 16 | 18* | 18* | 45 | 42
Karate | 12 | 8 | 8 | 11 | 9

Table 2: Comparison of Number of Edges Added to Tested Networks with Edge Restrictions Network | IP (nParts) | IP+ (nParts) | Alg2 | Alg2&3
---|---|---|---|---
Ciel | 1.58 (1) | 1.83 (1) | 0.09 | 0.13
Caviar6 | 39.92 (8) | 39.13 (5) | 4.95 | 10.68
Sawmill | 287.08 (14) | 2968.05 (8) | 4.83 | 14.25
Rhodes | 7200* (9) | 7200* (6) | 4.69 | 7.31
Dolphins small | 7200* (8) | 7200* (9) | 26.47 | 258.66
Montreal Gangs | 7200* (2) | 7200* (2) | 627.23 | 1256.49
Karate | 1725.09 (16) | 488.25 (12) | 3.80 | 13.82
Table 3: Run Time (in seconds) on Tested Networks with Edge Restrictions

Table 2 shows that the integer programming formulation is able to find a
relatively small number of edges as compared to the density of the network in
the case of the Karate club and Ciel networks. For some networks, such as
Caviar6 and Montreal Street Gangs, the layout of the network requires edges to
be added to distinguish a smaller cluster from a larger cluster that can
relatively easily absorb it, as well as a significant number of edges that
connect nodes that are clustered together in the ground truth, but do not have
any edges in the between them in the network. However, for the other networks,
the integer programming formulation is not able to solve the problem within
the time limit. This is due to the edge-anonymity property of modularity, as
that unintentionally creates a lot of symmetry in the branch-and-bound
procedure for the master problem. Integer programming problems with symmetry
are particularly difficult to solve [25]. Additionally, as the number of edges
being included increases, the number of symmetric solutions increases
exponentially. We observed that, in early iterations and in subsequent
iterations where the master problem objective value stays the same, the master
problem is solved to optimality relatively quickly. However, in later
iterations when the number of edges needed must be larger than in the previous
iteration, the master problem may not necessarily solve to optimality due to
the number of branches that need to be explored. We suspect that, in each
iteration of the master problem, the integer program identifies the optimal
solution with respect to the set of previously explored partitions, but does
not complete the branching procedure to verify that the identified solution is
indeed optimal.

Including disjunctive cuts and additional partitions improves the solve time
of the integer programming formulation for the Karate club network, and only
hinders the solve time of the Sawmill network. This indicates that the extra
information generated may have prevented a useful partition that was visited
by the base integer program from being visited by the augmented procedure. For
the Dolphins network, the best known lower bound increases from $24$ to $29$,
drastically improving the result, and improving the known quality of the
heuristic solution. For Montreal Gangs and Rhodes, no improvements are made by
including the disjunctive cuts and additional partitions. This indicates that
the edge-anonymity caused the solver to stall early on in this procedure.
Increasing the run time to $36000$ seconds on these instances shows that
little-to-no progress is made in improving these lower bounds for either
procedure. This demonstrates that problem difficulty is not only dependent on
the size of the network, but the density of the network and difference between
the ground truth community structure and initial optimal partition.

The heuristic approach is able to find reasonably high quality solutions as
compared to the integer program optimal solutions before post-processing is
implemented. With post processing, we are able to reduce the number of edges
needed in some networks. Furthermore, where the integer program may not solve
the problem in the time limit, the heuristic approach is able to find a
solution in under an hour in all instances. When the integer program does
solve the problem within the time limit, the heuristic approach is able to
solve the problem significantly faster. We note that, since some of the
allowable edges are determined by the heuristic without edge restrictions,
adding the same edge restrictions to the heuristic will result in the same
solution. However, allowing any edge currently not in the network to be added
to the integer programming formulation would drastically increase the run
time.

We observe that when the ground truth communities are relatively similar to
the initial optimal community structure, fewer edges are needed. Additionally,
if the ground truth community structure has clear boundaries between the
communities, fewer edges are needed. Figures 4 and 5 help distinguish this
difference. In Figure 5(b), there are very few edges between the two clusters,
and there are very few nodes that have edges that cross between clusters. This
results in fewer edges being needed to properly cluster these boundary nodes.
However, in Figure 4(b), the distinction between the clusters is significantly
less clear, resulting in a need for each cluster to have more edges added to
properly distinguish them. Another example of this is the difference between
the Caviar6 and Ciel networks. While both have only three initially
incorrectly classified nodes, the Caviar6 network has a cluster that is mostly
absorbed into another cluster, where the boundaries between the clusters in
the Ciel network are more distinct. This causes a need for more edges in the
Caviar6 network to properly distinguish this cluster.

###  5.2 Results on Weighted Networks

We also test our methods on weighted versions of each network used in our
unweighted test set. Networks that are already weighted have an upper bound on
the weight of edge set to the weight of the maximum weight edge. For networks
that are unweighted, we randomly generate weights between $1$ and $15$ for
each edge, with a maximum edge weight of $15$ allowed. We use the unweighted
optimal partition as the ground truth, and verify that the optimal partition
for the weighted network is not the same as that of the unweighted network. We
note that, unlike with unweighted networks, we cannot compare the heuristic
results directly against the integer programming results. The integer program
cannot find the optimal solution due to the inaccuracy in the McCormick
inequalities. Table 4 reports the weights on networks used for experiments,
using the same ordering as above.

Network | Total Weight
---|---
Ciel | 198
Caviar6 | 668
Sawmill | 481
Rhodes | 524
Dolphins small | 544
Montreal Gangs | 621
Karate | 654
Table 4: Total Weights of Edge Sets of Test Networks

As before, we restrict the set of edges to which weight can be added to in
each network. We restrict this to be the set of currently existing edges,
enforcing that no new edges can be added. Additionally, we restrict the total
weight that an edge can have. Table 5 reports the initial number of
incorrectly classified nodes and the amount of weight added to each network.
Table 6 reports the run time of each method.

Network | Init. #Misclassified | Weight Added (Alg2) | Weight Added(Alg2&3)
---|---|---|---
Ciel | 3 | 12 | 12
Caviar6 | 7 | 115 | 113
Sawmill | 16 | 30 | 26
Rhodes | 6 | 38 | 31
Dolphins small | 8 | 61 | 58
Montreal Gangs | 11 | 45 | 43
Karate | 13 | 45 | 44

Table 5: Comparison of Amount of Weight Added to Tested Networks with Edge Restrictions Network | Alg2 | Alg2&3
---|---|---
Ciel | 0.51 | 0.96
Caviar6 | 17.42 | 60.38
Sawmill | 11.02 | 48.58
Rhodes | 2.64 | 7.37
Dolphins small | 20.79 | 68.19
Montreal Gangs | 51.06 | 169.83
Karate | 46.15 | 120.11
Table 6: Run Time (in seconds) on Tested Networks with Edge Restrictions

For all networks except Caviar6, we are able to add less than $12\%$ of
initial total weight to the network in order to return to the desired
clustering. With Caviar6, we see that this percentage is much higher, at
$17\%$ of the initial total weight. This is likely due to the large disparity
in distribution of edge weights; there are a few nodes with a lot of weight on
incident edges, but many nodes with a small amount of weight on incident
edges. Being able to correctly classify these nodes with a small amount of
weight on incident edges takes much more weight. Similarly to the unweighted
case, the post processing is able to reduce the total weight added by a small
amount, but this decrease is nowhere near as significant as in the unweighted
case.

We next report results on weighted networks without edge restrictions. Table 7
reports the initial number of incorrectly classified nodes, as well as the
number of edges added to each network to reach the desired communities. Table
8 reports the run time of each method.

Network | Init. #Misclassified | Weight Added (Alg2) | Weight Added(Alg2&3)
---|---|---|---
Ciel | 3 | 12 | 12
Caviar6 | 7 | 115 | 113
Sawmill | 16 | 26 | 22
Rhodes | 6 | 28 | 26
Dolphins small | 8 | 62 | 58
Montreal Gangs | 11 | 48 | 46
Karate | 13 | 45 | 44

Table 7: Comparison of Amount of Weight Added to Tested Networks without Edge Restrictions Network | Alg2 | Alg2&3
---|---|---
Ciel | 0.50 | 0.96
Caviar6 | 17.57 | 60.15
Sawmill | 9.67 | 36.08
Rhodes | 1.72 | 6.93
Dolphins small | 21.22 | 52.27
Montreal Gangs | 34.32 | 90.14
Karate | 50.97 | 145.72
Table 8: Run Time (in seconds) on Tested Networks without Edge Restrictions

In this case, removing edge restrictions results in similar overall weight
being added to the networks, and has run times comparable to that of the
restricted edge case. An interesting result is that, less weight needs to be
added to the Dolphins and Montreal Gangs networks when edge restrictions are
enforced. Recall, that the proposed algorithms act in a greedy fashion by
focusing on local decisions at each point. These results indicate that adding
weight to any possible edge can result in adding weight to a locally good
choice, but creates global difficulties in later iterations.

###  5.3 Results for Sparsest Representations of Networks

Here, we present the results of finding the sparsest representation of
networks that share the same optimal partition maximizing modularity. The
edge-anonymity prevented the integer programming formulation from making
meaningful progress towards a solution, even with additional constraints
enforcing that clusters are connected components. We present the number of
edges remaining and run time on each network for all four rules in Table 9,
with run time being reported in Table 10. For rule $2$, we run $100$ tests,
and report the median, minimum and maximum of the number of edges remaining,
and the mean and standard deviation of the run times. We additionally report
the number of edge remaining after removing cross cluster edges as Pre, and
the number of edges determined by Theorem 2 in the column LB.

Network | #Edges | Rule 1 | Rule 2 Med (Min, Max) | Rule 3 | Rule 4 | LB
---|---|---|---|---|---|---
Ciel | 35 | 26 | 26 (25, 28) | 25 | 25 | 22
Caviar6 | 47 | 23 | 23 (23, 35) | 23 | 23 | 23
Sawmill | 62 | 46 | 34 (33, 46) | 34 | 34 | 32
Rhodes | 66 | 21 | 23 (21, 27) | 21 | 21 | 19
Dolphins small | 70 | 53 | 43 (38, 68) | 42 | 42 | 34
Montreal Gangs | 75 | 35 | 32 (29, 36) | 35 | 36 | 26
Karate | 78 | 31 | 33 (30, 36) | 31 | 31 | 30

Table 9: Comparison of Number of Edges Remaining in Tested Networks Network | Rule 1 | Rule 2 Mean (SD) | Rule 3 | Rule 4
---|---|---|---|---
Ciel | 2.51 | 2.11 (0.46) | 1.70 | 1.69
Caviar6 | 12.44 | 9.25 (2.90) | 11.54 | 11.59
Sawmill | 14.78 | 22.09 (8.40) | 23.96 | 23.76
Rhodes | 6.00 | 2.97 (0.84) | 2.72 | 2.64
Dolphins small | 98.06 | 35.17 (8.69) | 24.67 | 23.89
Montreal Gangs | 47.93 | 111.50 (57.73) | 172.28 | 177.58
Karate | 27.11 | 26.41 (13.37) | 23.06 | 23.09
Table 10: Run Time (in seconds) on Tested Networks

Network | #Edges | Pre | Rule 1 | Rule 2 Med (Min, Max) | Rule 3 | Rule 4 | LB
---|---|---|---|---|---|---|---
Ciel | 35 | 29 | 26 | 26 (25, 26) | 25 | 25 | 22
Caviar6 | 47 | 27 | 23 | 23 (23, 23) | 23 | 23 | 23
Sawmill | 62 | 51 | 46 | 34 (33, 43) | 34 | 34 | 32
Rhodes | 66 | 46 | 23 | 23 (21, 28) | 21 | 21 | 19
Dolphins small | 70 | 63 | 42 | 42 (38, 57) | 42 | 42 | 34
Montreal Gangs | 75 | 44 | 34 | 33 (29, 37) | 33 | 33 | 26
Karate | 78 | 57 | 31 | 33 (30, 40) | 30 | 30 | 30

Table 11: Comparison of Number of Edges Remaining in Tested Networks with Pre-Processing Network | Rule 1 | Rule 2 Mean (SD) | Rule 3 | Rule 4
---|---|---|---|---
Ciel | 1.46 | 1.49 (0.13) | 1.42 | 1.45
Caviar6 | 2.18 | 2.42 (0.37) | 2.50 | 2.31
Sawmill | 12.76 | 11.05 (1.39) | 10.78 | 10.75
Rhodes | 1.76 | 1.65 (0.30) | 1.61 | 1.51
Dolphins small | 18.04 | 20.63 (3.61) | 19.95 | 19.97
Montreal Gangs | 11.75 | 12.25 (1.40) | 11.29 | 11.74
Karate | 8.78 | 8.87 (0.93) | 10.03 | 9.86
Table 12: Run Time (in seconds) on Tested Networks with Pre-Processing

Here we see that for almost every data set, the rules ordering the edges based
on contribution to modularity either perform as well as or better than the
other rules. The minimums obtained by the second rule are at least as good as
these rules, but the medians rarely outperform the other rules. Additionally,
the high standard deviations in run time make it less appealing, even when it
is the fastest rule. We note that the first rule only happens to work well
when the ordering of the edges is already satisfactory. While it performs well
here, it is not guaranteed to do so in general, such as with the Dolphins data
set, where the third and fourth rule will be consistent regardless of initial
edge ordering. We note that, on these networks, the dynamic reordering impacts
execution time negligibly, but does not improve the results.

##  6 Conclusions and Future Work

We define two new problems regarding the modularity of a network. In the first
problem, we seek to find the minimum set of edges (or weights) to add to a
network to enforce that a given partition maximizes modularity. In the second
problem, we seek to find the minimum set of edges in the network such that the
partition optimizing modularity in the reduced network and the full network is
the same. We provide an integer programming framework for these problems, and
augment the method with the generation of disjunctive cuts and additional
partitions. We demonstrate how the wonky behavior of modularity proves
problematic for the branch-and-bound procedure and prevents proving nice
properties about how adding edges impacts modularity. The integer programming
formulation is able to solve the first problem when the ground truth
communities are relatively distinct or similar to the initial optimal
clusters, demonstrating that the difficulty of the problem is not just
dependent on the size of the network. We also devise heuristics that are able
to find high quality solutions to the problem significantly quicker than the
integer programming formulation, and finds solutions to the problem on
instances that the integer programming formulation is not able to solve. For
the second problem, our heuristics find that we are able to identify a small
number of edges necessary to maintain the original community structure, often
within a small percentage of edges away from identifying star networks within
each cluster.

Future work includes extending these methods to work for other optimization
based clustering measures, such as modularity density, as well as the
development of heuristics for these measures [26]. Additional future work
includes improving the components of the integer programming framework, such
as symmetry breaking techniques in the branch-and-bound procedure. Improving
the tightness of the McCormick inequalities in the weighted problem when
enforcing a given partition is optimal as compared to previously visited
partitions is also important to being able to apply the integer program to
weighted networks. This framework can also be improved upon to utilize more
advanced means of solving for the partition that optimizes modularity,
allowing for these problems to be solved more efficiently and on larger
networks.

## Acknowledgements

This work was supported by the U.S. Department of Homeland Security under
Grant Award Number 2017-ST061-CINA01. The views and conclusions contained in
this document are those of the authors and should not be interpreted as
necessarily representing the official policies, either expressed or implied,
of the U.S. Department of Homeland Security.

## References

  * [1] K. Taha, “Static and dynamic community detection methods that optimize a specific objective function: A survey and experimental evaluation,” IEEE Access, 2020.
  * [2] A. Bahulkar, B. K. Szymanski, N. O. Baycik, and T. C. Sharkey, “Community detection with edge augmentation in criminal networks,” in 2018 IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining (ASONAM), pp. 1168–1175, 2018.
  * [3] M. E. J. Newman and M. Girvan, “Finding and evaluating community structure in networks,” Physical Review E, vol. 69, Feb 2004.
  * [4] U. Brandes, D. Delling, M. Gaertler, R. Gorke, M. Hoefer, Z. Nikoloski, and D. Wagner, “On modularity clustering,” IEEE Transactions on Knowledge and Data Engineering, vol. 20, no. 2, pp. 172–188, 2007.
  * [5] V. Martínez, F. Berzal, and J.-C. Cubero, “A survey of link prediction in complex networks,” ACM Computing Surveys (CSUR), vol. 49, no. 4, pp. 1–33, 2016.
  * [6] M. Chen, K. Kuzmin, and B. K. Szymanski, “Community detection via maximization of modularity and its variants,” IEEE Transactions on Computational Social Systems, vol. 1, no. 1, pp. 46–65, 2014.
  * [7] G. Agarwal and D. Kempe, “Modularity-maximizing graph communities via mathematical programming,” The European Physical Journal B, vol. 66, no. 3, pp. 409–418, 2008.
  * [8] D. Aloise, S. Cafieri, G. Caporossi, P. Hansen, S. Perron, and L. Liberti, “Column generation algorithms for exact modularity maximization in networks,” Physical Review E, vol. 82, no. 4, p. 046112, 2010.
  * [9] M. Fischetti, I. Ljubic, M. Monaci, and M. Sinnl, “Interdiction games and monotonicity, with application to knapsack problems,” INFORMS Journal on Computing, vol. 31, pp. 390–410, 2019.
  * [10] M. E. Newman, “Equivalence between modularity optimization and maximum likelihood methods for community detection,” Physical Review E, vol. 94, no. 5, p. 052315, 2016.
  * [11] B. Serrano and T. Vidal, “Community detection in the stochastic block model by mixed integer programming,” submitted to arXiv, 2021.
  * [12] G. P. McCormick, “Computability of global solutions to factorable nonconvex programs: Part I—Convex underestimating problems,” Mathematical Programming, vol. 10, no. 1, pp. 147–175, 1976.
  * [13] J. H. Owen and S. Mehrotra, “On the value of binary expansions for general mixed-integer linear programs,” Operations Research, vol. 50, no. 5, pp. 810–819, 2002.
  * [14] W. W. Zachary, “An information flow model for conflict and fission in small groups,” Journal of Anthropological Research, pp. 452–473, 1977.
  * [15] CASOS, “Public datasets.” http://www.casos.cs.cmu.edu/tools/datasets/external/index.php.  Accessed Dec 2020.
  * [16] V. Batagelj and A. Mrvar, “Pajek datasets.” http://vlado.fmf.uni-lj.si/pub/networks/data/esna/default.htm, 2006.  Accessed Dec 2020.
  * [17] IBM ILOG CPLEX Optimization Studio, “Cplex users manual, version 12.7,” 2017\.
  * [18] G. Csardi and T. Nepusz, “The igraph software package for complex network research,” InterJournal, vol. Complex Systems, p. 1695, 2006.
  * [19] V. D. Blondel, J.-L. Guillaume, R. Lambiotte, and E. Lefebvre, “Fast unfolding of communities in large networks,” Journal of Statistical Mechanics: Theory and Experiment, vol. 2008, no. 10, p. P10008, 2008.
  * [20] C. Morselli, Inside criminal networks.  Springer, 2009.
  * [21] J. H. Michael and J. G. Massey, “Modeling the communication network in a sawmill,” Forest Products Journal, vol. 47, no. 9, p. 25, 1997.
  * [22] C. J. Rhodes and E. Keefe, “Social network topology: a bayesian approach,” Journal of the Operational Research Society, vol. 58, no. 12, pp. 1605–1611, 2007.
  * [23] D. Lusseau, K. Schneider, O. J. Boisseau, P. Haase, E. Slooten, and S. M. Dawson, “The bottlenose dolphin community of Doubtful Sound features a large proportion of long-lasting associations,” Behavioral Ecology and Sociobiology, vol. 54, no. 4, pp. 396–405, 2003.
  * [24] K. Descormiers and C. Morselli, “Alliances, conflicts, and contradictions in Montreal’s street gang landscape,” International Criminal Justice Review, vol. 21, no. 3, pp. 297–314, 2011.
  * [25] J. Ostrowski, J. Linderoth, F. Rossi, and S. Smriglio, “Orbital branching,” Mathematical Programming, vol. 126, no. 1, pp. 147–178, 2011.
  * [26] M. Chen, T. Nguyen, and B. K. Szymanski, “On measuring the quality of a network community structure,” in 2013 International Conference on Social Computing, pp. 122–127, IEEE, 2013.

##  Appendix A Appendix A: Computation of Change in Modularity

###  A.1 Change in Modularity from Including an Edge Within a Cluster

###### Proof.

Consider adding edge $e=(u,v)$ where $u,v\in C_{1}$. We compute the difference
$(Q_{T}(\\{e\\})-Q_{T}(\emptyset)$. To simplify fractions, we multiply by
$(2m)^{2}(2m+2)^{2}$.

$\displaystyle(2m)^{2}(2m+2)^{2}(Q_{T}(\\{e\\})-Q_{T}(\emptyset)$
$\displaystyle=(2m)^{2}(2m+2)A_{uu}-(2m)^{2}(d_{u}+1)^{2}-(2m)(2m+2)^{2}A_{uu}+(2m+2)^{2}d_{u}^{2}$
$\displaystyle+(2m)^{2}(2m+2)A_{vv}-(2m)^{2}(d_{v}+1)^{2}-(2m)(2m+2)^{2}A_{vv}+(2m+2)^{2}d_{v}^{2}$
$\displaystyle+2\left((2m)^{2}(2m+2)A_{uv}+(2m)^{2}(2m+2)-(2m)^{2}(d_{u}+1)(d_{v}+1)\right)$
$\displaystyle+\left(-(2m)(2m+2)^{2}A_{uv}+(2m+2)^{2}d_{u}d_{v}\right)$
$\displaystyle+2\left(\sum_{u,v\neq x\in
C_{1}}\left((2m)^{2}(2m+2)A_{ux}-(2m)^{2}(d_{u}+1)d_{x}-(2m)(2m+2)^{2}A_{ux}+(2m+2)^{2}d_{u}d_{x}\right)\right)$
$\displaystyle+2\left(\sum_{u,v\neq y\in
C_{1}}\left((2m)^{2}(2m+2)A_{vy}-(2m)^{2}(d_{v}+1)d_{y}-(2m)(2m+2)^{2}A_{vy}+(2m+2)^{2}d_{v}d_{y}\right)\right)$
$\displaystyle+\sum_{u,v\neq x,y\in
C_{1}}\left((2m)^{2}(2m+2)A_{xy}-(2m)^{2}d_{x}d_{y}-(2m)(2m+2)^{2}A_{xy}+(2m+2)^{2}d_{x}d_{y}\right)$
$\displaystyle+\sum_{C_{i}\in T,i\neq 1}\sum_{x,y\in
C_{i}}\left((2m)^{2}(2m+2)A_{xy}-(2m)^{2}d_{x}d_{y}-(2m)(2m+2)^{2}A_{xy}+(2m+2)^{2}d_{x}d_{y}\right)$
$\displaystyle=(-8m^{2}-8m)A_{uu}-8m^{2}d_{u}-4m^{2}+8md^{2}_{u}+4d^{2}_{u}$
$\displaystyle+(-8m^{2}-8m)A_{vv}-8m^{2}d_{v}-4m^{2}+8md^{2}_{v}+4d^{2}_{v}$
$\displaystyle+2\left((-8m^{2}-8m)A_{uv}+8m^{3}+8m^{2}-(4m^{2})(d_{u}d_{v}+d_{u}+d_{v}+1)+(4m^{2}+8m+4)d_{u}d_{v}\right)$
$\displaystyle+2\left(\sum_{u,v\neq x\in
C_{1}}\left((-8m^{2}-8m)A_{ux}-4m^{2}(d_{u}d_{x}+d_{x})+(4m^{2}+8m+4)d_{u}d_{x}\right)\right)$
$\displaystyle+2\left(\sum_{u,v\neq y\in
C_{1}}\left((-8m^{2}-8m)A_{vy}-4m^{2}(d_{v}d_{y}+d_{y})+(4m^{2}+8m+4)d_{v}d_{y}\right)\right)$
$\displaystyle+\sum_{u,v\neq x,y\in
C_{1}}\left((-8m^{2}-8m)A_{xy}+(8m+4)d_{x}d_{y}\right)$
$\displaystyle+\sum_{C_{i}\in T,i\neq 1}\sum_{x,y\in
C_{i}}\left((-8m^{2}-8m)A_{xy}+(8m+4)d_{x}d_{y}\right)$
$\displaystyle=(-8m^{2}-8m)A_{uu}+(8m+4)d^{2}_{u}-8m^{2}d_{u}-4m^{2}$
$\displaystyle+(-8m^{2}-8m)A_{vv}+(8m+4)d^{2}_{v}-8m^{2}d_{v}-4m^{2}$
$\displaystyle+2\left((-8m^{2}-8m)A_{uv}+(8m+4)d_{u}d_{v}+8m^{3}+4m^{2}-4m^{2}d_{u}-4m^{2}d_{v}\right)$
$\displaystyle+2\left(\sum_{u,v\neq x\in
C_{1}}\left((-8m^{2}-8m)A_{ux}+(8m+4)d_{u}d_{x}-4m^{2}d_{x}\right)\right)$
$\displaystyle+2\left(\sum_{u,v\neq y\in
C_{1}}\left((-8m^{2}-8m)A_{vy}+(8m+4)d_{v}d_{y}-4m^{2}dy\right)\right)$
$\displaystyle+\sum_{u,v\neq x,y\in
C_{1}}\left((-8m^{2}-8m)A_{xy}+(8m+4)d_{x}d_{y}\right)$
$\displaystyle+\sum_{C_{i}\in T,i\neq 1}\sum_{x,y\in
C_{i}}\left((-8m^{2}-8m)A_{xy}+(8m+4)d_{x}d_{y}\right)$
$\displaystyle=16m^{3}-\sum_{x\in C_{1}}16m^{2}d_{x}+\sum_{C_{i}\in
T}\sum_{x,y\in C_{i}}\left((-8m^{2}-8m)A_{xy}+(8m+4)d_{x}d_{y}\right)$

Thus, the increase in modularity by adding $(u,v)$ is

$\frac{16m^{3}-\sum_{x\in C_{1}}16m^{2}d_{x}+\sum_{C_{i}\in T}\sum_{x,y\in
C_{i}}\left((-8m^{2}-8m)A_{xy}+(8m+4)d_{x}d_{y}\right)}{(2m)^{2}(2m+2)^{2}}\text{.}$

∎

###  A.2 Change in Modularity from Including an Edge Between Clusters

###### Proof.

Consider adding edge $e=(u,v)$ where $u\in C_{1}$ and $v\in C_{2}$. We compute
the difference $(Q_{T}(\\{e\\})-Q_{T}(\emptyset)$. We simplify fractions by
multiplying by $(2m)^{2}(2m+2)^{2}$.

$\displaystyle(2m)^{2}(2m+2)^{2}(Q_{T}(\\{e\\})-Q_{T}(\emptyset)$
$\displaystyle=(2m)^{2}(2m+2)A_{uu}-(2m)^{2}(d_{u}+1)^{2}-(2m)(2m+2)^{2}A_{uu}+(2m+2)^{2}d_{u}^{2}$
$\displaystyle+2\left(\sum_{u\neq x\in
C_{1}}\left((2m)^{2}(2m+2)A_{ux}-(2m)^{2}(d_{u}+1)d_{x}-(2m)(2m+2)^{2}A_{ux}+(2m+2)^{2}d_{u}d_{x}\right)\right)$
$\displaystyle+\sum_{u\neq x,y\in
C_{1}}\left((2m)^{2}(2m+2)A_{xy}-(2m)^{2}d_{x}d_{y}-(2m)(2m+2)^{2}A_{xy}+(2m+2)^{2}d_{x}d_{y}\right)$
$\displaystyle+(2m)^{2}(2m+2)A_{vv}-(2m)^{2}(d_{v}+1)^{2}-(2m)(2m+2)^{2}A_{vv}+(2m+2)^{2}d_{v}^{2}$
$\displaystyle+2\left(\sum_{v\neq y\in
C_{2}}\left((2m)^{2}(2m+2)A_{vy}-(2m)^{2}(d_{v}+1)d_{y}-(2m)(2m+2)^{2}A_{vy}+(2m+2)^{2}d_{v}d_{y}\right)\right)$
$\displaystyle+\sum_{v\neq x,y\in
C_{2}}\left((2m)^{2}(2m+2)A_{xy}-(2m)^{2}d_{x}d_{y}-(2m)(2m+2)^{2}A_{xy}+(2m+2)^{2}d_{x}d_{y}\right)$
$\displaystyle+\sum_{C_{i}\in T,i\notin\\{1,2\\}}\sum_{x,y\in
C_{i}}\left((2m)^{2}(2m+2)A_{xy}-(2m)^{2}d_{x}d_{y}-(2m)(2m+2)^{2}A_{xy}+(2m+2)^{2}d_{x}d_{y}\right)$
$\displaystyle=(-8m^{2}-8m)A_{uu}-8m^{2}d_{u}-4m^{2}+8md^{2}_{u}+4d^{2}_{u}$
$\displaystyle+2\left(\sum_{u\neq x\in
C_{1}}\left((-8m^{2}-8m)A_{ux}-4m^{2}(d_{u}d_{x}+d_{x})+(4m^{2}+8m+4)d_{u}d_{x}\right)\right)$
$\displaystyle+\sum_{u\neq x,y\in
C_{1}}\left((-8m^{2}-8m)A_{xy}+(8m+4)d_{x}d_{y}\right)$
$\displaystyle+(-8m^{2}-8m)A_{vv}-8m^{2}d_{v}-4m^{2}+8md^{2}_{v}+4d^{2}_{v}$
$\displaystyle+2\left(\sum_{v\neq y\in
C_{2}}\left((-8m^{2}-8m)A_{vy}-4m^{2}(d_{v}d_{y}+d_{y})+(4m^{2}+8m+4)d_{v}d_{y}\right)\right)$
$\displaystyle+\sum_{v\neq x,y\in
C_{2}}\left((-8m^{2}-8m)A_{xy}+(8m+4)d_{x}d_{y}\right)$
$\displaystyle+\sum_{C_{i}\in T,i\notin\\{1,2\\}}\sum_{x,y\in
C_{i}}\left((-8m^{2}-8m)A_{xy}+(8m+4)d_{x}d_{y}\right)$
$\displaystyle=(-8m^{2}-8m)A_{uu}+(8m+4)d^{2}_{u}-8m^{2}d_{u}-4m^{2}$
$\displaystyle+2\left(\sum_{u\neq x\in
C_{1}}\left((-8m^{2}-8m)A_{ux}+(8m+4)d_{u}d_{x}-4m^{2}d_{x}\right)\right)$
$\displaystyle+\sum_{u\neq x,y\in
C_{1}}\left((-8m^{2}-8m)A_{xy}+(8m+4)d_{x}d_{y}\right)$
$\displaystyle+(-8m^{2}-8m)A_{vv}+(8m+4)d^{2}_{v}-8m^{2}d_{v}-4m^{2}$
$\displaystyle+2\left(\sum_{v\neq y\in
C_{2}}\left((-8m^{2}-8m)A_{vy}+(8m+4)d_{v}d_{y}-4m^{2}dy\right)\right)$
$\displaystyle+\sum_{v\neq x,y\in
C_{2}}\left((-8m^{2}-8m)A_{xy}+(8m+4)d_{x}d_{y}\right)$
$\displaystyle+\sum_{C_{i}\in T,i\notin\\{1,2\\}}\sum_{x,y\in
C_{i}}\left((-8m^{2}-8m)A_{xy}+(8m+4)d_{x}d_{y}\right)$
$\displaystyle=-8m^{2}-\sum_{x\in C_{1}}8m^{2}d_{x}-\sum_{y\in
C_{2}}8m^{2}d_{y}+\sum_{C_{i}\in T}\sum_{x,y\in
C_{i}}\left((-8m^{2}-8m)A_{xy}+(8m+4)d_{x}d_{y}\right)$

Thus, the increase in modularity by adding $(u,v)$ is

$\frac{-8m^{2}-\sum_{x\in C_{1}}8m^{2}d_{x}-\sum_{y\in
C_{2}}8m^{2}d_{y}+\sum_{C_{i}\in T}\sum_{x,y\in
C_{i}}\left((-8m^{2}-8m)A_{xy}+(8m+4)d_{x}d_{y}\right)}{(2m)^{2}(2m+2)^{2}}\text{.}$

∎