Department of Policy and Planning Sciences
Discussion Paper Series

No.1339

A doubly nonnegative relaxation for modularity density
maximization

by
Yoichi IZUNAGA, Tomomi MATSUI and Yoshitsugu YAMAMOTO

Feb 2016

UNIVERSITY OF TSUKUBA
Tsukuba, Ibaraki 305-8573
JAPANA DOUBLY NONNEGATIVE RELAXATIONFOR MODULARITY DENSITY MAXIMIZATIONYOICHI IZUNAGA, TOMOMI MATSUI, AND YOSHITSUGU YAMAMOTO
Department of Social Systems and Management Discussion Paper Series No.1339
Abstract. Modularity proposed by Newman and Girvan is the most commonly usedmeasure when the nodes of a graph are grouped into communities consisting of tightly
connected nodes. However, some authors pointed out drawbacks of the modularity, themain issue of which is resolution limit. Resolution limit refers to the sensitivity of themodularity to the total number of edges in the graph, which leaves small communities
not identified and hidden inside larger ones. To overcome this drawback, Li et al. haveproposed a new measure called modularity density. In this paper, introducing a variant of
the semidefinite programming called 0-1SDP, we show that the modularity density maximization can be modeled as 0-1SDP equivalently. Then we solve a relaxation problem of
0-1SDP to obtain an upper bound on the modularity density, and propose a lower bounding
algorithm based on the combination of spectral heuristics and dynamic programming.

1. Introduction
Network analysis has received a growing attention. One of the most important issues
in the network analysis is to find a meaningful structure, which often addresses identifying
or detecting community structure. Here communities are the sets of nodes consisting oftightly connected nodes, but loosely connected each other. Community detection is applied
to analyze the underlying relationship of diverse networks such as the social network, thebiological network, the world-wide-web, and VLSI network.
A variety of approaches to detect communities has been proposed by several researchers,
then a novel measure, called modularity, is proposed by Newman and Girvan [24]. The modularity was originally used as a stopping criterion of the hierarchical divisive algorithm [24],and then Newman [23] suggested the alternative approach of maximizing the modularity
directly since a high value of the modularity represents a good community structure. The
NP-hardness of the modularity maximization problem shown by Brandes et al. [6] turned
researchers’ attention to heuristic algorithms resulting in several efficient heuristics, while
exact algorithms have been proposed only in a few papers [2, 29, 3].
The modularity maximization became one of the central subjects of research, but someauthors pointed out its drawbacks, the main issue of which is resolution limit. Resolution
limit, brought up in Fortunato and Barthélemy [11], refers to the sensitivity of the modularity to the total number of edges in the graph, which leaves small communities not identified
and hidden inside larger ones. This narrows the application area of the modularity maximization since most of real-world networks may contain communities with different scales.
To overcome the resolution limit, there have been extensive studies so far [4, 22, 13, 28].
Date: January 18, 2016.
1

2

Y. IZUNAGA, T. MATSUI, AND Y. YAMAMOTO

Recently, Li et al. [17] have proposed a new measure for community detection, which iscalled modularity density, and the problem of maximizing the modularity density can bestraightforwardly formulated as a nonlinear binary programming.
As for the mathematical optimization approaches for the modularity density maximization, Costa [7] has presented some mixed integer linear programming formulations, MILP
for short, which enables an application of general-purpose solvers, e.g., CPLEX, Gurobi and
Xpress, to the problem. However, the number of communities must be fixed in advance,
and a difficult auxiliary problem need be solved in their formulations. More recently, ahierarchical divisive heuristics has been proposed by Costa et al. [8] to obtain a good lower
bound on the modularity density.
In this paper, for the modularity density maximization, we give a new formulation basedon a variant of the semidefinite programming called 0-1SDP. The advantage of this formulation is twofold: it does not require the number of communities be known, and its size is
independent of the number of edges of the graph in contrast to MILP formulations. In orderto obtain an upper bound on the modularity density, we propose to relax 0-1SDP to a semidefinite programming problem with non-negative constraints, and the resulting problem canbe solved in polynomial time. Moreover, we develop a method based on the combinationof spectral heuristics and dynamic programming to construct a feasible solution from thesolution obtained by the relaxation problem.
This paper is organized as follows. In Sections 2 and 3, giving the definition of themodularity density and a nonlinear binary programming formulation of the modularity
density maximization, we review some properties of the modularity density. In Section 4,
we present 0-1SDP formulation for the modularity density maximization and show theequivalence between both problems. In Section 5, we propose a method to solve a doubly
non-negative relaxation problem of 0-1SDP, and in Section 6 we explain a heuristic algorithm
which constructs a feasible solution by means of the solution of the relaxation problem. In
Section 7, we report the computational experiments. Finally, we give some conclusions andfurther research in Section 8.

