# A study on modularity density maximization:
Column generation acceleration and computational complexity analysis

Issey Sukeda<EMAIL_ADDRESS>Graduate School of
Information Science and Technology, University of Tokyo, Japan  Atsushi
Miyauchi<EMAIL_ADDRESS>(Corresponding author) Graduate School
of Information Science and Technology, University of Tokyo, Japan  Akiko
Takeda<EMAIL_ADDRESS>Graduate School of Information Science and
Technology, University of Tokyo, Japan  Center for Advanced Intelligence
Project, RIKEN, Japan

###### Abstract

Community detection is a fundamental network-analysis primitive with a variety
of applications in diverse domains. Although the modularity introduced by
Newman and Girvan (2004) has widely been used as a quality function for
community detection, it has some drawbacks. The modularity density introduced
by Li et al. (2008) is known to be an effective alternative to the modularity,
which mitigates one of the drawbacks called the resolution limit. A large body
of work has been devoted to designing exact and heuristic methods for
modularity density maximization, without any computational complexity
analysis. In this study, we investigate modularity density maximization from
both algorithmic and computational complexity aspects. Specifically, we first
accelerate column generation for the modularity density maximization problem.
To this end, we point out that the auxiliary problem appearing in column
generation can be viewed as a dense subgraph discovery problem. Then we employ
a well-known strategy for dense subgraph discovery, called the greedy peeling,
for approximately solving the auxiliary problem. Moreover, we reformulate the
auxiliary problem to a sequence of $0$–$1$ linear programming problems,
enabling us to compute its optimal value more efficiently and to get more
diverse columns. Computational experiments using a variety of real-world
networks demonstrate the effectiveness of our proposed algorithm. Finally, we
show the NP-hardness of a slight variant of the modularity density
maximization problem, where the output partition has to have two or more
clusters, as well as showing the NP-hardness of the auxiliary problem.

##  1 Introduction

Community detection, identifying community structure in networks, is a
fundamental network-analysis primitive with a variety of applications in
diverse domains [26, 28]. Roughly speaking, a community in a network is a
subset of vertices densely connected internally but only sparsely connected to
the rest of the network, which corresponds to a social group in friendship
networks (trivially), a molecular complex in protein–protein interaction
networks [59], a spam link farm in the Web graph [30], and so forth.

Newman and Girvan [52] introduced a quality function for community detection,
called the modularity. Let $G=(V,E)$ be an undirected graph consisting of
$n=|V|$ vertices and $m=|E|$ edges. For a partition $\mathcal{C}$ of $V$
(i.e., $\mathcal{C}=\\{C_{1},\dots,C_{k}\\}$ such that
$\bigcup_{i=1}^{k}C_{i}=V$ and $C_{i}\cap C_{j}=\emptyset$ for $i\neq j$), the
modularity is defined as

$\displaystyle
Q(\mathcal{C})=\sum_{C\in\mathcal{C}}\left(\frac{|E(C)|}{m}-\left(\frac{\sum_{v\in
C}\deg(v)}{2m}\right)^{2}\right),$

where $E(C)$ is the set of edges in the subgraph induced by $C\subseteq V$,
i.e., $E(C)=\\{\\{i,j\\}\in E\mid i,j\in C\\}$, and $\deg(v)$ is the degree of
$v\in V$ in $G$. Intuitively speaking, the modularity quantifies the actual
fraction of the number of edges within the clusters minus the expected value
of such a fraction assuming that the edges are put at random under the same
degree distribution.

The modularity has become the most popular quality function for community
detection; that is, community detection is often conducted through modularity
maximization [26]. Specifically, the modularity maximization problem asks to
find a partition $\mathcal{C}$ of $V$ that maximizes the modularity
$Q(\mathcal{C})$. Unlike most traditional clustering problems (e.g., the
$k$-means problem), the modularity maximization problem does not specify the
number of clusters in advance; it can be determined by algorithms
endogenously. The modularity maximization problem is known to be NP-hard [9].
To date, a large number of exact, approximation, and heuristic methods have
been developed for the problem (see Section 2).

However, it is also known that the modularity has some drawbacks as a quality
function for community detection, such as the degeneracy [35] and the
resolution limit [27]. The former means that there are numerous solutions that
have modularity values very close to the optimal value, causing the
significant difficulty of computing an optimal solution (even in exponential
time). The latter means that the modularity fails to detect a community even
if it is a relatively large clique connected by only a few edges with the rest
of the network, when the size of the whole network becomes large.

To overcome this issue, Li et al. [44] introduced a novel quality function,
which they refer to as the modularity density. For a partition $\mathcal{C}$
of $V$, the modularity density is defined as

$\displaystyle
D(\mathcal{C})=\sum_{C\in\mathcal{C}}\frac{2|E(C)|-|E(C,V\setminus C)|}{|C|},$

where $E(C,V\setminus C)$ is the set of edges connecting a vertex in $C$ and a
vertex in $V\setminus C$. The modularity density quantifies the following
value for each cluster and takes the summation: the average degree to the
inside of the cluster minus the average degree to the outside of the cluster.
The modularity density is known to be an effective alternative to the
modularity; in fact, it mitigates one of the aforementioned drawbacks, the
resolution limit. The modularity density maximization problem asks to find a
partition $\mathcal{C}$ of $V$ that maximizes $D(\mathcal{C})$, which again
does not specify the number of clusters in advance.

A large body of work has been devoted to designing exact and heuristic methods
for the modularity density maximization problem. Exact methods have been
actively studied from mathematical programming point of view [18, 20, 57].
Costa [18] and Costa et al. [20] designed exact methods using mixed-integer-
linear-programming reformulation techniques. Sato and Izunaga [57] developed a
branch-and-price method, i.e., a branch-and-bound method that incorporates a
column generation algorithm for solving a linear programming relaxation on
each node in the search tree. On the other hand, heuristic methods have been
designed based on various techniques such as local search [55], spectral
optimization [8], simulated annealing [65], mathematical programming [19, 21,
37], post processing [58], and others [54]. In particular, the algorithm
designed by de Santiago and Lamb [21] is the first column generation algorithm
for the modularity density maximization problem, which outputs an optimal
solution if the linear programming relaxation used has an integral optimal
solution (and it is often the case).

Unlike the modularity maximization problem, there are no computational
complexity results for the modularity density maximization problem; in
particular, it is not known whether the modularity density maximization
problem is NP-hard, as mentioned by Costa et al. [20]. Some literature (e.g.,
Li et al. [44]) states that the problem is NP-hard, but there is no valid
support for the statement.

###  1.1 Our contribution

In this study, we investigate modularity density maximization from both
algorithmic and computational complexity aspects.

Specifically, we first accelerate column generation for the modularity density
maximization problem. Column generation is a well-known strategy for solving
(linear programming relaxations of) combinatorial optimization problems. For
the modularity density maximization problem, as mentioned above, de Santiago
and Lamb [21] established a column generation framework and later Sato and
Izunaga [57] incorporated it into a branch-and-bound method to design a
branch-and-price method.

To describe our contribution, we review the mechanism of column generation in
a general manner. First the targeted combinatorial optimization problem is
formulated as an integer linear programming problem with exponential number of
variables, and then its linear programming relaxation is obtained by relaxing
the integral constraints. It is empirically known that such a relaxation tends
to have an integral optimal solution (e.g., see de Santiago and Lamb [21] for
the case of the modularity density maximization problem); therefore, in that
case, the goal can usually be achieved by solving the relaxation. However, as
the relaxation also has an exponential number of (continuous) variables, it is
still hard to solve (even store) the relaxation. To overcome this issue, the
dual problem of the relaxation is then considered, which has an exponential
number of constraints instead. Column generation first solves the dual problem
with only a part of constraints and obtains an optimal solution to the
subproblem. Then it solves an optimization problem called the auxiliary
problem, which finds a constraint in the (full) dual problem violated by the
incumbent solution to the most extent (if exists). If there is no violated
constraint found, it can be seen that the incumbent solution is an optimal
solution to the (full) dual problem and we are done. Otherwise we add the
constraint found to the subproblem, solve it to obtain an updated solution,
and repeat the above procedure. The auxiliary problem is often computationally
expensive; therefore, column generation algorithms usually solve the auxiliary
problem just approximately using some heuristic methods, and solve it exactly
only if the heuristic methods fail to find a violated constraint. In the
context of the modularity density maximization problem, de Santiago and Lamb
[21] devised a simple local search heuristic for solving the auxiliary problem
approximately and introduced $0$–$1$ quadratic programming formulations for
solving it exactly.

In this study, to design an effective heuristic method for the auxiliary
problem appearing in column generation for the modularity density maximization
problem, we point out that the auxiliary problem can be seen as a dense
subgraph discovery problem. Dense subgraph discovery, aiming to find a
(single) dense component in a network, is a well-studied task in the data
mining community (see Section 2). An important fact is that there is a well-
known strategy for designing effective dense subgraph discovery algorithms,
called the greedy peeling (see e.g., [7, 12, 13, 40, 41, 43, 48, 60, 61, 62]).
In this strategy, a vertex with the minimum contribution (e.g., the minimum
degree) in the current remaining graph is removed iteratively to produce a
sequence of vertex subsets from $V$ to a singleton, and then the best subset
among the sequence (in terms of the objective function at hand) is returned.
For the classic dense subgraph discovery problem, called the densest subgraph
problem, this strategy is known to produce a $1/2$-approximation algorithm
[12, 41]. We apply the greedy peeling to the auxiliary problem in a
sophisticated way.

Moreover, to solve the auxiliary problem exactly (when the above heuristic
method fails to find a violated constraint), we reformulate the auxiliary
problem to a sequence of $0$–$1$ linear programming problems, indexed by the
solution size (from $1$ to $n$), where one of the optimal solutions to the
problems (with the highest objective function value) corresponds to an optimal
solution to the auxiliary problem. This formulation has two advantages over
the above existing $0$–$1$ quadratic programming formulations. First, due to
the linearity of the objective functions, it can compute the optimal value to
the auxiliary problem more efficiently (despite the increase of the number of
problems to solve). Second, it computes not only an optimal solution to the
auxiliary problem but also quasi-optimal solutions with various sizes (from
$1$ to $n$), which often correspond to violated constraints in the (full) dual
problem.

Computational experiments using a variety of real-world networks demonstrate
the effectiveness of our proposed column generation algorithm. In particular,
our algorithm obtains an optimal solution within only two hours to an instance
with more than 100 vertices that cannot be solved by the state-of-the-art
methods in 24 hours.

Finally, we show the NP-hardness of a slight variant of the modularity density
maximization problem, where the output partition has to have two or more
clusters. Therefore, the variant rules out the trivial clustering consisting
of the single cluster containing all vertices. Although the NP-hardness of the
variant does not lead to the NP-hardness of the original problem (i.e., an
open problem mentioned by e.g, [20, 19, 15, 17]), ours is the first result to
verify the hardness of maximizing the modularity density function. In our
proof, we construct a polynomial-time reduction from the well-known maximum
cut problem (on $3$-regular graphs) to the above variant of the modularity
density maximization problem (on regular graphs). Our reduction is based on an
existing reduction from the maximum cut problem to the uniform sparsest cut
problem [6]; however, the difference from the existing reduction is not
totally trivial. As a by-product of our reduction, we also show that the
uniform sparsest cut problem is NP-hard even on regular graphs. In addition,
we also prove the NP-hardness of the above auxiliary problem appearing in
column generation, which justifies the use of effective heuristic methods for
it.

