

Department of Social Systems and Management
Discussion Paper Series

No.1309

A cutting plane algorithm for Modularity maximization
with heuristics for separation problem

byYoichi IZUNAGA and Yoshitsugu YAMAMOTO

July 2013

UNIVERSITY OF TSUKUBA
Tsukuba, Ibaraki 305-8573
JAPANA CUTTING PLANE ALGORITHM
FOR MODULARITY MAXIMIZATIONWITH HEURISTICS FOR SEPARATION PROBLEM
YOICHI IZUNAGA AND YOSHITSUGU YAMAMOTO
Abstract. The 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. We formulate the modularity maximization problem as a set partitioning problem, andpropose an algorithm for the problem based on the linear programming relaxation. We solve
the dual of the linear programming relaxation by using a cutting plane method. To mediate theslow convergence that cutting plane methods usually suffer, we propose a method for finding
and simultaneously adding multiple cutting planes which may complement well each other.

1. Introduction
As social network services grow, clustering on graphs has been attracting more attention,and since Newman and Girvan [11] proposed the modularity as a graph clustering measure,modularity maximization problem became one of the central subjects of research. The NP hardness of the modularity maximization problem shown by Brandes et al. [3] turned researchers’
attention to heuristic algorithms, which resulted in several efficient heuristic algorithms such as
the linear programming with rounding procedure [1] and the variable neighborhood search [10].
On the other hand, among exact algorithms two approaches should be mentioned. The first
approach is based on the formulation of the problem as a clique partitioning problem proposedby Grötschel and Wakabayashi [9]. In this formulation a binary variable corresponding to eachpair of nodes represents whether the two nodes belong to the same cluster. Then the numbers ofvariables and constraints amount to O(n2 ) and O(n3 ), respectively, both of which grow rapidly
with the number n of nodes. Based on this formulation Aloise et al. [2] solved instances up to
115 nodes by using the cutting plane algorithm proposed by Grötschel and Wakabayashi [9].
The second formulation is the set partitioning problem. Since this formulation has to take
into account all nonempty subsets of the node set, it has O(2n ) variables. We can hardly secure
the computational resource to hold the problem when n is large, say more than 200. Column
generation technique is a common trick to deal with such problems, but the auxiliary problem
of determining the entering column ends up as a non-convex quadratic programming in binary
variables, which is hard to solve exactly. Moreover, the column generation is known to suffer
slow convergence. To overcome this defect, du Merle et al. [8] have proposed an acceleration
method, called stabilized column generation method. Some computational results are reported
also in Aloise et al. [2].
In this paper, based on the set partitioning formulation, we propose a cutting plane algorithm.
Finding a “best” cutting plane in a sense is NP -hard, so that we propose to use a heuristic
algorithm to find a cutting plane to add. We also propose a method to find multiple cutting
planes to accelerate the convergence, and report on some computational experiments.
Date: July 20, 2013.
Key words and phrases. Modularity maximization, Set partitioning, Clique partitioning, Cutting planes, Integer programming.
The second author is grateful to the Okawa Foundation for Information and Telecommunication for theirfinancial support.
1

2

Y. IZUNAGA AND Y. YAMAMOTOThis paper is organized as follows. In Section 2, we give the definition of modularity andformulate the modularity maximization problem as a set partitioning problem. In Section 3, we
introduce the relaxation problem and its linear programming dual problem. In Section 4, after
reviewing cutting plane algorithms, we propose two versions of the cutting plane algorithm withheuristics for generating a cutting plane. In Section 5, we explain a local search algorithm which
is used in our proposed algorithm. In Section 6, we report the computational experiments ofthe proposed algorithm. Finally in Section 7 we give some conclusions.
2. Modularity maximization problem
Let G = (V, E) be an undirected graph with the set V = {1, 2, . . . , n} of n nodes and the set!
E of undirected m edges. We say that Π = {C1 , C2 , . . . , Ck } is a partition of V if V = kp=1 Cp ,Cp ∩ Cq = ∅ for any distinct p and q in {1, 2, . . . , k}, and Cp #= ∅ for any p ∈ {1, 2, . . . , k}. Eachmember Cp of a partition is called a cluster or a community. We denote the set of edges thathave one end-node in Cp and the other end-node in Cq by E(Cp , Cq ). When Cp = Cq we write
E(Cp , Cq ) simply as E(Cp ).
Modularity, denoted by µ(Π), of a partition Π is defined as