2. Definitions and Notation
Let G = (V, E) be an undirected graph with the set V of n nodes and the set E of medges. We assume that V has at least two nodes. We say that Π = { C1 , C2 , . . . , Ck } is apartition of V if V = ∪kp=1 Cp , Cp ∩ Cq = ∅ for any distinct pair p and q, and Cp $= ∅ forany p. Each member Cp of a partition is called a community. We denote the set of edges
that have one end-node in C and the other end-node in C ! by E(C, C ! ). When C = C ! , weabbreviate E(C, C ! ) to E(C) for the sake of simplicity. Then modularity density, denotedby D(Π), for a partition Π is defined asD(Π) =! " 2|E(C)| −

C∈Π

#

! $
C ! ∈Π |E(C, C )||C|,

where | · | denotes the cardinality of the corresponding set. We refer to each term of theabove summation as the contribution of community to the modularity density.

DOUBLY NONNEGATIVE RELAXATION FOR MODULARITY DENSITY MAXIMIZATION

3

Modularity density maximization problem, MD for short, is to find a partition Π of V
that maximizes the modularity density D(Π), which is formulated as%#! $! " 2|E(C)| −%C ! ∈Π |E(C, C )|% maximize
%
|C|
MD %
C∈Π
%
% subject to Π is a partition of V .

A nonlinear binary programming formulation for MD has been proposed in Li et al. [17].
Although the optimal number of communities is a priori unknown similarly to the modularity
maximization problem, we suppose it is known for the time being, and denote it by t, and let
T be the index set { 1, 2, . . . , t }. Introducing a binary variable xip indicating whether node ibelongs to community Cp , we have the following nonlinear binary programming formulation:###%! " 2 i∈V j∈V aij xip xjp − i∈V di xip $%% maximize#%%i∈V xipp∈T%!
% subject to
xip = 1 (i ∈ V ),
%
NLP %
p∈T
%
!
%
xip ≥ 1 (p ∈ T ),
%
%
i∈V
%
%
xip ∈ { 0, 1 } (i ∈ V, p ∈ T ),where aij is the (i, j) element of the adjacency matrix A of graph G, and di is the degreeof node i ∈ V . The first set of constraints states that each node belongs to exactly one
community, and the second set of constraints imposes that each community should be anonempty subset of V . The objective function in this problem is the sum of fractional
functions with a quadratic numerator and a linear denominator. One of the widely usedsolution approaches for the problem of this kind is a parametric algorithm by Dinkelbach [9].
Another approach is a branch-and-bound algorithm [5, 16] in global optimization area.
3. Some Properties of Modularity Density
Now suppose that there exist several isolated nodes in a graph G. After removing theisolated nodes from G, we find a partition Π! that maximizes the modularity density on the
reduced graph of G. If the contribution of a community is non-negative for any communityof Π! , then Π! ∪{C̄} is an optimal partition on the original graph G, where C̄ consists of theisolated nodes once deleted. If there exist some communities with a negative contribution,then the contribution increases by adding the isolated nodes to these communities since thedenominator of the contribution increases. Therefore we have the following lemma.
Lemma 3.1 (Costa [7], Lemma 1.). The isolated nodes can be assigned to communities aposteriori.
Due to Lemma 3.1, we have only to consider graphs with no isolated nodes.
Proposition 3.2 (Costa [7], Proposition 1, Corollary 1, and Corollary 2.). Let Π! be a partition with maximum modularity density, then the size of each community is between 2 andn − 2(|Π! | − 1), i.e.,2 ≤ |C| ≤ n − 2(|Π! | − 1) for any C ∈ Π! .

4

Y. IZUNAGA, T. MATSUI, AND Y. YAMAMOTO

4. Formulations
In this section, we first present a reformulation of the modularity density maximization,
which is based on MILP formulation according to Costa [7]. Next, we show that modularity
density maximization can be equivalently formulated as 0-1SDP, a variant of the semidefinite
programming.
Hereafter, Sn , Sn+ and Nn denote the set of n × n symmetric matrices, the positive semidefinite cone, and the symmetric non-negative cone, i.e.,
Sn = { Y ∈ Rn×n | Y = Y $ },Sn+ = { Y ∈ Sn | u$ Y u ≥ 0 for all u ∈ Rn },Nn = { Y = (yij ) ∈ Sn | yij ≥ 0 for all i, j ∈ {1, 2, . . . , n} }.

For a given vector u, Diag(u) is the diagonal matrix with ui as the i-th diagonal element,
and vec(U ) is the vector obtained by stacking columns of a given matrix U .
4.1. MILP formulation. We can rewrite the objective function in the problem MD as
follows:& #'#! 4 {i,j}∈E xip xjp − i∈V di xip#,
i∈V xipp∈T

due to the definition of adjacency matrix A = (aij )ij∈V . The quadratic term xip xjp can belinearized by replacing it with a new variable yijp and adding the following Fartet inequalities [10]:yijp ≤ xip ,yijp ≤ xjp ,xip + xjp ≤ yijp + 1