###  1.2 Paper organization

The rest of the paper is organized as follows: In Section 2, we review some
related work. In Section 3, we describe a general column generation framework
for the modularity density maximization problem. In Section 4, we present our
techniques for accelerating column generation. We conduct the computational
experiments and evaluate the performance of our proposed algorithm in Section
5. We show the NP-hardness of the variant of the modularity density
maximization problem and the auxiliary problem in Section 6. We conclude the
study and suggest future directions in Section 7.

##  2 Related Work

The modularity density is known to be an effective alternative of the original
modularity. de Santiago and Lamb [22] have recently conducted a systematic
comparison between the modularity and the modularity density, using exact
solutions for both quality functions. The modularity density itself also has
many variations. For example, Chen et al. [15] pointed out that the modularity
density still has the (milder) resolution limit, and developed a variant for
further mitigating it. Chen et al. [14] presented a variant suitable for
overlapping community detection, where each vertex may belong to more than one
cluster. Yan et al. [63] introduced a variant applicable to community
detection in signed networks [63], where each edge has a positive/negative
label.

For the modularity maximization problem, Aloise et al. [4] developed a column
generation algorithm. They used an effective heuristic method for the targeted
auxiliary problem, based on a variable neighborhood search, and succeeded in
solving an instance with more than 500 vertices in reasonable time. The design
of exact methods for the modularity maximization problem is based on a simple
integer linear programming formulation [1, 25]. This can be done due to the
fact that the modularity can be rewritten as a linear form using 0–1 variables
over the pairs of vertices. Such a property cannot be observed for the
modularity density, which makes it difficult to design such a simple exact
algorithm for the modularity density maximization problem.

There are several approximation algorithms for the modularity maximization
problem. Dinh et al. [24] developed a semidefinite-programming-based
approximation algorithm, and later Kawase et al. [39] improved the algorithm
and analysis. Heuristic methods for the modularity maximization problem are
based on various techniques such as greedy procedure [5, 16], spectral
optimization [51, 50, 53], simulated annealing [36, 46], and mathematical
programming [1, 10, 11, 49], some of which have inspired the design of some
heuristics for the modularity density maximization problem. The most popular
heuristic method is the so-called Louvain method [5], which works well in
practice in terms of both solution quality and computation time.

Dense subgraph discovery is a well-known graph-mining task with a variety of
applications [31]. The most well-studied optimization problem for dense
subgraph discovery would be the densest subgraph problem. For a vertex subset
$S\subseteq V$, the degree density (or simply density) is defined as
$d(S)=\frac{|E(S)|}{|S|}$, which is equal to half the average degree of the
subgraph induced by $S$. The densest subgraph problem asks to find $S\subseteq
V$ that maximizes $d(S)$ without any size constraint. Unlike most traditional
optimization problem for dense subgraph discovery such as the maximum clique
problem, the densest subgraph problem can be solved exactly in polynomial
time; indeed, there are exact algorithms based on a maximum-flow-based
technique [34] or a linear-programming-based technique [12]. As mentioned in
the introduction, the greedy peeling strategy produces a $1/2$-approximation
algorithm for the problem [12, 41], which can be implemented to run in linear
time, using some sophisticated data structure.

The densest subgraph problem has many other problem variations. The most well-
studied variants are size-restricted ones, e.g., the densest $k$-subgraph
problem, where given $G=(V,E)$ and positive integer $k$, we are asked to find
$S\subseteq V$ that maximizes $d(S)$ (or simply $|E(S)|$) under the size
constraint $|S|=k$. It is known that such a restriction makes the problem
harder to solve; in fact, the densest $k$-subgraph problem is NP-hard (the
maximum clique problem is trivially reducible). A variant recently introduced
by Miyauchi and Kakimura [48] is most related to this study. Their variant
asks to find $S\subseteq V$ that maximizes
$f_{\alpha}(S)=\frac{|E(S)|-\alpha|E(S,V\setminus S)|}{|S|}$, where
$\alpha\in[0,\infty)$ is a constant. Thus, setting $\alpha=1/2$, we see that
the optimization problem becomes the one-term variant of the modularity
density maximization problem, which has a strong connection to the auxiliary
problem appearing in column generation (see Section 3). For the above problem
with general $\alpha\in[0,\infty)$, the authors developed polynomial-time
exact and approximation algorithms, based on those for the original densest
subgraph problem. The acceleration techniques developed in our column
generation algorithm are inspired by their algorithm.

##  3 Column Generation Framework

In this section, we review a general column generation framework for the
modularity density maximization problem. The following framework is due to de
Santiago and Lamb [21], although there is a slight difference in the notations
used.

Let $G=(V,E)$ be an undirected graph consisting of $n=|V|$ vertices and
$m=|E|$ edges. Let $\mathcal{C}$ be a partition of $V$. We define the
contribution of $C\in\mathcal{C}$ as

$\displaystyle c(C)=\frac{2|E(C)|-|E(C,V\setminus C)|}{|C|},$

which is nothing but the term of $C\in\mathcal{C}$ in the modularity density
$D(\mathcal{C})$. Therefore, we can write
$D(\mathcal{C})=\sum_{C\in\mathcal{C}}c(C)$. In the column generation
framework, the modularity density maximization problem is formulated as a 0–1
linear programming problem with exponential number of variables, by viewing
the problem as the set partitioning problem. Let
$\mathcal{S}=2^{V}\setminus\emptyset$ be the set of all (nonempty) subsets of
$V$. For $v\in V$ and $S\in\mathcal{S}$, let $a_{vS}$ be a constant equal to 1
if $v\in S$ and 0 otherwise. Let $z_{S}$ be a 0–1 variable that takes 1 if
$S\subseteq V$ is used in the solution and 0 otherwise. Then the modularity
density maximization problem can be formulated as follows:

maximize $\displaystyle\sum_{S\in\mathcal{S}}c(S)z_{S}$ subject to
$\displaystyle\sum_{S\in\mathcal{S}}a_{vS}z_{S}=1$ $\displaystyle(v\in V),$
$\displaystyle z_{S}\in\\{0,1\\}$ $\displaystyle(S\in\mathcal{S}).$

The first constraint stipulates that every vertex belongs to exactly one
cluster. This formulation has an exponential number of 0–1 variables. By
relaxing the binary constraint, we can obtain the following linear programming
relaxation, which we refer to as the master problem (MP):

maximize $\displaystyle\sum_{S\in\mathcal{S}}c(S)z_{S}$ subject to
$\displaystyle\sum_{S\in\mathcal{S}}a_{vS}z_{S}=1,$ $\displaystyle(v\in V),$
$\displaystyle z_{S}\geq 0$ $\displaystyle(S\in\mathcal{S}).$

Note that in this relaxation, the constraint $z_{S}\leq 1$ ($S\in\mathcal{S}$)
is omitted because it is automatically satisfied due to the first constraint
together with the definition of $a_{vS}$. In general, the MP does not
necessarily admit an integral solution, but according to the computational
experiments of de Santiago and Lamb [21], the optimal solution tends to have
the integrality for most instances. Therefore, we aim to solve the MP instead
of the original problem. However, the MP in turn has an exponential number of
(continuous) variables; thus, it is still hard to solve (even store) the
problem.

To overcome this issue, the dual problem (DP) of the MP is then considered:

minimize $\displaystyle\sum_{v\in V}\lambda_{v}$ subject to
$\displaystyle\sum_{v\in V}a_{vS}\lambda_{v}\geq c(S)$
$\displaystyle(S\in\mathcal{S}),$ $\displaystyle\lambda_{v}\in\mathbb{R}$
$\displaystyle(v\in V).$

Notice that the DP has an exponential number of constraints rather than
variables. Column generation first solves the DP with only a part of
constraints. Specifically, taking a tiny subset
$\mathcal{S}^{\prime}\subseteq\mathcal{S}$, we solve the following relaxed
dual problem (RDP):

minimize $\displaystyle\sum_{v\in V}\lambda_{v}$ subject to
$\displaystyle\sum_{v\in V}a_{vS}\lambda_{v}\geq c(S)$
$\displaystyle(S\in\mathcal{S}^{\prime}\subseteq\mathcal{S}),$
$\displaystyle\lambda_{v}\in\mathbb{R}$ $\displaystyle(v\in V),$

and obtain its optimal solution $\bm{\lambda^{*}}=(\lambda^{*}_{v})_{v\in V}$.
Then we solve an optimization problem called the auxiliary problem (AP), which
can be formulated as follows:

minimize $\displaystyle\sum_{v\in V}a_{vS}\lambda^{*}_{v}-c(S)$ subject to
$\displaystyle S\in\mathcal{S}.$

We see that the AP is an optimization problem that finds a constraint in the
DP violated by the incumbent optimal solution $\bm{\lambda^{*}}$ to the most
extent (if exists). Indeed, if there exists a solution $S\in\mathcal{S}$ with
the objective function value less than zero, then the constraint in the DP
corresponding to $S$ is violated by $\bm{\lambda^{*}}$. Otherwise we see that
$\bm{\lambda^{*}}$ satisfies all constraints in the DP and it is an optimal
solution to the DP. However, the AP is a computationally expensive problem; in
fact, in Section 6, we prove that the AP is NP-hard even on $(n-4)$-regular
graphs (when regarding $\bm{\lambda^{*}}$ as a part of the input). Therefore,
the existing column generation algorithm by de Santiago and Lamb [21] employs
a simple local search heuristic for solving the AP approximately. Their
algorithm solves the AP exactly using $0$–$1$ quadratic programming
formulations only if the heuristic method fails to find a violated constraint.

##  4 The Proposed Method

In this section, we present our proposed techniques for accelerating column
generation for the modularity density maximization problem.

###  4.1 Greedy peeling algorithm for approximately solving the auxiliary
problem

Recalling that for $v\in V$ and $S\in\mathcal{S}$, $a_{vS}=1$ if and only if
$v\in S$ and multiplying the objective function of the AP by $-1$, we can
rewrite the AP as follows:

maximize $\displaystyle c(S)-\sum_{v\in S}\lambda^{*}_{v}$ subject to
$\displaystyle S\in\mathcal{S}.$

Here we wish to point out that the AP can be seen as a dense subgraph
discovery problem. The most well-studied optimization problem for dense
subgraph discovery would be the densest subgraph problem, which asks to find
$S\subseteq V$ that maximizes the density $d(S)=\frac{|E(S)|}{|S|}$. As
mentioned in Section 2, Miyauchi and Kakimura [48] introduced a generalization
of the density: $f_{\alpha}(S)=\frac{|E(S)|-\alpha|E(S,V\setminus S)|}{|S|}$,
where $\alpha\in[0,\infty)$ is a constant, and studied the optimization
problem that asks to maximize $f_{\alpha}(S)$. As the first term of the
objective function of the AP, i.e., $c(S)$, is nothing but $2\cdot
f_{1/2}(S)$, the AP can be seen as a variant of the optimization problem
studied by Miyauchi and Kakimura [48], where each vertex $v\in V$ has a weight
$\lambda^{*}_{v}\in\mathbb{R}$ that is subtracted if $v$ is chosen as an
element of $S$.