%
'2 
&k
k
"
|E(C
)|
+
|E(C
,
C
)|
p
p
q
q=1
 |E(Cp )| −,µ(Π) =m2m
p=1

where | · | denotes the cardinality of the corresponding set. See Newman and Girvan [11] andBrandes et al. [3]. For i, j ∈ V let eij be*
1 when {i, j} ∈ Eeij =0 otherwisei.e., the (i, j) element of the adjacency matrix of graph G, and di be the degree of node i, i.e.,
di = |{ j ∈ V | {i, j} ∈ E }|, and π(i) be the index of community which node i belongs to, i.e.,π(i) = p means i ∈ Cp . Then µ(Π) is rewritten as
+
,
di dj
1 ""
µ(Π) =eij −δ(π(i), π(j)),
2m2mi∈V j∈V

where δ is the Kronecker delta, i.e.,*
1 when p = q
δ(p, q) =0 when p #= q.

Modularity Maximization problem, MM for short, is the problem of finding a partition of V
that maximizes the modularity µ(Π). Denoting (eij − di dj /2m) by wij , then the problem isformulated as
1 ""wij δ(π(i), π(j))- maximize
2m
MM : -i∈V j∈V- subject to Π is a partition of V .

Let P denote the family of all nonempty subsets of V . Note that P is composed of 2n − 1
subsets of V . Introducing a binary variable zC for each C ∈ P, a partition Π is represented bythe (2n − 1)-dimensional binary vector z = (zC )C∈P defined as*
1 when C ∈ ΠzC =0 otherwise.

CUTTING PLANE FOR MODULARITY MAXIMIZATION

3This enables us to formulate problem MM as an integer programming problem. For each i ∈ V
and C ∈ P let aiC be defined by*
1 when i ∈ C
aiC =0 otherwise.
The column aC = (a1C , . . . , anC )" is the incidence vector of community C, i.e., C = { i ∈ V |aiC = 1 }. For each C ∈ P let fC be
1 ""
fC =
(2.1)
wij ,
2m
i∈C j∈C

which is rewritten as

1 ""wij aiC ajC
2mi∈V j∈V
"
1 " "
wij aiC ajC +wii aiC  .
=m=i∈V j∈V : i<ji∈VThe last equality is due to the symmetry of wij . The constant fC represents the contributionof community C to the objective function µ(Π) when community C is selected as a member ofthe partition Π. Thus MM is formulated as the integer programming P :"1 "fC zC +
wii
- maximize
m
C∈P
i∈V
"
P : -- subject to
aiC zC = 1(∀i ∈ V )C∈Pz ∈ {0, 1}(∀C ∈ P).
C

Since the first set of constraints states that the communities adopted form a partition of V , thisproblem is called a set partitioning problem. Due to its huge number of variables this problem
easily becomes computationally intractable as the number of nodes grows.
3. Relaxation problem and its dual problem

Not only the number of variables but also their integrality makes problem P a highly intractable problem. Then it would be a natural and clever strategy to consider relaxation problems for the useful information about the solution of P . See, for example [4, 5, 12, 13]. The first
choice to consider would be the Linear Programming relaxation, LP relaxation for short, where
the binary constraint zC ∈ {0, 1} is replaced by 0 ≤ zC ≤ 1. It is given as"
1 "fC zC +
wii
- maximize
m
C∈P
i∈V
"
RP : -- subject to
aiC zC = 1(∀i ∈ V )C∈Pz ≥0(∀C ∈ P).
CThe upper bound constraints zC ≤ 1 are redundant owing to the first set of constraints, hence
omitted. The linear programming dual problem of RP is given as the following RD:
"
1 "
λi +
wii
- minimize
m
i∈V
i∈V
"
RD : -- subject to
aiC λi ≥ fC(∀C ∈ P)i∈V
λ ∈R(∀i ∈ V ).
i