for p ∈ T .

(4.1)

Note that the last inequality in (4.1) is redundant owing to the objective function of maximizing with respect to the variable yijp , hence can be omitted. Next, we introduce acontinuous variable αp defined as:##
4 {i,j}∈E yijp − i∈V di xip#
αp =.
i∈V xip
This equality constraint can be relaxed to##
4 {i,j}∈E yijp − i∈V di xip#
αp ≤,i∈V xip(4.2)

without overlooking an optimal solution due to the objective function. Since the denominator in (4.2) is positive, we can rewrite it as!
!
!
αp
xip ≤ 4yijp −
di xip .
i∈V{i,j}∈E

i∈V

Finally, to linearize the product αp xip , we then introduce a continuous variable γip to replace
αp xip and make use of the following McCormick inequalities [19]:Lα xip ≤ γip ≤ Uα xip

for i ∈ V , p ∈ T ,αp − Uα (1 − xip ) ≤ γip ≤ αp − Lα (1 − xip ) for i ∈ V , p ∈ T ,DOUBLY NONNEGATIVE RELAXATION FOR MODULARITY DENSITY MAXIMIZATION

5

where Lα and Uα are lower and upper bounds of αp , respectively. From the above discussion,
MILP formulation is given as

MILP%% maximize%%%% subject to
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%
%

!

αp

p∈T
!

xip = 1
p∈T !
2≤

i∈V

(i ∈ V ),xip ≤ n − 2 (t − 1)y!
yijp ≤ xjp !
ijp ≤ xip ,!
γip ≤ 4yijp −
di xipi∈V{i,j}∈Ei∈V

Lα xip ≤ γip ≤ Uα xip
αp − Uα (1 − xip ) ≤ γip ≤ αp − Lα (1 − xip )xip ∈ {0, 1}
yijp ∈ R
Lα ≤ αp ≤ Uα
γip ∈ R(p ∈ T ),({i, j} ∈ E, p ∈ T ),(p ∈ T ),(i ∈ V, p ∈ T ),(i ∈ V, p ∈ T ),(i ∈ V, p ∈ T ),({i, j} ∈ E, p ∈ T ),(p ∈ T ),(i ∈ V, p ∈ T ).

From the inequality (4.2) and Proposition 3.1, a valid lower bound on the variable αp is
attained when the corresponding community consists of only two nodes with the largest
degree, thus Lα is given as Lα = −(dmax1 + dmax2 )/2 where dmax1 and dmax2 are the largest
and the second largest degrees, respectively. On the other hand, in order to obtain theupper bound on αp , we need to solve the following auxiliary problem:##%
4 {i,j}∈E xi xj − i∈V di xi%% maximize#%%i∈V xi!
AP %% subject to 2 ≤xi ≤ n − 2(t − 1),
%
i∈V
%
%
xi ∈ {0, 1} (i ∈ V ).

The problem AP is a problem of maximizing a nonlinear objective function with binary variables, thus difficult to optimize globally. Using a nonlinear programming solver SCIP [1],Costa solved the continuous relaxation problem of AP. In our experiment which is done forthe purpose of comparing the solutions obtained by MILP formulation and 0-1SDP formulation introduced in Subsection 4.2, we solve the problem AP to optimality by Dinkelbach’s
parametric algorithm.
4.2. 0-1SDP formulation. Now we consider the following problem:%% maximize Tr((2A − D)Z)%% subject to Zen = en ,0-1SDP %%Z 2 = Z,%%Z ∈ Nn ,where D is a diagonal matrix whose i-th diagonal entry is the degree of node i, i.e., D =Diag(d1 , d2 , . . . , dn ) ∈ Rn×n , and en is the n-dimensional vector of 1’s. We call the problem
0-1SDP because Peng and Xia stated “we call it 0-1SDP owing to the similarity of theconstraint Z 2 = Z to the classical 0-1 requirement in integer programming” in [25].
Henceforth we will show the equivalence between the problems MD and 0-1SDP.

6

Y. IZUNAGA, T. MATSUI, AND Y. YAMAMOTO

Lemma 4.1. For any feasible solution Z = (zij ) of 0-1SDP, there is a matrix X which
satisfies
Xe# = en ,(4.3)

X $ en ≥ e# ,X ∈ {0, 1}

(4.4)
(4.5)

n×#

,

$

−1

Z = X(X X)

$

(4.6)

X ,

for some positive integer #.
Proof. Let Z be a feasible solution of 0-1SDP. Then it clearly satisfies the positive semidefiniteness due to the idempotence constraint. Then there exists an index i1 ∈ V which
satisfies zi1 i1 = max{ zij | i, j ∈ V }, which is positive owing to the non-negativity of Z and
the constraint Zen = en . Let us define the index set I1 = { j ∈ V | zi1 j > 0 }, then wereadily obtain the following equalities!
!
(zi1 j )2 = zi1 i1 andzi1 j = 1,
j∈I1