The greedy peeling is a well-known strategy for designing effective algorithms
for dense subgraph discovery. In this strategy, a vertex with the minimum
contribution (e.g., the minimum degree) in the current remaining graph is
removed iteratively to produce a sequence of vertex subsets from $V$ to a
singleton, and then the best subset among the sequence (in terms of the
objective function at hand) is returned. For $S\subseteq V$ and $v\in S$, let
$\deg_{S}(v)$ be the degree of $v$ in the induced subgraph $G[S]=(S,E(S))$.
Then for the densest subgraph problem, the algorithm based on the greedy
peeling strategy can be written as Algorithm 1.

Input : $G=(V,E)$

Output : $S\subseteq V$

1 $S_{n}\leftarrow V$, $i\leftarrow n$;

2 while  _$i >1$_ do

3 $v_{\text{min}}\leftarrow\mathop{\rm argmin}\\{\deg_{S_{i}}(v)\mid v\in
S_{i}\\}$;

4 $S_{i-1}\leftarrow S_{i}\setminus\\{v_{\text{min}}\\}$, $i\leftarrow i-1$;

5

6$S_{\mathrm{max}}\leftarrow\mathop{\rm argmax}\\{d(S)\mid
S\in\\{S_{1},\dots,S_{n}\\}\\}$;

7 return $S_{\mathrm{max}}$;

Algorithm 1 Greedy peeling algorithm for the densest subgraph problem

It is known that Algorithm 1 is a $1/2$-approximation algorithm for the
densest subgraph problem [12, 41]. Miyauchi and Kakimura [48] employed the
greedy peeling for their generalized problem and demonstrated that the
algorithm performs well in practice. As the AP is a variant of their
optimization problem, it would be reasonable to apply the greedy peeling to
the AP.

From now on, we describe our greedy peeling algorithm for the AP. Let $g$ be
the objective function of the AP to be maximized, i.e.,

$\displaystyle g(S)=c(S)-\sum_{v\in S}\lambda^{*}_{v}\quad\text{ for
}S\subseteq V.$

When solving the AP approximately, we wish to obtain not only a single
solution but also diverse solutions corresponding to various constraints
violated by the incumbent optimal solution $\bm{\lambda^{*}}$ to the RDP. To
this end, we introduce a more general objective function $g_{p}$ ($p\in[0,1]$)
by taking a convex combination of the first and second terms of $g$, i.e., for
$S\subseteq V$,

$\displaystyle g_{p}(S)$ $\displaystyle=p\cdot c(S)+(1-p)\left(-\sum_{v\in
S}\lambda^{*}_{v}\right)=\frac{p(2|E(S)|-|E(S,V\setminus
S)|)-(1-p)|S|\sum_{v\in S}\lambda^{*}_{v}}{|S|}.$ (1)

If we take a large value of $p$, the function emphasizes the first term
$c(S)$, whereas if we take a small value of $p$, the function emphasizes the
second term $-\sum_{v\in S}\lambda^{*}_{v}$. It is worth mentioning that
considering a general objective function like the above might lead to a better
solution even in terms of the original objective function $g$, because the
greedy peeling algorithm is not an exact algorithm (even for the case of the
densest subgraph problem).

The greedy peeling strategy starts with the whole graph and iteratively
removes a vertex with the minimum contribution in the current remaining graph.
Our greedy peeling algorithm introduces two measures for computing
contribution of vertices. Let $S\subseteq V$. When removing a vertex from $S$,
the change of the denominator of (1) does not depend on the choice of a
vertex. Therefore, it suffices to take into account only the numerator of (1).
Our first measure comes from the following transformation:

$\displaystyle\text{The numerator of }\eqref{eq:objective}=\sum_{v\in
S}\left(p\left(\deg_{S}(v)-\deg_{V\setminus(S\setminus\\{v\\})}(v)\right)-(1-p)|S|\lambda^{*}_{v}\right).$

Note that $\deg_{V\setminus(S\setminus\\{v\\})}(v)$ represents the degree of
$v\in S$ to the outside of $G[S]$. It would be reasonable to regard the term
of $v\in S$ as a contribution of $v$, leading to our first measure:

$\displaystyle\textsf{cont}^{\textsf{sum}}_{p}(v)=p\left(\deg_{S}(v)-\deg_{V\setminus(S\setminus\\{v\\})}(v)\right)-(1-p)|S|\lambda^{*}_{v}\quad\text{
for }v\in S.$

Our second measure focuses on the actual difference between the objective
function values before and after removing a vertex. For any $u\in S$, when we
remove $u$ from $S$, the numerator of (1) decreases by

$\displaystyle p(2|E(S)|-|E(S,V\setminus S)|)-(1-p)|S|\sum_{v\in
S}\lambda^{*}_{v}$
$\displaystyle\quad-\left(p(2|E(S\setminus\\{u\\})|-|E(S\setminus\\{u\\},V\setminus(S\setminus\\{u\\}))|)-(1-p)|S\setminus\\{u\\}|\sum_{v\in
S\setminus\\{u\\}}\lambda^{*}_{v}\right)$
$\displaystyle=p\left(3\deg_{S}(u)-\deg_{V\setminus(S\setminus\\{u\\})}(u)\right)-(1-p)\left((|S|-1)\lambda^{*}_{u}+\sum_{v\in
S}\lambda^{*}_{v}\right).$

Noticing that the last term $\sum_{v\in S}\lambda^{*}_{v}$ is independent of
the choice of $u\in S$, we can obtain our second measure:

$\displaystyle\textsf{cont}^{\textsf{diff}}_{p}(v)=p\left(3\deg_{S}(v)-\deg_{V\setminus(S\setminus\\{v\\})}(v)\right)-(1-p)(|S|-1)\lambda^{*}_{v}\quad\text{
for }v\in S.$

To increase the diversity of solutions to be obtained, we introduce a more
general measure by again taking a convex combination of the above two
measures:

$\displaystyle\textsf{cont}_{p,q}(v)=q\cdot\textsf{cont}^{\textsf{sum}}_{p}(v)+(1-q)\cdot\textsf{cont}^{\textsf{diff}}_{p}(v)\quad\text{
for }v\in S.$

Our greedy peeling algorithm runs as follows: Let $P$ and $Q$ be the sets of
candidates for the above parameters $p$ and $q$, respectively. For each $p\in
P$ and each $q\in Q$, the algorithm starts with the whole graph and
iteratively removes a vertex with the minimum contribution
$\textsf{cont}_{p,q}(v)$ in the current remaining graph to obtain a sequence
of vertex subsets from $V$ to a singleton. Unlike the original greedy peeling
strategy, which focuses only on the best subset among the sequence, we store
all subsets $S\subseteq V$ among the sequence that satisfy $g(S)>0$ to find as
many constraints as possible in the DP violated by $\bm{\lambda^{*}}$. For
reference, the procedure is described in Algorithm 2.

Input : $G=(V,E)$, $\bm{\lambda^{*}}=(\lambda^{*}_{v})_{v\in V}$

Output : A family $\mathcal{S}$ of vertex subsets

1 $\mathcal{S}\leftarrow\emptyset$;

// The sets $P$ and $Q$ of candidates for the parameters $p$ and $q$,
respectively, are fixed in advance.

2 for  _each $p\in P$_ do

3 for  _each $q\in Q$_ do

4 $S_{n}\leftarrow V$, $i\leftarrow n$;

5 while  _$i >1$_ do

6 if  _$g(S_{i}) >0$_ then

7 $\mathcal{S}\leftarrow\mathcal{S}\cup\\{S_{i}\\}$;

8 $v_{\text{min}}\leftarrow\mathop{\rm argmin}\\{\textsf{cont}_{p,q}(v)\mid
v\in S_{i}\\}$;

9 $S_{i-1}\leftarrow S_{i}\setminus\\{v_{\text{min}}\\}$, $i\leftarrow i-1$;

10

11

12

13return $\mathcal{S}$;

Algorithm 2 $\textsf{Peeling}(G,\bm{\lambda^{*}})$

It is easy to see that Algorithm 2 runs in $O(|P||Q|n^{2})$ time, that is,
each iteration for a pair of $p$ and $q$, corresponding to a usual greedy
peeling algorithm, consumes $O(n^{2})$ time. On the other hand, the most
similar algorithm to ours, i.e., the greedy peeling algorithm by Miyauchi and
Kakimura [48], runs in $O(m+n\log n)$ time. This difference is due to the
change of the definition of contribution. In the algorithm by Miyauchi and
Kakimura [48], when a vertex is removed, the contribution values are changed
only for its neighbors. Therefore, using the Fibonacci heap, we can obtain the
running time of $O(m+n\log n)$. On the other hand, in our algorithm, the
contribution values are changed for all vertices, due to the term of the
current subset size $|S|$, which makes it difficult to reduce the running time
of $O(n^{2})$.

###  4.2 Reformulation of the auxiliary problem

Here we provide a new reformulation of the AP. We deal with the original
minimization version of the AP:

minimize $\displaystyle\sum_{v\in V}a_{vS}\lambda^{*}_{v}-c(S)$ subject to
$\displaystyle S\in\mathcal{S}.$

Note that for any $S\subseteq V$, we have

$\displaystyle c(S)=\frac{2|E(S)|-|E(S,V\setminus
S)|}{|S|}=\frac{4|E(S)|-\sum_{v\in S}\deg(v)}{|S|}.$

For each $e\in E$, let $x_{e}$ be a $0$–$1$ variable that takes $1$ if both of
the endpoints of $e$ are included in the solution and 0 otherwise. For each
$v\in V$, let $y_{v}$ be a $0$–$1$ variable that takes $1$ if $v$ is included
in the solution and $0$ otherwise. Then the AP can be formulated as the
following $0$–$1$ fractional programming problem:

minimize $\displaystyle\sum_{v\in V}\lambda_{v}^{*}y_{v}-\frac{\sum_{e\in
E}4x_{e}-\sum_{v\in V}\deg(v)y_{v}}{\sum_{v\in V}y_{v}}$ subject to
$\displaystyle x_{e}\leq y_{u},\ x_{e}\leq y_{v}$
$\displaystyle(e=\\{u,v\\}\in E),$ $\displaystyle x_{e}\in\\{0,1\\}$
$\displaystyle(e\in E),$ $\displaystyle y_{v}\in\\{0,1\\}$ $\displaystyle(v\in
V).$

The fractionality appears in the second term of the objective function,
corresponding to the contribution of the solution in terms of the modularity
density. To remove this, we consider fixing the size of solutions $\sum_{v\in
V}y_{v}$ to each $k=1,2,\dots,n$. Specifically, for each $k=1,2,\dots,n$, we
introduce the following $0$–$1$ linear programming problem ($\text{AP}(k)$):