4

Y. IZUNAGA AND Y. YAMAMOTO

Note that fC = 0 when C is a singleton from (2.1). Thus n of inequality constraints imply
the non-negativity of λi ’s. Therefore the problem RD is equivalent to the problem with nonnegativity constraints of variables:
"
1 "
λi +
wii
- minimize
m
i∈V
i∈V
"
RD : -- subject to
aiC λi ≥ fC(∀C ∈ P)i∈V
λi ≥ 0(∀i ∈ V ).

Now let us denote the feasible region and the optimal value of an optimization problem,
say Q, by F(Q) and ω(Q), respectively. Since RP is a relaxation problem of P , it holds that
F(P ) ⊆ F(RP ), hence ω(P ) ≤ ω(RP ). Applying the linear programming duality theorem tothe primal dual pair RP and RD, we see ω(P ) ≤ ω(RD). Namely, solving RD we will obtain an
upper bound of ω(P ). The optimal solution of RP often provides a clue as to possibly a good
feasible solution of P with aid of the information collected in solving its dual RD. Although
RD has only n variables, its exponentially large number of constraints makes it intractable.
4. Cutting plane algorithm with heuristics
The constraints of problem RD far outnumber the variables, hence most of them should not bebinding at an optimal solution. The cutting plane algorithm is one of commonly used methods
for LP problems of this kind. We will give a brief review of cutting plane algorithms, and thenpropose an algorithm for RD.
4.1. Prototype of cutting plane algorithm. The key idea of the cutting plane algorithm isto deal with a small subfamily C of P, and instead of RD to solve the following problem withfewer constraints:
"
1 "
λi +
wii
- minimize
m
i∈V
i∈V
"
RD(C) : -- subject to
aiC λi ≥ fC(∀C ∈ C)i∈V
λi ≥ 0(∀i ∈ V ).
&Let λ∗ (C) denote an optimal solution of RD(C). Since the constraints i∈V aiC λi ≥ fC for
C ∈ P \ C are not considered, it is not necessarily a feasible solution of RD. To check thefeasibility of λ∗ (C), we define a measure of violation γC (C) of the constraint corresponding to C
as
"
(4.1)
γC (C) =aiC λ∗i (C) − fC .
i∈V

Note that γC (C) ≥ 0 for all C ∈ C. When γC (C) ≥ 0 for all C ∈ P \ C, λ∗ (C) is a feasible solutionof problem RD, hence an optimal solution of problem RD. WhenγC (C) < 0
holds for some C ∈ P \ C, adding this C to C can lead to an improvement of the optimal valueof problem RD(C), i.e., ω(RD(C ∪ {C})) > ω(RD(C)). Substituting (2.1) for fC of (4.1) yields"
1 " "γC (C) =aiC λ∗i (C) −wij aiC ajC ,mi∈Vi∈V j∈V ; i<jCUTTING PLANE FOR MODULARITY MAXIMIZATION

5hence the problem of minimizing γC (C) over P is formulated as the problem AP (C) with aquadratic objective function in binary variables:"
1 " "λ∗i (C)yi −wij yi yj- minimizem
AP (C) : -i∈Vi∈V j∈V ; i<j- subject to y ∈ {0, 1}(∀i ∈ V ).
i

An optimal solution y ∗ of this problem provides the incidence vector of a community thatminimizes γC (C) over P. Since y = 0 is a feasible solution of this problem, the optimal value isnon-positive. Finding y ∗ with negative optimal value, we have only to add the constraint"yi∗ λi ≥ f ∗i∈V

to RD(C), where

f∗ =

1 ""
wij yi∗ yj∗ .
2mi∈V j∈V

From the above discussion, a prototype of the cutting plane algorithm can be given as follows.
Prototype of the Cutting Plane Algorithm
Step 0Let C be an initial family of nonempty subsets of V .
Step 1
Solve RD(C) to obtain an optimal solution λ∗ (C) and the optimal value ω(RD(C)).
Step 2
Solve AP (C) and set y ∗ be an optimal solution.
Step 3
If ω(AP (C)) ≥ 0 thenSet C ! ← C and ω ! ← ω(RD(C)). Output C ! and ω ! , and terminate.
else
Set C ← { i ∈ V | yi∗ = 1 } and increment C ← C ∪ {C}. Return to Step 1.
end if