j∈I1

due to the constraints Z 2 = Z, Zen = en and the symmetry of Z. Since zi1 i1 is positive,
the first equality reduces to! " zi j $1
zi1 j = 1.
zi1 i1j∈I1

By using the second equality, this yields! " zi j $zi1 j =

1j∈I1

which is rewritten as

zi1 i1! " zi j
1j∈I1

zi1 i1!

zi1 j ,j∈I1$
− 1 zi1 j = 0.

From the non-negativity of zij and the maximality of zi1 i1 , we have (zi1 j /zi1 i1 − 1) = 0,
which implies zi1 j = zi1 i1 for any j ∈ I1 . Then we havezi1 i1 = zi1 j =

1
|I1 |

for any j ∈ I1 .

(4.7)

For an index j ∈ I1 , we consider the (i1 , j) element, denoted by Zi21 j , of the matrix Z 2 ,which is given as


!
!
!
Zi21 j =zi1 k zkj =zi1 k zkj = zi1 i1 
zkj  .
k∈V

k∈I1

k∈I1

Note that the last equality is due to (4.7). The above equality, Z 2 = Z and (4.7) yield!
zi1 i1 = zi1 i1 
zkj  ,k∈I1

DOUBLY NONNEGATIVE RELAXATION FOR MODULARITY DENSITY MAXIMIZATION

which is reduced to
1=

!

7

zkj .

k∈I1

This implies that zjk = 1/|I1 | for any j, k ∈ I1 owing to the maximality of zi1 i1 and theconstraint Zen = en . Denote the sub-matrix (zij )i,j∈I1 by ZI1 . By permuting the rows andcolumns of Z simultaneously, we obtain
,
ZI1 O
$
P ZP =O ZĪ1
where P is an appropriate permutation matrix, and Ī1 = V \ I1 . Then it is clear that the
sub-matrix ZĪ1 satisfies the following:
ZĪ1 e|Ī1 | = e|Ī1 | , ZĪ21 = ZĪ1 and ZĪ1 ∈ N|Ī1 | .

Thus repeating the process described above, we can convert Z to a block diagonal matrix
as follows:


ZI1


ZI2


$
P ZP = 
,
..


.
ZI!

where each block diagonal element ZIp is the |Ip |×|Ip | matrix whose elements are all 1/|Ip |.
Now, we construct a matrix X = (xip ) such that4
1 (if i ∈ Ip ),xip =0 (otherwise),then X clearly satisfies (4.3), (4.4) and (4.5), and we can confirm Z = X(X $ X)−1 X $by a simple calculation. Note that the inverse of X $ X exists since the matrix X $ X is anonsingular diagonal matrix whose diagonal entry is the number of 1’s of each column of X
by (4.3), (4.4) and (4.5).
!

We say that the matrix X = (xip ) which satisfies the conditions (4.3), (4.4) and (4.5) isan incidence matrix of a partition Π = { C1 , C2 , . . . , Ck }, that is, each member Cp of Π isrepresented as Cp = { i ∈ V | xip = 1 }.
Theorem 4.2. The problem 0-1SDP is equivalent to the problem MD.
Proof. Let Z be an optimal solution of 0-1SDP, then we obtain an incidence matrix X
representing a partition Π which satisfies Z = X(X $ X)−1 X $ due to Lemma 4.1, and let us
denote the optimal value of an optimization problem by ω. Since the objective function in
the problem MD is written as Tr((2A − D)Z) by means of the matrix which satisfies (4.6),
we have the following inequality
Tr((2A − D)Z) = Tr((2A − D)X(X $ X)−1 X $ )#! $! " 2|E(C)| −C ! ∈Π |E(C, C )|
=
|C|
C∈Π

≤ ω(MD).

(4.8)

8

Y. IZUNAGA, T. MATSUI, AND Y. YAMAMOTO

Next, we show the optimality of the partition Π obtained from Z for the problem MD.
For an optimal solution Π̂ of MD, let X̂ be an incidence matrix corresponding to Π̂, and let
Ẑ = X̂(X̂ $ X̂)−1 X̂ $ . Clearly Ẑ satisfies
Ẑen = Ẑ X̂et = X̂et = enfrom (4.3) and (4.6). Moreover, it is symmetric, non-negative and idempotent, i.e., Ẑ ∈ Nn
and Ẑ 2 = Ẑ due to (4.6), hence Ẑ is feasible to 0-1SDP. Then we have#
! $! " 2|E(C)| −C ! ∈Π̂ |E(C, C )|ω(MD) =|C|C∈Π̂

= Tr((2A − D)X̂(X̂ $ X̂)−1 X̂ $ )= Tr((2A − D)Ẑ)≤ Tr((2A − D)Z) = ω(0-1SDP).