minimize $\displaystyle\sum_{v\in V}\lambda_{v}^{*}y_{v}-\frac{\sum_{e\in
E}4x_{e}-\sum_{v\in V}\deg(v)y_{v}}{k}$ subject to $\displaystyle\sum_{v\in
V}y_{v}=k,$ $\displaystyle x_{e}\leq y_{u},\ x_{e}\leq y_{v}$
$\displaystyle(e=\\{u,v\\}\in E),$ $\displaystyle x_{e}\in\\{0,1\\}$
$\displaystyle(e\in E),$ $\displaystyle y_{v}\in\\{0,1\\}$ $\displaystyle(v\in
V).$

An optimal solution to $\text{AP}(k)$ is an optimal solution to the AP with
the size constraint $|S|=k$. Therefore, the best solution among the optimal
solutions to $\text{AP}(k)$ ($k=1,2,\dots,n$) is optimal to the AP, and its
objective value is equal to the optimal value of the AP. Owing to this
reformulation, the AP becomes a sequence of $0$–$1$ linear programming
problems.

Our reformulation has the following advantages: It is expected that an optimal
solution to the AP can be found more quickly than the $0$–$1$ quadratic
programming formulations by de Santiago and Lamb [21]. In fact, $\text{AP}(k)$
is no longer nonlinear and hence easy to handle using MIP solvers, despite of
the increase of the number of problems from one to $n$. Moreover, the
reformulation is useful also in terms of finding diverse constraints in the DP
violated by $\bm{\lambda^{*}}$. Indeed, in the process of solving the
reformulation, we can obtain optimal solutions to $\text{AP}(k)$
($k=1,2,\dots,n$), which can be seen as quasi-optimal solutions to the AP,
perhaps corresponding to some violated constraints in the DP.

Our algorithm for solving the AP exactly using the reformulation is described
in Algorithm 3. In general, MIP solvers employ a branch-and-bound method to
find an optimal solution, where many feasible solutions (with good objective
function values) are found. To obtain diverse constraints, we collect all
vertex subsets $S\subseteq V$ that appeared as feasible solutions for
$\text{AP}(k)$ ($k=1,2,\dots,n$) with negative objective function values.

Input : $G=(V,E)$, $\bm{\lambda^{*}}=(\lambda^{*}_{v})_{v\in V}$

Output : A family $\mathcal{S}$ of vertex subsets

1 $\mathcal{S}\leftarrow\emptyset$;

2 for  _each $k=1,2,\dots,n$_ do

3 Solve $\text{AP}(k)$ and let $\mathcal{S}_{k}$ be the family of vertex
subsets that appear as feasible solutions for $\text{AP}(k)$ with negative
objective function values;

4 $\mathcal{S}\leftarrow\mathcal{S}\cup\mathcal{S}_{k}$;

5

6return $\mathcal{S}$;

Algorithm 3 $\textsf{Exact}(G,\bm{\lambda^{*}})$

###  4.3 Summary: Column generation algorithm using our proposed techniques

We are now ready to present our column generation algorithm for the modularity
density maximization problem. To make the RDP have an optimal solution (i.e.,
bounded), we have to set an initial $\mathcal{S}^{\prime}$ in the RDP
appropriately. One trivial way is to use the singletons, i.e.,
$\mathcal{S}^{\prime}=\\{\\{v\\}\mid v\in V\\}$. Note that the output of any
algorithm for the modularity density maximization problem can be used as an
initial $\mathcal{S}^{\prime}$. Our algorithm first solves the RDP with
$\mathcal{S}^{\prime}$ and obtains its optimal solution $\bm{\lambda^{*}}$.
Then it performs the greedy peeling procedure (Algorithm 2). If new violated
constraints are found, they are added to the RDP and the process is repeated
from the beginning; otherwise the algorithm conducts the exact procedure
(Algorithm 3). For reference, the entire procedure is summarized in Algorithm
4.

Input : $G=(V,E)$

Output : An optimal solution $\bm{\lambda^{*}}=(\lambda^{*}_{v})_{v\in V}$ to
the DP

1 $\mathcal{S}^{\prime}\leftarrow\\{\\{v\\}\mid v\in V\\}$;

// Initialization of the set of constraints in the RDP.

2 while  _TRUE_ do

3 Solve the RDP with $\mathcal{S}^{\prime}$ and obtain its optimal solution
$\bm{\lambda^{*}}=(\lambda^{*}_{v})_{v\in V}$;

4
$\mathcal{S}_{\text{peeling}}\leftarrow\textsf{Peeling}(G,\bm{\lambda^{*}})$;

5 if  _$\mathcal{S}_{\mathrm{peeling}}=\emptyset$_ then

6 $\mathcal{S}_{\text{exact}}=\textsf{Exact}(G,\bm{\lambda}^{*})$;

7 if  _$\mathcal{S}_{\mathrm{exact}}=\emptyset$_ then

8 return $\bm{\lambda}^{*}$;

9

10
$\mathcal{S^{\prime}}\leftarrow\mathcal{S^{\prime}}\cup\mathcal{S}_{\text{exact}}$;

11 continue ;

12

13
$\mathcal{S^{\prime}}\leftarrow\mathcal{S^{\prime}}\cup\mathcal{S}_{\text{peeling}}$;

14

Algorithm 4 Our column generation algorithm

##  5 Experimental Evaluation

In this section, we conduct computational experiments using a variety of real-
world graphs and evaluate the performance of our proposed method.

###  5.1 Experimental setup

Here we explain our experimental setup.

####  5.1.1 Instances

Table 1: Real-world graphs used in our experiments.

ID | Name | $n$ | $m$ | Best $D(\mathcal{C})$ | Description
---|---|---|---|---|---
1 | Strike | 24 | 38 |  $8.86111^{*}$ [18] | Social network [47]
2 | Karate | 34 | 78 |  $7.8451^{*}$ [18] | Social network [64]
3 | Dolphins | 62 | 159 |  $12.1252^{*}$ [21] | Animal social network [45]
4 | Les Misérables | 77 | 254 |  $24.5474^{*}$ [57] | Co-appearance network [45]
5 | Polbooks | 105 | 441 |  $21.9652^{*}$ [21] | Co-purchased network [42]
6 | Adjnoun | 112 | 425 | 7.651 [54] | Word adjacency network [50]
7 | Football | 115 | 613 | 44.340 [54] | Sports game network [32]
8 | Jazz | 198 | 2,742 | 49.716 [54] | Social network [33]

Table 1 summarizes the instances on which our experiments were conducted. The
first and second columns indicate IDs and names, respectively, of the
instances. The third and fourth columns present the number of vertices and the
number of edges, respectively. The fifth column reports the best previously-
known modularity density value for the instances, where the value is
associated with “$*$” if it is known to be the optimal value. The last column
presents the description of the instances.

####  5.1.2 Baseline methods

We employ powerful existing algorithms designed by de Santiago and Lamb [21]
and Sato and Izunaga [57], as baseline methods.

The algorithms by de Santiago and Lamb [21] are the first series of column
generation algorithms for the modularity density maximization problem. The
design of the algorithms follows the column generation framework reviewed in
Section 3, but there are some additional techniques for accelerating the
computation. Among a variety of algorithms they developed, the algorithm
called CGII+ILS is found to be the best algorithm. CGII+ILS solves the AP
using the formulation called AP-II. AP-II is a 0–1 convex quadratic
programming formulation derived by linearizing the nonconvex terms of the
previous formulation called AP-I, which is a 0–1 nonconvex quadratic
programming formulation produced by multiplying the objective function of the
AP by the size term $\sum_{v\in V}y_{v}$. It should be noted that the
transformation from the AP to AP-I (and AP-II) does change the ordering of
solutions in terms of the objective function value; however, it does not
change the sign of the objective function values of solutions and thus AP-I
and AP-II can be used for deciding if the incumbent solution
$\bm{\lambda^{*}}$ of the RDP is feasible for the DP. The term ILS of CGII+ILS
means that the algorithm employs a heuristic method called the iterative local
search for approximately solving the AP. The method starts with a randomly
chosen $S\subseteq V$. Then it repeats the following process for appropriate
number of times: Along with a random ordering of vertices produced, it moves a
vertex from or to the cluster if the movement improves the objective function
value, and then perturbs the cluster if such a movement does not happen. To
construct an initial $\mathcal{S}^{\prime}$ in the RDP, CGII+ILS uses the
method called the hybrid local search, which can be seen as a modularity
density counterpart of the Louvain method [5] followed by an additional local
search procedure.

The algorithm by Sato and Izunaga [57] is the state-of-the-art exact algorithm
for the modularity density maximization problem. The algorithm is based on a
branch-and-price method, i.e., a branch-and-bound method that incorporates a
column generation algorithm for solving a linear programming relaxation for
each node in the search tree. The authors incorporated two existing techniques
into the method, called the set-packing relaxation [56] and the multiple-
cutting-plane-at-a-time [38]. As an initial $\mathcal{S}^{\prime}$ in the RDP,
the authors just employed the singletons, i.e.,
$\mathcal{S}^{\prime}=\\{\\{v\\}\mid v\in V\\}$. An important fact is that
unlike the algorithms by de Santiago and Lamb [21] and even our proposed
algorithm, the algorithm by Sato and Izunaga [57] is guaranteed to output an
optimal solution to the modularity density maximization problem (rather than
its relaxation) at the termination. Therefore, comparing the performance of
those two types of algorithms seems unfair. However, the algorithm by Sato and
Izunaga [57] can be a baseline in the sense that our proposed algorithm needs
to have greater scalability than their algorithm.

####  5.1.3 Implementation and machine

Our proposed algorithm was implemented in Python 3. Throughout the
experiments, we set $P=\\{0.0,0.1,\dots,1.0\\}$ and $Q=\\{0.0,0.5,1.0\\}$ in
the greedy peeling (Algorithm 2). The algorithm CGII+ILS by de Santiago and
Lamb [21] was also implemented in Python 3; however, to perform a fair
comparison, the algorithm was modified so that an initial
$\mathcal{S}^{\prime}$ of the RDP is set to be the singletons, consistent with
our proposed algorithm and the algorithm by Sato and Izunaga [57]. Note that
the choice of initial $\mathcal{S}^{\prime}$ is critical to the performance of
the algorithms; indeed, if we set it to an optimal solution to the modularity
density maximization problem, the algorithms are expected to terminate
shortly. The purpose of our experiments is to evaluate the performance of the
algorithms in terms of exactly or approximately solving the AP, irrespective
of the choice of initial $\mathcal{S}^{\prime}$. The parameter called “factor”
in CGII+ILS, which controls the degree of cluster perturbation, was set to
0.7, as recommended in de Santiago and Lamb [21]. As for the algorithm by Sato
and Izunaga [57], we used the implementation provided by the authors.