The initial subfamily C could be empty, but a clever choice may enhance the efficiency of thealgorithm. In our experiments, as reported in Section 6, we collected all singletons of V to makethe initial family C. When the algorithm terminates, we have solved RD, hence also RP , whichusually admits a fractional optimal solution. The final family C ! however would yield a set ofprimal variables that are likely to be positive at an optimal solution of problem P . Then we
propose to solve the following problem P (C ! ) of variables zC with C ∈ C ! .
"1 "fC z C +wii- maximize
m
!
C∈C
i∈V
"
P (C ! ) : -- subject to
aiC zC = 1(∀i ∈ V )C∈C !
zC ∈ {0, 1}(∀C ∈ C ! ).

This problem is expected to have much fewer variables than problem P does, so that it could be
solved within a reasonable time by an IP solver, e.g., CPLEX, Gurobi, Xpress. Lacking variables
zC with C not in C ! , P (C ! ) provides a lower bound of ω(P ). Then
ω(P (C ! )) ≤ ω(P ) ≤ ω(RD(C ! )).

6

Y. IZUNAGA AND Y. YAMAMOTO

The value ω(RD(C ! )) − ω(P (C ! )) provides an upper bound of the difference between ω(P (C ! ))and ω(P ), hence the quality of the solution of P (C ! ) given by an IP solver.
4.2. Proposed algorithm. Since problem AP (C) to be solved in Step 2 of the prototype cutting
plane algorithm is an NP -hard non-convex quadratic programming with binary variables, we
propose to apply a heuristic algorithm. We put off the description of the heuristic algorithm
until Section 5, and we first discuss how the stopping criterion should be modified accordingly.
Let α(AP (C)) denote the objective function value of a solution provided by a heuristic algorithm for AP (C) in Step 2. The non-negativity of α(AP (C)) does not imply the non-negativity ofω(AP (C)), hence we cannot claim that we have reached a solution of RD even if α(AP (C)) ≥ 0
holds. Therefore, we propose to stop the algorithm when α(AP (C)) continues to be non-negative
for a predetermined number of successive iterations. Our proposed algorithm is described asfollows, where kmax is the predetermined number of iterations. Note that C is incremented by asingle set C determined by y ∗ found in Step 2. We will call this algorithm the Single-CuttingPlane-at-a-Time Algorithm, SCP for short.
Single-Cutting-Plane-at-a-Time Algorithm (SCP )Step 0Let C be an initial family of nonempty subsets of V .
Determine a natural number kmax , and set k ← 0.
Step 1
Solve RD(C) to obtain an optimal solution λ∗ (C) and the optimal value ω(RD(C)).
Step 2
Apply the heuristic algorithm to AP (C) (see Section 5 for the algorithm).
Set y ∗ and α(AP (C)) be a solution obtained and its objective function value, respectively.
Step 3if α(AP (C)) ≥ 0 thenk ← k + 1.
else
Set C ← { i ∈ V | yi∗ = 1 }, increment C ← C ∪ {C}, and set k ← 0. Go to Step 1.
end ifStep 4
if k > kmax thenSet C ! ← C and ω ! ← ω(RD(C)). Output C ! and ω ! , and terminate.
else
Go to Step 2.
end ifCarrying out some preliminary experiments by SCP, we frequently observed that the optimalvalue ω(RD(C)) stays constant for many iterations even when C is repeatedly incremented. We
show in Table 1 how slowly ω(RD(C)) increases as the algorithm SCP progresses. Concerning“Dolphins” and “Jazz” the optimal value ω(RD(C)) did not change at all after 500th iteration.
Here “Dolphins”, “Football” and “Jazz” are the benchmark instances available from the site:
http://www.cc.gatech.edu/dimacs10/archive/clustering.shtml
The size as well as the known optimal value ω(P ) of each instance is given in Table 2.
The slow convergence we observed may arise from a particular structure of RD(C) that all
coefficients of the objective function are one and all coefficients of the constraints are either
zero or one. This makes the contour of the objective function and some face of F(RD(C)) beparallel, and the whole face be optimal. As a consequence, the optimal value ω(RD(C)) stays