(4.9)

The last inequality is due to the optimality of Z for the problem 0-1SDP. By the inequalities (4.8) and (4.9), we conclude that ω(MD) = ω(0-1SDP), which implies the optimality ofthe partition Π obtained from Z for the problem MD.
!
Note that the size of 0-1SDP depends on neither the number of edges nor the numberof communities. Moreover, we need not solve the auxiliary problem unlike the case of
MILP formulation. These features make 0-1SDP more attractive than MILP formulation.
Although the objective function in 0-1SDP is linear with respect to the matrix Z, theidempotence constraint makes the problem difficult. We will discuss how to deal with thisdifficult part in the next section.
5. Doubly Nonnegative Relaxation
As stated in Subsection 4.2, what makes 0-1SDP difficult to solve is the idempotence
constraint, which imposes the condition that each eigenvalue of Z, denoted by λi , is either
0 or 1. Then it would be a natural strategy to relax the constraint to a more tractable
constraint. The first choice is to relax the binary constraint to 0 ≤ λi ≤ 1, which isexpressed as Z ∈ Sn+ and I − Z ∈ Sn+ , where I is the identity matrix. The latter constraint
I − Z ∈ Sn+ which represents the upper bound constraint of λi is redundant owing to thepresence of Z ∈ Nn and Zen = en , hence can be omitted. Then we obtain the following
relaxation problem over the doubly non-negative cone Nn ∩ Sn+ :%% maximize Tr((2A − D)Z)%
DNN %% subject to Zen = en ,%Z ∈ Nn ∩ Sn+ .

The optimization problems over a symmetric cone are solved efficiently, e.g., linear programming, second-order cone programming, and semidefinite programming problems. Indeed, the primal-dual-interior-point method solves the problems in polynomial time. On
the other hand, since the doubly non-negative cone is not symmetric, we cannot directly
apply the primal-dual-interior-point method to solve the problem DNN. Representing the
doubly non-negative cone as a symmetric cone embedded in a higher dimension, we could

DOUBLY NONNEGATIVE RELAXATION FOR MODULARITY DENSITY MAXIMIZATION

9apply the primal-dual-interior-point method to the embedded problem which is describedas follows:%% maximize Tr((2A − D)Z)%% subject to Zen = en ,
%
,
%
Z
O
+
%
∈ Sn+n2.
%
O Diag(vec(Z))
Although the above problem can be solved in polynomial time theoretically, it is a quite large
optimization problem over the positive semidefinite cone and computationally expensive tosolve in practice. Nevertheless it is worthwhile to solve the problem DNN due to the factthat the doubly non-negative relaxation often provides significantly tight bound for somecombinatorial optimization problems.
Now, we introduce valid inequalities for 0-1SDP in order to strengthen the bound obtainedby the relaxation problem. From the proof of Lemma 4.1, any feasible solution Z of 0-1SDP
can be transformed to a block diagonal matrix. It is easy to see that the maximum value inthe i-th row of Z is located on the (i, i) element for each i ∈ V , hence we have the following
result.
Lemma 5.1. The following inequalities are valid for 0-1SDP.
zii ≥ zij

for i, j ∈ V .

We denote the problem DNN with the above valid inequalities added by DNN.
6. Heuristics based on Dynamic Programming
In this section, we will develop an algorithm to construct a feasible solution for MD
from the solution obtained by solving the relaxation problems DNN and DNN presented at
the end of Section 5. The algorithm we propose is based on the combination of spectral
heuristics and dynamic programming.
6.1. Permutation based on spectrum. From the proof of Lemma 4.1, we have seen thatan optimal solution of 0-1SDP forms a matrix with block diagonal structure by applying anappropriate simultaneous permutation of the rows and columns, and each block corresponds
to a community. Unless otherwise stated, we refer to the simultaneous permutation of therows and columns simply as a permutation. The optimal solution of the problem DNN or
DNN is not necessarily transformed to a block diagonal matrix by any permutation since
we relaxed some constraints. The solution however may provide a clue as to possibly a goodsolution of the original problem MD. Thus, it would be helpful to transform the solution
matrix to a matrix which is close to a block diagonal one. To this end, we exploit the
spectrum of the optimal solution.
Let λ1 , λ2 , . . . , λn ∈ R be the eigenvalues of an optimal solution Z ∗ for the relaxation
problem and u1 , u2 , . . . , un ∈ Rn be their corresponding eigenvectors. We assume thatthe eigenvalues are sorted in the non-increasing order, that is, 1 = λ1 ≥ λ2 ≥ · · · ≥λn ≥ 0. Focusing on the eigenvector u2 = (u21 , u22 , . . . , u2n )$ corresponding to the second
largest eigenvalue, we permute the rows and columns of the matrix Z ∗ consistent with thenon-decreasing order of values of components of u2 . Take a benchmark instance Karate
for example, we show an optimal solution Z ∗ of DNN and the matrix obtained in themanner described above in Figures 1 (a) and (b). The color-code in the figures indicates the