It should be noted that when implementing the above three algorithms, to
mitigate the effect of numerical errors, we have to modify them slightly so
that only the constraints that are violated by the incumbent solution
$\bm{\lambda^{*}}$ with some tolerance $\epsilon>0$ are added to the RDP (or
the counterpart as for the algorithm by Sato and Izunaga [57]). Without such a
modification, due to numerical errors, an optimal solution to the DP might be
decided to be infeasible for the DP, and the algorithm might not terminate. In
our experiments, the error parameter $\epsilon$ was set to $10^{-6}$.

For all algorithms, we used Gurobi Optimizer 9.5.0 as a mathematical
programming solver, with the default parameter settings. All the experiments
were conducted on a machine equipped with an Apple M1 chip and 16GB RAM.

###  5.2 Experiments and results

Table 2: Performance of our algorithm and the state-of-the-art methods.

| CGII+ILS [21] | Sato and Izunaga [57] | Ours (Algorithm 4) |
---|---|---|---|---
ID | time(s) | #constr. | time(s) | #constr. | time(s) | #constr. | Best $D(\mathcal{C})$ updated
1 | 5.9$\pm$0.3 | 479.4$\pm$53.8 | 0.4 | 62 | 0.3 | 780 | —
2 | 43.7$\pm$1.9 | 1056.0$\pm$45.5 | 0.7 | 80 | 2.4 | 3157 | —
3 | 3261.4$\pm$316.5 | 7567.6$\pm$603.2 | 16.4 | 185 | 42.8 | 21999 | —
4 | 13294.0$\pm$1991.1 | 18309.8$\pm$2036.0 | 57.9 | 205 | 119.7 | 59862 | —
5 |  $>\text{24h}$ (20.7943) | (70160) | 365.4 | 380 | 1886.3 | 184659 | —
6 |  $>\text{24h}$ (7.5892) | (2374) | $>\text{24h}$ | (11440) | 5681.1 | 121894 | $7.8250^{*}$
7 |  $>\text{24h}$ (34.0075) | (51181) | 70195.1 | 426 | 5860.8 | 359002 | $44.3879^{*}$
8 |  $>\text{24h}$ (32.7393) | (8068) | $>\text{24h}$ | (571) |  $>\text{24h}$ (47.2153) | (1265072) | —

First we compare the performance of our algorithm with the aforementioned
baseline methods. Time limit was set to 86,400s (24 hours). As explained
above, CGII+ILS has a stochastic behavior: For each instance, if the algorithm
does not terminate in the time limit at the first run, we do not run the
algorithm for the instance anymore. If it is not the case, we execute the
algorithm five times and report its statistics.

The results are shown in Table 2, where we report the computation time needed
and the number of constraints added to the RDP (or the counterpart as for the
algorithm by Sato and Izunaga [57]). For CGII+ILS, we report the average and
standard deviation of the computation time and the number of constraints over
the five trials. When an algorithm exceeds the time limit, we report the lower
bound on the optimal value of the DP (only for CGII+ILS and our proposed
algorithm) and the number of constraints at the time limit (in parentheses).
As the algorithm by Sato and Izunaga [57] solves the dual problems of
relaxations of the MP (using the set-packing relaxation) and tightens the
relaxation gradually, the lower bound kept in the algorithm is not necessarily
a lower bound on the optimal value of the DP; therefore, it is omitted in the
table. For each instance, the best result among the algorithms (in terms of
the computation time) is written in bold. We do not apply such a rule to the
number of constraints because the number of constraints is not necessarily a
reasonable indicator of the performance of algorithms; indeed, it may be a
good strategy to add constraints to the RDP actively to terminate quickly.

As can be seen, our proposed algorithm outperforms the baseline methods.
Indeed, for all instances except for the largest one (Jazz), only our
algorithm terminates in reasonable time (within two hours), and moreover, the
optimal solution to the MP derived by the optimal solution to the DP computed
by the algorithm has an integrality for each of those instances; thus, it
turns out that the algorithm obtains an optimal solution to the modularity
density maximization problem for those instances. In particular, the optimal
solutions to Adjnoun (ID 6) and Football (ID 7) are obtained for the first
time; the best previously-known modularity density values 7.651 and 44.340 are
improved to 7.8250 and 44.3879, respectively, with the optimality guarantee.
Although the algorithm by Sato and Izunaga [57] computes an optimal solution
faster for small instances, it has poorer performance for large instances.
Indeed, the algorithm does not even terminate for Adjnoun and consumes more
than 70,000s for Football. The performance of CGII+ILS is not comparable to
that of our proposed algorithm or the algorithm by Sato and Izunaga [57].

Table 3: Performance of the variants of our algorithm.

| Ours w/o greedy peeling | Ours w/o reformulation
---|---|---
ID | time(s) | #constr. | time(s) | #constr.
1 | 2.5$\pm$0.3 | 542.4$\pm$44.3 | 1.0 | 781
2 | 13.7$\pm$1.7 | 1359.0$\pm$51.8 | 4.7 | 3087
3 | 487.3$\pm$58.1 | 8660.8$\pm$713.7 | 1052.2 | 23605
4 | 1194.6$\pm$101.4 | 20146.6$\pm$1555.7 | 1230.0 | 62038
5 | 12491.4$\pm$1445.2 | 61678.6$\pm$5240.7 |  $>\text{10h}$ (21.5496) | (214623)
6 | 10185.9$\pm$847.4 | 46127.6$\pm$2092.3 |  $>\text{10h}$ (7.7835) | (118258)
7 |  $>\text{10h}$ (44.1272) | (42358) |  $>\text{10h}$ (44.3879) | (361468)
8 |  $>\text{10h}$ (31.3547) | (4052) |  $>\text{10h}$ (46.4529) | (850531)

Finally we investigate how much each of the proposed techniques, the greedy
peeling algorithm and the reformulation, contributes to the significant
performance of our algorithm. To this end, we implement and evaluate the
following two variants of our proposed algorithm, one replacing the greedy
peeling algorithm with the above heuristic method called the iterative local
search designed by de Santiago and Lamb [21] and the other replacing the
reformulation with AP-II by de Santiago and Lamb [21]. Time limit was set to
36,000s (10 hours). Note that the first variant has a stochastic behavior due
to the heuristic method by de Santiago and Lamb [21]; thus, we run the
algorithm and report the results in the same way as CGII+ILS.

The results are shown in Table 3. As can be seen, both algorithms have poorer
performance than our proposed algorithm but have better performance than
CGII+ILS (see also Table 2). Therefore, we see that both of the proposed
techniques do contribute to the significant performance of our algorithm.
However, we cannot conclude which of the techniques affects the performance
more significantly; indeed, only the algorithm using the reformulation
technique succeeds in computing an optimal solution to Polbooks (ID 5) and
Adjnonu (ID 6), but the algorithm using the greedy peeling obtains a better
lower bound for the larger instances Football (ID 7) and Jazz (ID 8).

##  6 Computational Complexity Analysis

In this section, we show the NP-hardness of a variant of the modularity
density maximization problem and the AP appearing in column generation.

###  6.1 NP-hardness of a variant of the modularity density maximization
problem

Here we study the following variant of the modularity density maximization
problem: Given an undirected graph $G=(V,E)$, we are asked to find a partition
$\mathcal{C}$ of $V$ that maximizes the modularity density $D(\mathcal{C})$
under the constraint that $|\mathcal{C}|\geq 2$. Therefore, this variant rules
out the trivial clustering $\mathcal{C}=\\{V\\}$. We call the variant the
modularity density maximization problem with at least two clusters. We prove
the following theorem:

######  Theorem 1.

The modularity density maximization problem with at least two clusters is NP-
hard even on regular graphs.

We first introduce the decision version of the problem:

######  Problem 1 (Modularity-Density($|\mathcal{C}|\geq 2$)).

Given an undirected graph $G=(V,E)$ and a real number $r\in\mathbb{R}$, does
there exist a partition $\mathcal{C}$ of $V$ with $|\mathcal{C}|\geq 2$ that
satisfies $D(\mathcal{C})\geq r$?

To prove Theorem 1, it suffices to show the NP-completeness of Modularity-
Density($|\mathcal{C}|\geq 2$) on regular graphs. For a given partition
$\mathcal{C}$ of $V$ with $|\mathcal{C}|\geq 2$ and a real number $r$, we can
check if $D(\mathcal{C})\geq r$ in polynomial time; thus, Modularity-
Density($|\mathcal{C}|\geq 2$) is in the class NP.

To show the NP-completeness, we construct a polynomial-time reduction from the
decision version of the maximum cut problem on $3$-regular graphs, which is
known to be NP-complete [3]. Let $G=(V,E)$ be an undirected graph. A
bipartition $\\{X,Y\\}$ of $V$ (i.e., a partition of $V$ with size two) is
called a cut of $G$. The value of cut $\\{X,Y\\}$ is defined as the number of
edges between $X$ and $Y$, i.e., $|\\{\\{u,v\\}\in E\mid u\in X,\ v\in Y\\}|$.
We denote by $\textsf{val}(\\{X,Y\\})$ the value of cut $\\{X,Y\\}$. Then we
can introduce the problem:

######  Problem 2 (3-Regular-Max-Cut).

Given a $3$-regular graph $G=(V,E)$ and a positive integer
$k\in\mathbb{Z}_{>0}$, does there exist a cut $\\{X,Y\\}$ of $G$ that
satisfies $\textsf{val}(\\{X,Y\\})\geq k$?

The following lemma guarantees the existence of the desired reduction:

######  Lemma 1.

There exists a polynomial-time reduction from 3-Regular-Max-Cut to Modularity-
Density($|\mathcal{C}|\geq 2$) on regular graphs.

Our reduction is based on an existing reduction from the maximum cut problem
to the (uniform) sparsest cut problem developed by Bonsma et al. [6]; however,
the difference from theirs is not totally trivial. Before starting the proof
of the lemma, we present a useful fact about the modularity density value for
regular graphs:

######  Lemma 2.

Let $G=(V,E)$ be a $d$-regular graph and $\mathcal{C}$ a partition of $V$.
Then the modularity density can be rewritten as follows:

$\displaystyle
D(\mathcal{C})=4\sum_{C\in\mathcal{C}}\frac{|E(C)|}{|C|}-d|\mathcal{C}|\quad\text{and}\quad
D(\mathcal{C})=d|\mathcal{C}|-2\sum_{C\in\mathcal{C}}\frac{|E(C,V\setminus
C)|}{|C|}.$

In particular, if $|\mathcal{C}|=2$, say $\mathcal{C}=\\{C,V\setminus C\\}$,
then the above second expression can be rewritten as

$\displaystyle D(\mathcal{C})=2d-2n\frac{|E(C,V\setminus C)|}{|C||V\setminus
C|}.$

###### Proof of Lemma 2.

We have

$\displaystyle
D(\mathcal{C})=\sum_{C\in\mathcal{C}}\frac{2|E(C)|-|E(C,V\setminus
C)|}{|C|}=\sum_{C\in\mathcal{C}}\frac{2|E(C)|-(d|C|-2|E(C)|)}{|C|}=4\sum_{C\in\mathcal{C}}\frac{|E(C)|}{|C|}-d|\mathcal{C}|$

and