CUTTING PLANE FOR MODULARITY MAXIMIZATION

7

Table 1. Plateau situation of SCP
iteration
1
500
1000
2000
3000
4000
5000
6000

ω(RD(C))
Dolphins Football
Jazz-0.0213 -0.0087 -0.0070
0.42410.4430 0.3276
0.42410.4470 0.3276
0.42410.4549 0.3276
0.42410.4587 0.3276
0.42410.4605 0.3276
0.42410.4605 0.3276
0.42410.4609 0.3276

Table 2. Solved instancesnamenm ω(P )
Dolphins 62 159 0.5285
Football 115 613 0.6045
Jazz
198 2742 0.4448
constant although lots of cutting planes are added. To cut off such a face entirely, we propose
to simultaneously add multiple cutting planes which may complement well each other. The first
cutting plane is the same as the one defined by y ∗ and f ∗ obtained from problem AP (C) in (4.1).
We then fix the variables yi to zero for all i with yi∗ = 1, and consider AP (C). More precisely,
we let V (1) = { i ∈ V | yi∗ = 1 } and approximately solve the problem
"
1 " "
∗
- minimizeλ(C)y−wij yi yjiimi∈Vi∈V j∈V ; i<jAP (C, V (1) ) : - subject to yi ∈ {0, 1}(∀i∈ V \ V (1) )yi = 0(∀i ∈ V (1) )to obtain y (1) and f (1) , i.e., the second cutting plane. In a general step, with V (h) = { i ∈ V |(l)yi = 1 for some l < h } we approximately solve
"
1 " "
∗
- minimizeλ(C)y−wij yi yjiimi∈Vi∈V j∈V ; i<j(h)AP (C, V ) : - subject to yi ∈ {0, 1}(∀i ∈ V \ V (h) )yi = 0(∀i ∈ V (h) ),where y (0) = y ∗ and V (0) = ∅. We call these problems Restricted Auxiliary problems. As longas α(AP (C, V (h) )) is negative, we keep on generating cutting planes. When α(AP (C, V (h) ))becomes non-negative, we add all the cutting planes obtained so far to C. The Step 3 of thealgorithm SCP should be modified as follows. We will call the algorithm with this modification
the Multiple-Cutting-Planes-at-a-Time Algorithm, MCP for short.
Step 3 of the Multiple-Cutting-Planes-at-a-Time Algorithm (MCP )Step 3if α(AP (C)) ≥ 0 thenk ← k + 1.
else

8

Y. IZUNAGA AND Y. YAMAMOTO

Generate cutting planes until α(AP (C, V (h) )) becomes non-negative.
Add all the cutting planes generated to C. Set k ← 0 and go to Step 1.
end if

5. Heuristic algorithm based on the noising method
Now we discuss the heuristic algorithm for AP (C) used in Step 2, which is based on theNoising method originally proposed by Charon and Hudry [6, 7].
We start with the explanation of the conventional local search. Given a feasible solution y of
AP (C), we define the neighborhood N (y) of y as
N (y) = { y $ | ,y $ − y,1 ≤ 1 }.

Then each vector in N (y) is obtained by replacing a component, say yq of y by its complement
1 − yq . We denote the thus obtained vector by y q . The variation, denoted by v(y, y q ), of theobjective function of AP (C) when y changes to y q is given as"
1v(y, y q ) = (1 − 2yq ) λ∗q (C) −wiq yi  .
m
i∈V \{q}The above variation (5) is computed for each solution in N (y). If they are all non-negative,then the local search stops and return y as a local optimal solution. Otherwise, we move y to
the solution which minimizes v(y, y q ) over N (y). The algorithm is described as follows.