10

Y. IZUNAGA, T. MATSUI, AND Y. YAMAMOTO

magnitude of elements: elements whose value is greater than or equal to 1/n are marked
with gray, while other elements (with value of less than 1/n) are marked with white. From(a) Original matrix(b) Permuted matrix

Figure 1: Comparison with two matrices

these figures, we observe that the permuted matrix is considerably close to a block diagonal
matrix. However, there is no theoretical validity of using the eigenvector corresponding to
the second largest eigenvalue. Here still remains room for further research.
6.2. Dynamic programming. Next, we discuss how to construct a feasible solution of
MD from the permuted matrix. Let V̄ be a sequence of nodes consistent with the nondecreasing order of components of u2 . For the sake of notational simplicity, we renumber
the nodes in V and denote the sequence V̄ by (1, 2, . . . , n). Now, we try to find a partition
with maximum modularity density of V̄ under the constraint that each member consistsof consecutive indices of V̄ . For this problem, we propose an algorithm using the dynamic
programming. We define q(k, l) by###
2 li=k lj=k aij − li=k diq(k, l) =l − (k − 1)for k and l of V̄ with k ≤ l. The value q(k, l) represents the contribution of the community(k, . . . , l) of V̄ to the modularity density when it is selected as a member of the partition.
For any index s of V̄ , let µ(s) be the maximum modularity density that is achieved bypartitioning the sequence (1, . . . , s) into several consecutive subsequences. If we defineµ(0) = 0 for notational convenience, then µ(s) satisfies the recursive equationµ(s) = max{ µ(h) + q(h + 1, s) | h ∈ {0, 1, . . . , s − 1} }.

(6.1)

Owing to (6.1), starting from µ(1) = q(1, 1), we first obtainµ(2) = max{ q(1, 2), µ(1) + q(2, 2) },then obtain µ(3) by means of µ(1) and µ(2), and so on. Our algorithm is given as follows.
Algorithm DP
Step 1 : Solve the relaxation problem to obtain an optimal solution Z ∗ .

DOUBLY NONNEGATIVE RELAXATION FOR MODULARITY DENSITY MAXIMIZATION

11

Step 2 : Find the eigenvector u2 corresponding to the second largest eigenvalue of Z ∗ .
Let V̄ = (1, 2, . . . , n) be a sequence of nodes obtained by rearranging them in nondecreasing order of corresponding components in u2 .
Step 3 : Set µ(0) := 0.
for s = 1 to n do
Compute µ(s) according to (6.1).
end-do7. Computational Experiments
To evaluate the lower and upper bounds obtained by our algorithm, we conducted thecomputational experiments on a computer with Intel Core i7, 3.70 GHz processor and 32.0 GB
of memory. Using SeDuMi 1.2 as an SDP solver, we implemented the algorithm in MATLAB 2010.
In the experiments, we used seven instances; Michael’s strike dataset [20], Zachary’s
karate dataset [30], Gil-Mendieta and Schmidt’s Mexico dataset [12], Michael and Massey’s
sawmill dataset [21], Lusseau’s dolphins dataset [18], Hugo’s Les Misérables dataset [14], andKrebs’ books dataset [15]. Details of each instance are listed in Table 1, where the columns“ LB! ” and “ t! ” represent the known best lower bound and the corresponding numberof communities, respectively. Since the first four instances in the table were solved to
optimality in Costa [7], we list the optimal values and the optimal numbers of communities,
for the remaining instances, the lower bounds and the number of communities reported in
Costa et al. [8] since the instances were not solved so far to our knowledge.

Table 1: Instances
ID
name
n
1
Strike24
2Karate34
3
Mexico35
4
Sawmill365Dolphins
62
6 Les Misérables 77
7
Books
105

m
LB! t!
38 8.8611 4
78 7.8451 3
117 8.7180 3
62 8.6233 4
159 12.1252 5
254 24.5339 9
441 21.9652 7

Table 2 shows the computational results of the algorithm described in Subsection 6.2,
where the columns “UB”, “LB”, “gap” and “time” represent the optimal value of the relaxation problem, the lower bound obtained by the algorithm DP, the duality gap defined by100(UB − LB)/LB, and the computation time in seconds, respectively. For each instance,
we observed that the predominant portion of the computation time was spent for solving
the relaxation problem, and the remaining parts of the algorithm require a fraction of time,
specifically less than one second.

12

Y. IZUNAGA, T. MATSUI, AND Y. YAMAMOTO

Table 2: Computational results of our algorithm
ID
1
2
3
4
5
6
7