$\displaystyle
D(\mathcal{C})=\sum_{C\in\mathcal{C}}\frac{2|E(C)|-|E(C,V\setminus
C)|}{|C|}=\sum_{C\in\mathcal{C}}\frac{d|C|-2|E(C,V\setminus
C)|}{|C|}=d|\mathcal{C}|-2\sum_{C\in\mathcal{C}}\frac{|E(C,V\setminus
C)|}{|C|}.$

The case of $|\mathcal{C}|=2$ is trivial. ∎

###### Proof of Lemma 1.

Let $(G,k)$ be an instance of 3-Regular-Max-Cut, where $G=(V,E)$ is a
$3$-regular graph and $k$ is a positive integer. Let $n=|V|$. Note that $n\geq
4$. From this instance, we construct an instance $(\overline{G^{*}},r^{*})$ of
Modularity-Density($|\mathcal{C}|\geq 2$) on regular graphs as follows: For
each $v\in V$, we construct two sets $I_{v}=\\{v_{0},v_{1},\dots,v_{M-1}\\}$
and
$I^{\prime}_{v}=\\{v^{\prime}_{0},v^{\prime}_{1},\dots,v^{\prime}_{M-1}\\}$,
each of which has size $M=n^{3}$. Note that for example, for a vertex named
$u$, we have $I_{u}=\\{u_{0},u_{1},\dots,u_{M-1}\\}$ and
$I^{\prime}_{u}=\\{u^{\prime}_{0},u^{\prime}_{1},\dots,u^{\prime}_{M-1}\\}$.
For each $e=\\{u,v\\}\in E$, we create two edges $\\{u_{0},v_{0}\\}$ and
$\\{u^{\prime}_{0},v^{\prime}_{0}\\}$. Hence, the subgraphs induced by
$\\{v_{0}\mid v\in V\\}$ and $\\{v^{\prime}_{0}\mid v\in V\\}$, respectively,
are isomorphic to $G$. For each $v\in V$, we connect two subsets $I_{v}$ and
$I^{\prime}_{v}$ as follows: add an edge connecting $v_{0}$ and each of
$\\{v^{\prime}_{j}\mid j=4,5,\dots,M-1\\}$; add an edge connecting $v_{i}$
($i=1,2,3$) and each of $\\{v^{\prime}_{j}\mid j=1,2,\dots,M-1\\}$; add an
edge connecting $v_{i}$ ($i=4,5,\dots,M-1$) and each of $\\{v^{\prime}_{j}\mid
j=0,1,\dots,M-1,\ j\neq i\\}$ (see Figure 1). We denote by
$G^{*}=(V^{*},E^{*})$ the graph constructed above. Note that $G^{*}$ is an
$(M-1)$-regular graph defined on $|V^{*}|=2Mn$ vertices. Let
$\overline{G^{*}}=(V^{*},\overline{E^{*}})$ be a complement graph of $G^{*}$.
In addition, set $r^{*}=2M-4-\frac{12}{M}+\frac{8k}{Mn}$. As
$\overline{G^{*}}$ is a $(2Mn-M)$-regular graph and $r^{*}$ is a positive real
number, $(\overline{G^{*}},r^{*})$ is a valid instance of Modularity-
Density($|\mathcal{C}|\geq 2$). Obviously we can construct
$(\overline{G^{*}},r^{*})$ from $(G,r)$ in polynomial time. Thus, to prove the
lemma, it suffices to show the following: $G$ has a cut $\\{X,Y\\}$ that
satisfies $\textsf{val}(\\{X,Y\\})\geq k$ if and only if $\overline{G^{*}}$
has a partition $\mathcal{C}$ of $V^{*}$ with $|\mathcal{C}|\geq 2$ that
satisfies $D(\mathcal{C})\geq r^{*}$.

Figure 1: The connection between $I_{v}$ and $I^{\prime}_{v}$ (on $G^{*}$).

Assume that $G$ has a cut $\\{X,Y\\}$ that satisfies
$\textsf{val}(\\{X,Y\\})\geq k$. Now consider the partition
$\mathcal{C}=\\{C,V^{*}\setminus C\\}$ of $V^{*}$ such that
$C=\left(\bigcup_{v\in X}I_{v}\right)\cup\left(\bigcup_{v\in
Y}I^{\prime}_{v}\right)$. Then using Lemma 2, we can evaluate the modularity
density value of $\mathcal{C}$ as follows:

$\displaystyle D(\mathcal{C})$
$\displaystyle=2(2Mn-M)-4Mn\frac{|\overline{E^{*}}(C,V^{*}\setminus
C)|}{(Mn)^{2}}$
$\displaystyle=2(2Mn-M)-4Mn\left(1-\frac{|E^{*}(C,V^{*}\setminus
C)|}{(Mn)^{2}}\right)$ $\displaystyle\geq
2(2Mn-M)-4Mn\left(1-\frac{n(M^{2}-M-3)+2k}{(Mn)^{2}}\right)$
$\displaystyle=2M-4-\frac{12}{M}+\frac{8k}{Mn}$ $\displaystyle=r^{*},$

where the inequality follows from the fact that for each $v\in V$, two subsets
$I_{v}$ and $I^{\prime}_{v}$ are connected by $M^{2}-M-3$ edges on $G^{*}$ and
the fact that $\textsf{val}(\\{X,Y\\})\geq k$. Therefore, we see that
$\mathcal{C}$ is a partition that we desired.

Assume that $\overline{G^{*}}$ has a partition $\mathcal{C}$ of $V^{*}$ with
$|\mathcal{C}|\geq 2$ that satisfies $D(\mathcal{C})\geq r^{*}$. First we show
that such a partition $\mathcal{C}$ has size exactly equal to two. Suppose for
contradiction that $|\mathcal{C}|\geq 3$. Then using Lemma 2, we have

$\displaystyle D(\mathcal{C})$
$\displaystyle=4\sum_{C\in\mathcal{C}}\frac{|\overline{E^{*}}(C)|}{|C|}-(2Mn-M)|\mathcal{C}|$
$\displaystyle\leq
4\sum_{C\in\mathcal{C}}\frac{|C|-1}{2}-(2Mn-M)|\mathcal{C}|$
$\displaystyle=4Mn-2|\mathcal{C}|-(2Mn-M)|\mathcal{C}|$
$\displaystyle\leq-2Mn+3M-6$ $\displaystyle<r^{*},$

where the last inequality follows from $n\geq 4$. This contradicts to
$D(\mathcal{C})\geq r^{*}$. Thus, we can represent $\mathcal{C}$ as
$\mathcal{C}=\\{C,V^{*}\setminus C\\}$.

Next we show that $|C|=|V^{*}\setminus C|=Mn$. Suppose for contradiction that
$|C|\leq Mn-1$. Then using Lemma 2, we have

$\displaystyle D(\mathcal{C})$
$\displaystyle=2(2Mn-M)-4Mn\frac{|\overline{E^{*}}(C,V^{*}\setminus
C)|}{|C||V^{*}\setminus C|}$
$\displaystyle=2(2Mn-M)-4Mn\left(1-\frac{|E^{*}(C,V^{*}\setminus
C)|}{|C||V^{*}\setminus C|}\right)$ $\displaystyle\leq
2(2Mn-M)-4Mn\left(1-\frac{(M-1)|C|}{|C|(2Mn-|C|)}\right)$
$\displaystyle=-2M+\frac{4Mn(M-1)}{2Mn-|C|}$
$\displaystyle\leq-2M+\frac{4Mn(M-1)}{Mn+1}$
$\displaystyle=2M-4-\frac{12}{M}-\frac{4M-12n-4-\frac{12}{M}}{Mn+1}$
$\displaystyle<2M-4-\frac{12}{M}+\frac{8k}{Mn}$ $\displaystyle=r^{*},$

where the first inequality follows from the fact that $G^{*}$ is
$(M-1)$-regular and the last inequality follows from $n\geq 4$. This
contradicts to $D(\mathcal{C})\geq r^{*}$. Hence, we can represent
$\mathcal{C}$ as $\mathcal{C}=\\{C,V^{*}\setminus C\\}$ in which
$|C|=|V^{*}\setminus C|=Mn$ holds.

Finally we show that for any $v\in V$, none of $I_{v}$ and $I^{\prime}_{v}$ is
divided by $\mathcal{C}$. Suppose for contradiction that for some $v\in V$, it
holds that $I_{v}\cap C\neq\emptyset$ and $I_{v}\cap(V\setminus
C)\neq\emptyset$. The following discussion also holds for $I^{\prime}_{v}$.
Recall that the number of edges between $I_{v}$ and $I^{\prime}_{v}$ on
$G^{*}$ is equal to $M^{2}-M-3$. An important fact is that among those
$M^{2}-M-3$ edges, at most $M^{2}-2M+2$ edges are included in
$E^{*}(C,V^{*}\setminus C)$, as $I_{v}$ is divided by $\mathcal{C}$. To see
this, let $a=|I_{v}\cap C|$ and $b=|I^{\prime}_{v}\cap C|$. Then we have

$\displaystyle|E^{*}((I_{v}\cup I^{\prime}_{v})\cap C,(I_{v}\cup
I^{\prime}_{v})\cap(V^{*}\setminus C))|$
$\displaystyle=a(M-b)+b(M-a)-|\overline{E^{*}}((I_{v}\cup I^{\prime}_{v})\cap
C,(I_{v}\cup I^{\prime}_{v})\cap(V^{*}\setminus C))|$
$\displaystyle=(a+b)M-2ab-|\overline{E^{*}}((I_{v}\cup I^{\prime}_{v})\cap
C,(I_{v}\cup I^{\prime}_{v})\cap(V^{*}\setminus C))|.$

Consider the case where $a+b=M$ holds. Then as $I_{v}$ is divided by
$\mathcal{C}$, we have $a\neq 0$ and $b\neq 0$, and so $ab\geq M-1$. Thus, we
have

$\displaystyle|E^{*}((I_{v}\cup I^{\prime}_{v})\cap C,(I_{v}\cup
I^{\prime}_{v})\cap(V^{*}\setminus C))|\leq M^{2}-2(M-1)=M^{2}-2M+2.$

Consider the other case where $a+b=\ell$ ($\ell=1,2,\dots,M-1$) holds. Then we
have

$\displaystyle|E^{*}((I_{v}\cup I^{\prime}_{v})\cap C,(I_{v}\cup
I^{\prime}_{v})\cap(V^{*}\setminus C))|$ $\displaystyle=\ell
M-2ab-|\overline{E^{*}}((I_{v}\cup I^{\prime}_{v})\cap C,(I_{v}\cup
I^{\prime}_{v})\cap(V^{*}\setminus C))|$ $\displaystyle\leq\begin{cases}\ell
M-0-\ell&(\text{if }b=0)\\\ \ell M-2(\ell-1)-0&(\text{otherwise})\end{cases}$
$\displaystyle\leq\ell(M-1)+1$ $\displaystyle\leq M^{2}-2M+2.$

Therefore, using the above fact, we can evaluate the modularity density value
of $\mathcal{C}$ as follows:

$\displaystyle D(\mathcal{C})$
$\displaystyle=2(2Mn-M)-4Mn\left(1-\frac{|E^{*}(C,V^{*}\setminus
C)|}{(Mn)^{2}}\right)$ $\displaystyle\leq
2(2Mn-M)-4Mn\left(1-\frac{3n+(n-1)(M^{2}-M-3)+M^{2}-2M+2}{(Mn)^{2}}\right)$
$\displaystyle=2M-4-\frac{4}{n}+\frac{20}{Mn}$
$\displaystyle<2M-4-\frac{12}{M}+\frac{8k}{Mn}$ $\displaystyle=r^{*},$

where the second inequality follows from $n\geq 4$. This contradicts to
$D(\mathcal{C})\geq r^{*}$. Moreover, by a similar argument (where
$M^{2}-2M+2$ in the above calculation is replaced by $0$), we can also show
that for any $v\in V$, it is impossible that both $I_{v}$ and $I^{\prime}_{v}$
are in the same cluster in $\mathcal{C}$. Therefore, we see that $\mathcal{C}$
can be represented as $\mathcal{C}=\\{C,V^{*}\setminus C\\}$ such that

$\displaystyle C=\left(\bigcup_{v\in X}I_{v}\right)\cup\left(\bigcup_{v\in
Y}I^{\prime}_{v}\right)$

for some cut $\\{X,Y\\}$ of $G$. Then we can evaluate the modularity density
value of $\mathcal{C}$ as follows:

$\displaystyle
D(\mathcal{C})=2M-4-\frac{12}{M}+\frac{8\cdot\textsf{val}(\\{X,Y\\})}{Mn}.$

As $D(\mathcal{C})\geq r^{*}$, we have $\textsf{val}(\\{X,Y\\})\geq k$.
Therefore, $\\{X,Y\\}$ is a cut that we desired, and we are done. ∎

As a by-product of our reduction, we can see that the uniform sparsest cut
problem is NP-hard even on regular graphs. This is a stronger result of the
NP-hardness of the uniform sparsest cut problem on general graphs, proved by
Bonsma et al. [6]. In the uniform sparsest cut problem, given an undirected
graph $G=(V,E)$, we are asked to find a cut $\\{C,V\setminus C\\}$ of $G$ that
minimizes $\frac{|E(C,V\setminus C)|}{|C||V\setminus C|}$. On the gadget
constructed above (i.e., the instance of (the decision version of) the
modularity density maximization problem with at least two clusters constructed
from an instance of (the decision version of) the maximum cut problem on
$3$-regular graphs), the modularity density maximization problem with at least
two clusters is essentially the same as the uniform sparsest cut problem. In
fact, any optimal solution $\mathcal{C}$ satisfies $|\mathcal{C}|=2$ and hence
by Lemma 2, the objective is equivalent to that of the uniform sparsest cut
problem.

######  Theorem 2.

The uniform sparsest cut problem is NP-hard even on regular graphs.

###  6.2 NP-hardness of the auxiliary problem

We investigate the computational complexity of the AP appearing in column
generation for the modularity density maximization problem. Specifically, we
study the following problem: Given an undirected graph $G=(V,E)$ and
$\bm{\lambda}=(\lambda_{v})_{v\in V}\in\mathbb{R}^{V}$, we are asked to find a
vertex subset $S\subseteq V$ that maximizes $\frac{2|E(S)|-|E(S,V\setminus
S)|}{|S|}-\sum_{v\in S}\lambda_{v}$. When solving the AP in a column
generation algorithm, $\bm{\lambda}$ is determined by an optimal solution
$\bm{\lambda^{*}}$ to the RDP. Here, to investigate the computational
complexity in general, we suppose that $\bm{\lambda}$ is also a part of the
input. Our goal is to prove the following theorem:

######  Theorem 3.

The AP (with $\bm{\lambda}$ input) is NP-hard even for $(n-4)$-regular graphs
and $\bm{\lambda}>\bm{0}$.

We again introduce the decision version of the problem:

######  Problem 3 (Auxiliary-Problem).

Given an undirected graph $G=(V,E)$, $\bm{\lambda}=(\lambda_{v})_{v\in
V}\in\mathbb{R}^{V}$, and a real number $r$, does there exist a vertex subset
$S\subseteq V$ that satisfies $\frac{2|E(S)|-|E(S,V\setminus
S)|}{|S|}-\sum_{v\in S}\lambda_{v}\geq r$?

To prove Theorem 3, it suffices to show the NP-completeness of Auxiliary-
Problem for $(n-4)$-regular graphs and $\bm{\lambda}>\bm{0}$. For a given
vertex subset $S\subseteq V$, $\bm{\lambda}\in\mathbb{R}^{V}$, and a real
number $r$, we can check if $\frac{2|E(S)|-|E(S,V\setminus
S)|}{|S|}-\sum_{v\in S}\lambda_{v}\geq r$ in polynomial time; thus, Auxiliary-
Problem is in the class NP.

To show the NP-completeness, we construct a polynomial-time reduction from the
decision version of the maximum clique problem on $(n-4)$-regular graphs,
which is known to be NP-hard [2, 29]. Let $G=(V,E)$ be an undirected graph. A
vertex subset $S\subseteq V$ is called a clique if every pair of vertices in
the subset has an edge in the graph. Then we can introduce the problem:

######  Problem 4 ($(n-4)$-Regular-$k$-Clique).

Given an $(n-4)$-regular graph $G=(V,E)$ and a positive integer
$k\in\mathbb{Z}_{>0}$, does there exist a clique $S\subseteq V$ of size $k$?

The following lemma guarantees the existence of the desired reduction:

######  Lemma 3.

There exists a polynomial-time reduction from $(n-4)$-Regular-$k$-Clique to
Auxiliary-Problem for $(n-4)$-regular graphs and $\bm{\lambda}>\bm{0}$.

Before starting the proof of the lemma, we borrow the following useful fact
from extremal graph theory:

######  Lemma 4 ((A corollary of) Turán’s theorem (Theorem 7.1.1 in Diestel
[23])).

Let $G=(V,E)$ be an undirected graph and $k$ a positive integer. If $G$ does
not contain a clique of size $k$, then
$|E|\leq\left(1-\frac{1}{k-1}\right)\frac{|V|^{2}}{2}$ holds.

###### Proof of Lemma 3.

Let $(G,k)$ be an instance of $(n-4)$-Regular-$k$-Clique, where $G=(V,E)$ is
an $(n-4)$-regular graph and $k$ is a positive integer. From this instance, we
construct an instance $(G^{*},\bm{\lambda^{*}},r^{*})$ of Auxiliary-Problem as
follows: We just use the same graph, i.e., $G^{*}=G$, and set
$\lambda^{*}_{v}=\frac{2(k-1)}{k}$ for every $v\in V$, and $r^{*}=-(n-4)$. To
prove the theorem, it suffices to show that $G$ has a clique $S\subseteq V$ of
size $k$ if and only if $G^{*}$ has a vertex subset $S\subseteq V$ that
satisfies $\frac{2|E(S)|-|E(S,V\setminus S)|}{|S|}-\sum_{v\in
S}\lambda^{*}_{v}\geq r^{*}$. Note that as $G$ is $(n-4)$-regular, we have
$\frac{2|E(S)|-|E(S,V\setminus S)|}{|S|}=\frac{4|E(S)|}{|S|}-(n-4)$.

Assume that $G$ has a clique $S\subseteq V$ of size $k$. Then we have

$\displaystyle\frac{4|E(S)|}{|S|}-(n-4)-\sum_{v\in
S}\lambda^{*}_{v}=2(|S|-1)-(n-4)-\frac{2(k-1)}{k}|S|=r^{*},$

meaning that $S$ is a vertex subset that we desired.

Assume that $G$ does not have a clique of size $k$. Then for any $S\subseteq
V$, $G[S]$ does not contain a clique of size $k$. Thus, by Lemma 4, we have
that for any $S\subseteq V$,

$\displaystyle|E(S)|\leq\left(1-\frac{1}{k-1}\right)\frac{|S|^{2}}{2}<\frac{(k-1)|S|^{2}}{2k},$

and hence

$\displaystyle\frac{|E(S)|}{|S|}<\frac{(k-1)|S|}{2k}.$

Therefore, we have

$\displaystyle\frac{4|E(S)|}{|S|}-(n-4)-\sum_{v\in
S}\lambda^{*}_{v}<\frac{2(k-1)}{k}|S|-(n-4)-\frac{2(k-1)}{k}|S|=-(n-4),$

which completes the proof. ∎

##  7 Conclusion

In this study, we have investigated modularity density maximization from both
algorithmic and computational complexity aspects. Specifically, we first
accelerated column generation for the modularity density maximization problem.
We first pointed out that the AP, i.e., the auxiliary problem appearing in
column generation, can be viewed as a dense subgraph discovery problem. Then
we employed a well-known strategy for dense subgraph discovery, called the
greedy peeling, for approximately solving the AP. Moreover, we reformulate the
AP to a sequence of $0$–$1$ linear programming problems, enabling us to
compute its optimal value more efficiently and to get more diverse columns.
Computational experiments using a variety of real-world networks demonstrate
the effectiveness of our proposed algorithm. In particular, our algorithm
obtains an optimal solution within only two hours to an instance with more
than 100 vertices that cannot be solved by the state-of-the-art methods in 24
hours. Finally, we showed the NP-hardness of a slight variant of the
modularity density maximization problem, where the output partition has to
have two or more clusters, as well as showing the NP-hardness of the AP.

There are several directions for future research. The most interesting one is
to design a more scalable exact or column generation algorithm for the
modularity density maximization problem. Although our proposed algorithm is
scalable than the state-of-the-art methods, it is still difficult to handle
large-scale graphs. Another direction is to apply the techniques similar to
that used in the present paper to other important optimization problems. Our
observation suggests that when designing a column generation algorithm for an
optimization problem based on the set partitioning formulation, we can see the
nature of dense subgraph discovery. Finally investigating the computational
complexity of the modularity density maximization problem rather than its
variants is a trivial future work. Is it possible to prove the NP-hardness of
the modularity density maximization problem by constructing a more
sophisticated polynomial-time reduction from some problem?

## Acknowledgments

The authors would like to thank Rafael de Santiago and Luís C. Lamb for
sharing the code of their algorithm, which makes it easier to implement their
algorithm in our experimental environment. The authors would also like to
thank Keisuke Sato and Yoichi Izunaga for sharing the code of their algorithm,
which is directly used in our experiments. A.M. is supported by Grant-in-Aid
for Early-Career Scientists (No. 19K20218). A.T. is supported by Grant-in-Aid
for Scientific Research (B) (No. 19H04069).