Conventional Local Search
Step 0Let y be an initial feasible solution of AP (C).
Step 1for q = 1 to n do
Compute v(y, y q ).
Step 2v ∗ ← min{v(y, y q )} and y ∗ ← argmin{v(y, y q )}.
if v ∗ ≥ 0 thenOutput y, and terminate.
else
Set y ← y ∗ , and go to Step 1.
end if

The above naive local search is notorious for that it is liable to get trapped in a local optimal
solution. To overcome this drawback Charon and Hudry [6, 7] proposed the Noising method,which is based on the local search but differs in that the objective function is perturbed byadding a random value, which they call noise.
First, we give an initial feasible solution y, and find the best solution y ∗ in N (y) by the
local search. Then we calculate the following perturbed variation ṽ(y, y ∗ ) by adding a numberρ randomly chosen from the interval [−r, r]:ṽ(y, y ∗ ) = v(y, y ∗ ) + ρ.

CUTTING PLANE FOR MODULARITY MAXIMIZATION9

If ṽ(y, y ∗ ) is negative, we move y to y ∗ . The interval [−r, r] is reduced to [−r + d, r − d] after
a predetermined number of iterations, and the Noising method terminates when the interval
shrinks to {0} or vanishes. Noising method is described as follows.
Noising Method
Step 0Let y be an initial feasible solution.
Set r ← 100, d ← 1, FixedIter ← 10, and t ← 0.
Step 1
Set t ← t + 1 and call Local Search to obtain a local optimal solution y ∗ in N (y).
Step 2
Draw the noise ρ randomly from the interval [−r, r] and compute ṽ(y, y ∗ ).
if ṽ(y, y ∗ ) < 0 then set y ← y ∗ .
Step 3
if t ≡ 0 (mod FixedIter) then set r ← r − d and go to Step 1.
Step 4
if r = 0 thenOutput y and terminate.
else
Go to Step 1.
end if

6. Computational experiments
We report the computational experiment with MCP. The experiment was performed on a
PC with an Intel Core2 Duo, 3.06GHz processor and 4.0GB of memory. We implemented thealgorithm in Java, and used CPLEX 12.3 as the LP solver, and solved the three benchmark
instances introduced in Subsection 4.2.
Since the result may change due to the random noise ρ used in the algorithm, we executed
the algorithm five times for each instance. We set C initially to the family of all singletons, i.e.,
C = {{1}, {2}, . . . , {n}}, and the parameter kmax to 30. The statistics collected are given in
Table 3, and Table 4 shows the results of 15 executions.
Table 3. Statistics|C ! |
cardinality of the final family of subsets C !
!
ω(RD(C )) optimal value of RD(C ! ) obtained at the end of the algorithm
ω(P (C ! )) optimal value of P (C ! )
.
/
gap
time

relative gap defined by gap =computation time in seconds

ω(P )−ω(P (C ! ))
ω(P )

× 100

From Table 4, we see that MCP solved the instances of Dolphins and Football to optimality,
although it does not provide a proof of optimality. At four out of five executions MCP failed tosolve Jass, but the gap was less than 0.5%. Aloise et al. [2] reported that the stabilized column
generation method solved Dolphins in approximately seven seconds on a PC with Intel Pentium,
3.20GHz processor and 3.0GB of memory, and CPLEX 10.110 as the LP solver. However they
could not solve Football after the computation of more than 100,000 seconds.

10

Y. IZUNAGA AND Y. YAMAMOTO

Table 4. Computational results of MCP
instanceexec.
1
2
Dolphins
3
4
5
1
2
Football
3
4
5
1
2
Jazz
3
4
5

|C ! | ω(RD(C ! )) ω(P (C ! )) gap (%) time (s)797
0.52850.52850.000
11982
0.52850.52850.000
23676
0.52850.5285
0.000
9
1,081
0.52850.52850.000
23998
0.52850.5285
0.000
16
5,617
0.60450.6045
0.000
326
5,046
0.60450.6045
0.000
291
5,248
0.60450.6045
0.000
312
4,982
0.60450.6045
0.000
265
5,029
0.60450.60450.000
294
26,8040.44390.44360.269
25,002
26,8900.44480.44480.000
19,78026,8720.44450.44450.067
19,18826,5090.44450.44450.067
21,238
27,1030.44340.44290.427
25,001