DNN
UB
LB
gap(%) time(sec.)
9.5808 8.8611
8.122
1.05
8.9548 7.8424 14.184
5.83
10.3151 8.5580 20.532
7.64
10.5048 7.0486 49.034
7.75
15.0218 9.8286 52.838
316.6128.0957 22.2680 35.279
658.2826.5387 20.2470 31.075
4626.11

DNN
UB
LB
gap(%) time(sec.)
9.3049 8.8611
5.008
3.54
8.4141 7.8451
7.253
36.04
9.9570 8.5227 16.829
43.48
10.0311 7.3587 36.316
54.21
14.3552 11.4610 25.253
1681.8127.4276 23.3416 17.505
7018.0324.7749 20.3150 21.953 60437.45

From Table 2, we observe that the upper bounds UB provided by DNN are tighter thanthose provided by DNN for all instances, which indicates the effectiveness of the valid inequalities in Lemma 5.1. Moreover, we also confirm that the lower bounds obtained for DNN
are equal to or larger than those obtained for DNN for all instances except Mexico (ID=3).
However, solving the problem DNN requires a rather long computation time compared withsolving DNN as the instance size grows. To be specific, for the instance Books (ID=7), DNN
takes approximately 4,600 seconds, whereas DNN takes over 60,000 seconds, approximately
16 hours.
Table 3: Computational results of the branch-and-bound algorithm for MILP formulation
ID
1
2
3
4
5
6
7

UB
8.8611
7.84518.71808.623313.8466
61.630254.6234

MILP
LB
gap(%) time(sec.)
8.86110
0.50
7.84510
0.74
8.71800
7.84
8.62330
6.10
12.1252 14.196 86400.00
24.5339 151.204 86400.00
21.6803 151.949 86400.00

Next, we solve the problem MILP by the branch-and-bound algorithm in Gurobi 6.0.0
to compare the quality of the solutions obtained by our algorithm for 0-1SDP formulation.
The computational results are given in Table 3. In our experiments, we impose a time limitof 86,400 seconds, 24 hours, on the computation time of the branch-and-bound algorithm.
In Table 3, we see that the first four instances are solved to optimality within a shortcomputation time, while the remaining instances cannot be solved within the time limit.
Especially, for the instances Les Misérables (ID=6) and Books (ID=7), a quite large duality
gap remains. Figures 2 (a) and (b) show the upper and lower bounds vs. the elapsed time.
From the figures, we observe the following behavior: (i) the branch-and-bound algorithm
gives a good lower bound in early stage, and (ii) the improvement of upper bound rarely

DOUBLY NONNEGATIVE RELAXATION FOR MODULARITY DENSITY MAXIMIZATION

80

13

60

50
60
40
40
30

obj.val.

obj.val.

20

0

20

10

0
-20
-10
-40
-20
UB
LB

-60
0

10000

20000

30000

40000
50000
elapsed time (sec.)

60000

70000

80000

UB
LB

-30
90000

0

10000

(a) Instance Les Misérables

20000

30000

40000
50000
elapsed time (sec.)

60000

70000

80000

90000

(b) Instance Books

Figure 2: The upper and lower bounds vs. elapsed time in the branch-and-boundoccurs throughout the computation. Owing to the latter of these, the duality gap still
remains large even though a good feasible solution has been found. This suggests thatderiving a tight upper bound enables us to estimate the accuracy of an incumbent solution
obtained by the branch-and-bound algorithm more precisely.
8. Conclusion
In this paper, we presented 0-1SDP formulation originally introduced by Peng and Xia [25]for minimum sum-of-squares clustering problem, and showed that the equivalence between
the problem 0-1SDP and the modularity density maximization. The problem 0-1SDP has
the big advantage that its size is independent of the number of both edges and communities.
Then, we proposed to solve a doubly non-negative relaxation of the problem 0-1SDP inorder to obtain an upper bound on the modularity density. In addition, we developed alower bounding algorithm based on the combination of spectral heuristics and dynamic
programming.
The relaxation problem for our formulation was numerically compared to MILP formulation, and the results showed that the upper bounds provided by our formulation are muchtighter than those provided by MILP formulation.
We observed that the solution matrix showed a block-diagonal-like structure when its
rows and columns are rearranged simultaneously in accordance with the magnitude of thecomponents of the eigenvector corresponding to the second largest eigenvalue. Theoretical
study should be carried out about whether this procedure functions effectively, and why
if it does. The alternative to form a block-diagonal-like matrix is to develop a heuristics
based on the numerical linear algebraic computation such as the algorithms of Sargent andWesterberg [26], and Tarjan [27].
References[1] T. Achterberg. SCIP: Solving constraint integer programs. Mathematical Programming Computation,
1(1):1–41, 2009.
[2] G. Agarwal and D. Kempe. Modularity-maximizing graph communities via mathematical programming.
The European Physical Journal B, 66(3):409–418, 2008.
[3] 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(4):046112, 2010.

14