## References

  * [1] G. Agarwal and D. Kempe.  Modularity-maximizing graph communities via mathematical programming.  The European Physical Journal B, 66(3):409–418, 2008.
  * [2] C. C. Aggarwal, R. K. Ahuja, J. Hao, and J. B. Orlin.  Diagnosing infeasibilities in network flow problems.  Mathematical Programming, 81(3):263–280, 1998.
  * [3] P. Alimonti and V. Kann.  Hardness of approximating problems on cubic graphs.  In CIAC ’97: Proceedings of the 3rd Italian Conference on Algorithms and Complexity, pages 288–298, 1997.
  * [4] 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, 82:046112, 2010.
  * [5] 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, 2008(10):P10008, 2008.
  * [6] P. Bonsma, H. Broersma, V. Patel, and A. Pyatkin.  The complexity of finding uniform sparsest cuts in various graph classes.  Journal of Discrete Algorithms, 14:136–149, 2012.
  * [7] D. Boob, Y. Gao, R. Peng, S. Sawlani, C. E. Tsourakakis, D. Wang, and J. Wang.  Flowless: Extracting densest subgraphs without flow computations.  In WWW ’20: Proceedings of The Web Conference 2020, pages 573–583, 2020.
  * [8] F. Botta and C. I. del Genio.  Finding network communities using modularity density.  Journal of Statistical Mechanics: Theory and Experiment, 2016(12):123402, 2016.
  * [9] 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, 20(2):172–188, 2007.
  * [10] S. Cafieri, A. Costa, and P. Hansen.  Reformulation of a model for hierarchical divisive graph modularity maximization.  Annals of Operations Research, 222(1):213–226, 2014.
  * [11] S. Cafieri, P. Hansen, and L. Liberti.  Locally optimal heuristic for modularity maximization of networks.  Physical Review E, 83:056105, 2011.
  * [12] M. Charikar.  Greedy approximation algorithms for finding dense components in a graph.  In APPROX ’00: Proceedings of the 3rd International Workshop on Approximation Algorithms for Combinatorial Optimization, pages 84–95, 2000.
  * [13] C. Chekuri, K. Quanrud, and M. R. Torres.  Densest subgraph: Supermodularity, iterative peeling, and flow.  In SODA ’22: Proceedings of the 33rd Annual ACM–SIAM Symposium on Discrete Algorithms, pages 1531–1555, 2022.
  * [14] M. Chen, K. Kuzmin, and B. K. Szymanski.  Extension of modularity density for overlapping community structure.  In ASONAM ’14: Proceedings of the 2014 IEEE/ACM International Conference on Advances in Social Networks Analysis and Mining, pages 856–863, 2014.
  * [15] T. Chen, P. Singh, and K. E. Bassler.  Network community detection using modularity density measures.  Journal of Statistical Mechanics: Theory and Experiment, 2018(5):053406, 2018.
  * [16] A. Clauset, M. E. J. Newman, and C. Moore.  Finding community structure in very large networks.  Physical Review E, 70:066111, 2004.
  * [17] A. Costa.  Some remarks on modularity density, 2014.
  * [18] A. Costa.  MILP formulations for the modularity density maximization problem.  European Journal of Operational Research, 245(1):14–21, 2015.
  * [19] A. Costa, S. Kushnarev, L. Liberti, and Z. Sun.  Divisive heuristic for modularity density maximization.  Computers & Operations Research, 71:100–109, 2016.
  * [20] A. Costa, T. S. Ng, and L. X. Foo.  Complete mixed integer linear programming formulations for modularity density based clustering.  Discrete Optimization, 25:141–158, 2017.
  * [21] R. de Santiago and L. C. Lamb.  Exact computational solution of modularity density maximization by effective column generation.  Computers & Operations Research, 86:18–29, 2017.
  * [22] R. de Santiago and L. C. Lamb.  A ground truth contest between modularity maximization and modularity density maximization.  Artificial Intelligence Review, 53(6):4575–4599, 2020.
  * [23] R. Diestel.  Graph Theory, volume 173 of Graduate Texts in Mathematics.  Springer-Verlag Berlin Heidelberg, 5th edition, 2016.
  * [24] T. N. Dinh, X. Li, and M. T. Thai.  Network clustering via maximizing modularity: Approximation algorithms and theoretical limits.  In ICDM ’15: Proceedings of the 2015 IEEE International Conference on Data Mining, pages 101–110, 2015.
  * [25] T. N. Dinh and M. T. Thai.  Toward optimal community detection: From trees to general weighted networks.  Internet Mathematics, 11(3):181–200, 2015.
  * [26] S. Fortunato.  Community detection in graphs.  Physics Reports, 486(3):75–174, 2010.
  * [27] S. Fortunato and M. Barthélemy.  Resolution limit in community detection.  Proceedings of the National Academy of Sciences of the United States of America, 104(1):36–41, 2007.
  * [28] S. Fortunato and D. Hric.  Community detection in networks: A user guide.  Physics Reports, 659:1–44, 2016.
  * [29] M. R. Garey and D. S. Johnson.  Computers and Intractability: A Guide to the Theory of NP-Completeness.  W. H. Freeman & Co., NY, 1979.
  * [30] D. Gibson, R. Kumar, and A. Tomkins.  Discovering large dense subgraphs in massive graphs.  In VLDB ’05: Proceedings of the 31st International Conference on Very Large Data Bases, pages 721–732, 2005.
  * [31] A. Gionis and C. E. Tsourakakis.  Dense subgraph discovery: KDD 2015 Tutorial.  In KDD ’15: Proceedings of the 21st ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 2313–2314, 2015.
  * [32] M. Girvan and M. E. J. Newman.  Community structure in social and biological networks.  Proceedings of the National Academy of Sciences of the United States of America, 99(12):7821–7826, 2002.
  * [33] P. M. Gleiser and L. Danon.  Community structure in jazz.  Advances in Complex Systems, 06(04):565–573, 2003.
  * [34] A. V. Goldberg.  Finding a maximum density subgraph.  University of California Berkeley, 1984.
  * [35] B. H. Good, Y.-A. de Montjoye, and A. Clauset.  Performance of modularity maximization in practical contexts.  Physical Review E, 81:046106, 2010.
  * [36] R. Guimerà and L. A. Nunes Amaral.  Functional cartography of complex metabolic networks.  Nature, 433(7028):895–900, 2005.
  * [37] Y. Izunaga, T. Matsui, and Y. Yamamoto.  A doubly nonnegative relaxation for modularity density maximization.  Discrete Applied Mathematics, 275:69–78, 2020.
  * [38] Y. Izunaga and Y. Yamamoto.  A cutting plane algorithm for modularity maximization problem.  Journal of the Operations Research Society of Japan, 60(1):24–42, 2017.
  * [39] Y. Kawase, T. Matsui, and A. Miyauchi.  Additive approximation algorithms for modularity maximization.  Journal of Computer and System Sciences, 117:182–201, 2021.
  * [40] Y. Kawase and A. Miyauchi.  The densest subgraph problem with a convex/concave size function.  Algorithmica, 80(12):3461–3480, 2018.
  * [41] G. Kortsarz and D. Peleg.  Generating sparse 2-spanners.  Journal of Algorithms, 17(2):222–236, 1994.
  * [42] V. Krebs.  Orgnet.  http://www.orgnet.com/.
  * [43] Y. Kuroki, A. Miyauchi, J. Honda, and M. Sugiyama.  Online dense subgraph discovery via blurred-graph feedback.  In ICML ’20: Proceedings of the 37th International Conference on Machine Learning, pages 5522–5532, 2020.
  * [44] Z. Li, S. Zhang, R.-S. Wang, X.-S. Zhang, and L. Chen.  Quantitative function for community detection.  Physical Review E, 77:036109, 2008.
  * [45] 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, 54(4):396–405, 2003.
  * [46] C. P. Massen and J. P. K. Doye.  Identifying communities within energy landscapes.  Physical Review E, 71:046101, 2005.
  * [47] J. H. Michael.  Labor dispute reconciliation in a forest products manufacturing facility.  Forest Products Journal, 47(11/12):41, 1997.
  * [48] A. Miyauchi and N. Kakimura.  Finding a dense subgraph with sparse cut.  In CIKM ’18: Proceedings of the 27th ACM International Conference on Information and Knowledge Management, pages 547–556, 2018.
  * [49] A. Miyauchi and Y. Miyamoto.  Computing an upper bound of modularity.  The European Physical Journal B, 86(7):1–7, 2013.
  * [50] M. E. J. Newman.  Finding community structure in networks using the eigenvectors of matrices.  Physical Review E, 74:036104, 2006.
  * [51] M. E. J. Newman.  Modularity and community structure in networks.  Proceedings of the National Academy of Sciences of the United States of America, 103(23):8577–8582, 2006.
  * [52] M. E. J. Newman and M. Girvan.  Finding and evaluating community structure in networks.  Physical Review E, 69:026113, 2004.
  * [53] T. Richardson, P. J. Mucha, and M. A. Porter.  Spectral tripartitioning of networks.  Physical Review E, 80:036111, 2009.
  * [54] R. Santiago and L. C. Lamb.  Efficient modularity density heuristics for large graphs.  European Journal of Operational Research, 258(3):844–865, 2017\.
  * [55] R. Santiago and L. C. Lamb.  Efficient quantitative heuristics for graph clustering.  In GECCO 2017: Proceedings of the Genetic and Evolutionary Computation Conference Companion, page 117–118, 2017.
  * [56] K. Sato and N. Fukumura.  Real-time freight locomotive rescheduling and uncovered train detection during disruption.  European Journal of Operational Research, 221(3):636–648, 2012\.
  * [57] K. Sato and Y. Izunaga.  An enhanced MILP-based branch-and-price approach to modularity density maximization on graphs.  Computers & Operations Research, 106:236–245, 2019.
  * [58] R. Shang, W. Zhang, L. Jiao, R. Stolkin, and Y. Xue.  A community integration strategy based on an improved modularity density increment for large-scale networks.  Physica A: Statistical Mechanics and its Applications, 469:471–485, 2017.
  * [59] V. Spirin and L. A. Mirny.  Protein complexes and functional modules in molecular networks.  Proceedings of the National Academy of Sciences of the United States of America, 100(21):12123–12128, 2003.
  * [60] C. E. Tsourakakis.  The k-clique densest subgraph problem.  In WWW ’15: Proceedings of the 24th International Conference on World Wide Web, pages 1122–1132, 2015.
  * [61] C. E. Tsourakakis, F. Bonchi, A. Gionis, F. Gullo, and M. Tsiarli.  Denser than the densest subgraph: Extracting optimal quasi-cliques with quality guarantees.  In KDD ’13: Proceedings of the 19th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, pages 104–112, 2013.
  * [62] N. Veldt, A. R. Benson, and J. Kleinberg.  The generalized mean densest subgraph problem.  In KDD ’21: Proceedings of the 27th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, page 1604–1614, 2021.
  * [63] C. Yan and Z. Chang.  Modularized convex nonnegative matrix factorization for community detection in signed and unsigned networks.  Physica A: Statistical Mechanics and its Applications, 539:122904, 2020.
  * [64] W. W. Zachary.  An information flow model for conflict and fission in small groups.  Journal of Anthropological Research, 33(4):452–473, 1977.
  * [65] S. Zhang, X.-M. Ning, C. Ding, and X.-S. Zhang.  Determining modular organization of protein interaction networks by maximizing modularity density.  BMC Systems Biology, 4(2):1–12, 2010.