We also observe that the number of generated constraints |C ! | is much smaller than that of
the original problem. Take Dolphins with 62 nodes for instance, the generated constraints areless than 1/1015 of the original constraints totalling 4.6 × 1018 .
Figures 1, 2 and 3 provide a comparison of SCP with MCP. We stopped SCP when the numberof cutting planes generated amounts to the same number of those generated by MCP. The figures
show ω(RD(C)) vs. the number of iterations. For both algorithms ω(RD(C)) rapidly increases
at an early stage, and then increases slowly or stays constant as the algorithm progresses.
0.6

0.5

0.4

obj.val

0.3

0.2

0.1

0
SCP
MCP
!(P)

-0.1
0

100

200

300

400
iteration

500

600

700

Figure 1. ω(RD(C)) vs. iterations for Dolphins

800

CUTTING PLANE FOR MODULARITY MAXIMIZATION

11

0.7

0.6

0.5

obj.val

0.4

0.3

0.2

0.1

0

SCP
MCP
!(P)

-0.1
0

1000

2000

3000
iteration

4000

5000

6000

Figure 2. ω(RD(C)) vs. iterations for Football
0.5

0.4

obj.val

0.3

0.2

0.1

0
SCP
MCP
!(P)

-0.1
0

5000

10000

15000
iteration

20000

25000

30000

Figure 3. ω(RD(C)) vs. iterations for Jazz
7. Conclusion
One of the key points in developing a good algorithm for modularity maximization problems
would be generating deep cutting planes. We proposed to solve the auxiliary problem AP (C)by a heuristic algorithm based on the noising method in this paper, however here still remains
room for further research. The method of multiple cutting planes that we proposed in this paper
performed fairly well, however it should need further investigation from both theoretical andcomputational view points.

12

Y. IZUNAGA AND Y. YAMAMOTO

References[1] G. Agarwal and D. Kempe, “Modularity-maximizing graph communities via mathematical programming,”
The European Physical Journal, B.66, pp.409-418, 2008.
[2] D. Aloise, S. Cafieri, G. Caporossi, P. Hansen, L. Liberti, and S. Pellon, “Column generation algorithms forexact modularity maximization in networks,” Physical Review, E.82, 2012.
[3] 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, pp.172-188, 2008.
[4] A. Caprara, M. Fischetti, and P. Toth, “Heuristic method for the set covering problem,” Operations Research,47, pp.730-743, 1999.
[5] A. Caprara, P. Toth, and M. Fischetti, “Algorithms for the set covering problem,” Annals of Operations
Research, 98, pp.353-371, 2000.
[6] I. Charon and O. Hudry, “Application of the noising method to the travelling salesman problem,” European
Journal of Operations Research, 125, pp.266-277, 2000.
[7] I. Charon and O. Hudry, “The noising methods: A generalization of some metaheuristics,” European Journalof Operations Research, 135, pp.86-101, 2001.
[8] O. du Merle, D. Villeneuve, J. Desrosiers, and P. Hansen, “Stabilized column generation,” Discrete Mathematics, 194, pp.229-237, 1999.
[9] M. Grötschel and Y. Wakabayashi, “A cutting plane algorithm for a clustering problem,” Mathematical
Programming, 45, pp.59-96, 1989.
[10] P. Hansen and N. Mladenovic, “Variable neighborhood search: Principles and applications,” European Journlof Operational Research, 130, pp.449-467, 2001.
[11] M. E. J. Newman and M. Girvan, “Finding and evaluating community structure in networks,” PhysicalReview, E.69, 2004.
[12] C. R. Reeves (ed.): Modern heuristic techniques for combinatorial problem. Blackwell Scientific Publications,1993.
[13] S. Umetani and M. Yagiura, “Relaxation heuristics for the set covering problem,” Journal of the Operations
Research Society of Japan, 50, pp.350-375, 2007.
(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(Y. Yamamoto) Faculty of Engineering, Information and Systems, University of Tsukuba, Tsukuba,
Ibaraki 305-8573, Japan
E-mail address: yamamoto@sk.tsukuba.ac.jp