Y. IZUNAGA, T. MATSUI, AND Y. YAMAMOTO[4] A. Arenas, A. Fernández, and S. Gómez. Analysis of the structure of complex networks at different
resolution levels. New Journal of Physics, 10(5):053039, 2008.
[5] H.P. Benson. Global optimization algorithm for the nonlinear sum of ratios problem. Journal of Optimization Theory and Applications, 112(1):1–29, 2002.
[6] 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, 2008.
[7] A. Costa. MILP formulations for the modularity density maximization problem. European Journal ofOperational Research, 245(1):14–21, 2015.
[8] A. Costa, S. Kushnarev, L. Liberti, and Z. Sun. Divisive heuristic for modularity density maximization.
Optimization Online, 2015. http://www.optimization-online.org/DB HTML/2015/09/5116.html.
[9] W. Dinkelbach. On nonlinear fractional programming. Management Science, 13(7):492–498, 1967.
[10] R. Fortet. Applications de lalgebre de boole en recherche opérationelle. Revue Française de RechercheOpérationelle, 4(14):17–26, 1960.
[11] S. Fortunato and M. Barthélemy. Resolution limit in community detection. Proceedings of the National
Academy of Sciences, 104(1):36–41, 2007.
[12] J. Gil-Mendieta and S. Schmidt. The political network in mexico. Social Networks, 18(4):355–381, 1996.
[13] C. Granell, S. Gomez, and A. Arenas. Hierarchical multiresolution method to overcome the resolution
limit in complex networks. International Journal of Bifurcation and Chaos, 22(7):1250171, 2012.
[14] D.E. Knuth. The Stanford GraphBase: A platform for combinatorial computing, volume 37. AddisonWesley Reading, 1993.
[15] V Krebs. http://www.orgnet.com/(last visited october 13, 2015.).
[16] T. Kuno. A branch-and-bound algorithm for maximizing the sum of several linear ratios. Journal of
Global Optimization, 22(1-4):155–174, 2002.
[17] Z. Li, S. Zhang, R.S. Wang, X.S. Zhang, and L. Chen. Quantitative function for community detection.
Physical Review E, 77(3):036109, 2008.
[18] 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.
[19] G.P. McCormick. Computability of global solutions to factorable nonconvex programs: Part I convex
underestimating problems. Mathematical Programming, 10(1):147–175, 1976.
[20] J.H. Michael. Labor dispute reconciliation in a forest products manufacturing facility. Forest Products
Journal, 47(11/12):41–45, 1997.
[21] J.H. Michael and J.G. Massey. Modeling the communication network in a sawmill. Forest Products
Journal, 47(9):25–30, 1997.
[22] S. Muff, F. Rao, and A. Caflisch. Local modularity measure for network clusterizations. Physical Review
E, 72(5):056107, 2005.
[23] M.E. Newman. Fast algorithm for detecting community structure in networks. Physical Review E,
69(6):066133, 2004.
[24] M.E. Newman and M. Girvan. Finding and evaluating community structure in networks. Physical Review
E, 69(2):026113, 2004.
[25] J. Peng and Y. Xia. A new theoretical framework for k-means-type clustering. In W. Chu and T.Y.
Lin, editors, Foundations and Advances in Data Mining, volume 180 of Studies in Fuzziness and Soft
Computing, pages 79–96. Springer, 2005.
[26] R.W.H. Sargent and A.W. Westerberg. Speed-up in chemical engineering design. Transactions of theInstitution of Chemical Engineers, 42:190–197, 1964.
[27] R. Tarjan. Depth-first search and linear graph algorithms. SIAM Journal on Computing, 1(2):146–160,
1972.
[28] V.A. Traag, P. Van Dooren, and Y. Nesterov. Narrow scope for resolution-limit-free community detection. Physical Review E, 84(1):016114, 2011.
[29] G. Xu, S. Tsoka, and L.G. Papageorgiou. Finding community structures in complex networks using
mixed integer optimisation. The European Physical Journal B, 60(2):231–239, 2007.
[30] W.W. Zachary. An information flow model for conflict and fission in small groups. Journal of Anthropological Research, 33(4):452–473, 1977.

DOUBLY NONNEGATIVE RELAXATION FOR MODULARITY DENSITY MAXIMIZATION

15(Y. Izunaga) Graduate School of Systems and Information Engineering, University of Tsukuba,
Tsukuba, Ibaraki 305-8573, Japan
E-mail address: s1130131@sk.tsukuba.ac.jp(T. Matsui) Department of Social Engineering, Graduate School of Decision Science andTechnology, Tokyo Institute of Technology, Meguro-ku, Tokyo 152-8550, Japan
E-mail address: matsui.t.af@m.titech.ac.jp(Y. Yamamoto) Faculty of Engineering, Information and Systems, University of Tsukuba,
Tsukuba, Ibaraki 305-8573, Japan
E-mail address: yamamoto@sk.tsukuba.ac.jp