

Department of Social Systems and Management
Discussion Paper Series

No.1291

Lagrangian Relaxation and Pegging Testfor the Clique Partitioning Problem
byNoriyoshi SUKEGAWA, Yoshitsugu YAMAMOTO, and Liyuan ZHANG

May 2012

UNIVERSITY OF TSUKUBA
Tsukuba, Ibaraki 305-8573
JAPANLAGRANGIAN RELAXATION AND PEGGING TESTFOR THE CLIQUE PARTITIONING PROBLEMNORIYOSHI SUKEGAWA, YOSHITSUGU YAMAMOTO, AND LIYUAN ZHANG
Abstract. The clique partitioning problem is an NP -hard combinatorial optimization
problem with applications to data analysis such as clustering. Though a binary integer
linear programming formulation has been known for years, one needs to deal with ahuge number of variables and constraints when solving a large instance. In this paper,
we propose a size reduction algorithm which is based on the Lagrangian relaxation andthe pegging test, and verify its validity through numerical experiments. We modify
the conventional subgradient method in order to manage the high dimensionality of the
Lagrangian multipliers, and also make an improvement on the ordinary pegging test by
taking advantage of the structural property of the clique partitioning problem.

1. Introduction
Given a set of objects with a number of attributes, cluster analysis aims to groupthem into clusters according to similarity or dissimilarity of the attributes. When theattributes are quantitative, plenty of choices of dissimilarity measurement are available
such as Euclidean, Minkowsky and Mahalanobis distances. On the other hand, we havefew choices for the qualitative attributes as mentioned by Brusco and Köhn [3]. An
approach would be the qualitative-quantitative data conversion, i.e., assigning numericalnumbers to the qualitative attributes. This is simple hence popular among practitioners,
but the clustering can be greatly affected by the conversion scheme.
Suppose that a qualitative attribute consists of three distinct symbols. Then it defines
an equivalence relation on the set of objects and forms three equivalence classes, i.e., twoobjects having the same symbol are equivalent. Thus if there are q attributes, we have q
different equivalence relations. Grötschel and Wakabayashi [8] proposed to aggregate thoseequivalence relations for a clustering of the objects, and formulated it to a combinatorial
optimization problem called Clique Partitioning Problem, CPP for short. They proposeda cutting plane algorithm based on the simplex method and solved real-world instances
up to 158 objects to optimality. For subsequent researches including heuristic methods,
see [1, 3, 11, 13]. In recent years, CPP has broadened its application range to otherreal-world problems such as the flight-to-gate assignment problem [5] and microarray data
analysis [9].
The CPP for n objects, when formulated as a binary integer linear programming, hasO(n2 ) binary variables, O(n2 ) equality constraints and O(n3 ) inequality constraints, allof which grow very rapidly as n grows and will soon outstrip the computing capability ofcurrently available optimizers. Therefore we need to devise means of reducing the problem
size. If we have the information that a binary variable is zero at any optimal solution,
we can fix it to zero, delete it, and reduce the number of variables. Such information
is obtained by comparing the lower bounds of the optimal value of the problem and theDate: May 6, 2012.
Key words and phrases. Clique partitioning, Clustering, NP -hard, Binary integer programming, Lagrangian dual problem, Subgradient method.
This work is partly supported by Grant-in-Aid for Scientific Research (C) No.22510136.
1

2

SUKEGAWA, YAMAMOTO, AND ZHANG

upper bound of the optimal value of the problem with a variable temporarily fixed to one.
If the latter is less than the former, we can conclude that the variable should be zeroat any optimal solution. This idea, called the pegging test, was proposed years ago, andfavorable results have been reported on the knapsack problem [12, 14] and the set covering
problem [2]. One of the aims of this paper is to propose the pegging test which is tailored
to exploit the problem structure of the CPP. To obtain the upper bound of the problem,
we use the Lagrangian relaxation problem instead of the conventional linear programming
relaxation since the problem size remains intractable even when the binary constraints arerelaxed. We combine the subgradient method for solving the Lagrangian dual problem
with some heuristics and local search.
This paper is organized as follows. In Section 2, we give the definition of CPP and
its binary integer linear programming formulation. In Section 3, we first summarize theconventional Lagrangian relaxation problem and then develop the subgradient method
for the Lagrangian dual problem so that it can handle a huge number of Lagrangian
multipliers. In Section 4, we develop a heuristic method and introduce an efficient andnumerically stable method. In Section 5, we improve the ordinary pegging test to utilize
the structural property of CPP. Describing our algorithm in Section 6, we report thenumerical experiments of our method in Section 7. Some comments are made in Section 8.
2. Clique partitioning problem
We first introduce the clique partitioning problem and its formulation as a binary integer
linear programming problem according to Grötschel and Wakabayashi [8].
2.1. Definitions and notations. We denote an undirected graph G with vertex set N
and edge set E by G = (N, E). An edge e with endnodes u and v is denoted by {u, v}.
For a set S of vertices, we denote the set of edges in G with both endnodes in S by E(S),that is, E(S) = {{u, v} ∈ E | u, v ∈ S}. For a set of subsets S1 , S2 , . . . , Sk of N , let
E(S1 , S2 , . . . , Sk ) :=k!

E(Si ).

i=1

A graph is called complete if every pair of its nodes is linked by an edge. A clique is acomplete subgraph.
We say that Γ = {W1 , W2 , . . . , Wk } is a partition of N if Wi ∩ Wj = ∅ for 1 ≤ i < j ≤ k,
V = W1 ∪W2 ∪· · ·∪Wk , and Wi &= ∅ for 1 ≤ i ≤ k. A set A of edges in a graph G = (N, E)
is called a clique partitioning of G if there is a partition Γ = {W1 , W2 , . . . , Wk } of N suchthat A = E(W1 , W2 , . . . , Wk ) and the subgraph (Wi , E(Wi ) ∩ A) induced by Wi is a cliquefor 1 ≤ i ≤ k. In case G is complete, every partition of N induces a clique partitioning.
The clique partitioning problem (CPP) is formally defined as follows [8]. Given acomplete graph K = (N, E) with weights ce ∈ R for all e ∈ E, find a clique partitioningA ⊆ E such that"c(A) :=
ceis as large as possible.

e∈A

2.2. Binary integer linear programming formulation. In this section, we show astandard formulation of CPP as a binary integer linear programming problem.
For the sake of simplicity we assume that N = {1, 2, . . . , n} in what follows. For aclique partitioning A let binary variables xij for {i, j} ∈ E be defined as#
1 {i, j} ∈ A,xij =0 otherwise,CLIQUE PARTITIONING PROBLEMS

3

then CPP is formulated as"$$ maximizec̄ij xij$$
2(i,j)∈N$!=$ subject to x ∈ {0, 1}
for all {i, j} ∈ E(binary)$ij$
xij − xji = 0
for all {i, j} ∈ E(symmetry)$$xij + xjk − xik ≤ 1 for all i, j, k ∈ N, i &= j, j &= k, k &= i, (transitivity)wherec̄ij = c̄ji :=c{i,j}
2 .

The transitivity constraints guarantee that the relation defined by those variables xij ’s
is transitive (i.e., if i and j are in the same clique and j and k are in the same clique, theni and k must be in the same clique), and this assures that the feasible solutions correspondto the clique partitioning. The point is that this formulation has n(n − 1) binary variables,n(n − 1)/2 equality constraints and n(n − 1)(n − 2)/2 inequality constraints, all of which
grow very rapidly as n grows. Substituting xij for xji for all i, j ∈ N with i < j halves
the decision variables and yields the following equivalent problem(P ) :

where$$ maximize$$$$ subject to
$
$
$
$
$
$

"

2
(i,j)∈N<cij xijxij ∈ {0, 1}
xij + xjk − xik ≤ 1
xij − xjk + xik ≤ 1−xij + xjk + xik ≤ 1

for all (i, j) ∈ N<2 ,for all (i, j, k) ∈ N<3 (type U)
for all (i, j, k) ∈ N<3 (type V)
for all (i, j, k) ∈ N<3 , (type W)

cij := c{i,j} ,N<2 := { (i, j) | i, j ∈ N, i < j },N<3 := { (i, j, k) | i, j, k ∈ N, i < j < k }.
The inequality constraints separated into three types, which we call type U, V and W,
respectively. We will denote the optimal objective function value of the problem (P ) byω(P ).
3. Lagrangian Relaxation and Subgradient Method
We propose to solve the Lagrangian relaxation problem of (P ) instead of the linear
programming relaxation in order to obtain an upper bound of ω(P ). The conventional
Lagrangian relaxation problem needs to handle Lagrangian multipliers as many as theconstraints of (P ), which can easily exceed the manipulable number. Then we propose amodified Lagrangian relaxation problem and modified subgradient method.
3.1. Lagrangian relaxation problem. One of the common tricks to deal with the problem (P ) would be the Lagrangian relaxation problem. Namely, introducing a nonnegative
multiplier uijk for each constraint of type U, vijk for type V and wijk for type W, we
consider the following integer linear programming problem

4

SUKEGAWA, YAMAMOTO, AND ZHANG(LR(u, v, w)) :"$$ maximize
cij xij$$
2(i,j)∈N
<
$
$
$
$
$
$
$
$
$
$
$ subject to x ∈ {0, 1}
ij

+

"

3
(i,j,k)∈N<+"

3(i,j,k)∈N<+"

3(i,j,k)∈N<

uijk (1 − xij − xjk + xik )vijk (1 − xij + xjk − xik )wijk (1 + xij − xjk − xik )for all (i, j) ∈ N<2

where u, v and w denote multiplier vectors (uijk )(i,j,k)∈N<3 , (vijk )(i,j,k)∈N<3 and (wijk )(i,j,k)∈N<3 ,respectively. The inequality constraints are relaxed and plugged into the objective function as a penalty. Denoting the coefficient of variable xij in the objective function byr(u, v, w)ij , the objective function of (LR(u, v, w)) is rewritten as""(uijk + vijk + wijk ) +r(u, v, w)ij xij ,3(i,j,k)∈N<

2(i,j)∈N<and the coefficient r(u, v, w)ij is given as"r(u, v, w)ij = cij −uijk −"

3
k:(i,j,k)∈N<

(3.1)

−
+

"

3
k:(i,j,k)∈N<"

3
k:(i,j,k)∈N<

3
k:(k,i,j)∈N<"

vijk +

3
k:(k,i,j)∈N<

wijk −""

ukij +uikj

3
k:(i,k,j)∈N<

vkij −

3
k:(k,i,j)∈N<

"

vikj

3
k:(i,k,j)∈N<

wkij −"

wikj .

3
k:(i,k,j)∈N<

Due to the simple constraint of (LR(u, v, w)), its optimal solution can be obtained simplyby#
1 if r(u, v, w)ij > 0,
(3.2)
x(u, v, w)ij =0 if r(u, v, w)ij ≤ 0.

The optimal objective function value ω(LR(u, v, w)) of the problem (LR(u, v, w)) provides an upper bound of ω(P ), and it clearly depends on the choice of the Lagrangian
multiplier vector (u, v, w). We propose in Section 3.2 to apply the subgradient method
to find a multiplier vector that provides the smallest value of ω(LR(u, v, w)).
Even if the optimal solution x(u, v, w) of the problem (LR(u, v, w)) is feasible to the
problem (P ), it is not necessarily optimal for the problem (P ). The following theorem iswell known. See e.g., Geoffrion [7].
Theorem 3.1. Let (ū, v̄, w̄) := ((ūijk )(i,j,k)∈N<3 , (v̄ijk )(i,j,k)∈N<3 , (w̄ijk )(i,j,k)∈N<3 ) be a Lagrangian multiplier vector, and let x be an optimal solution of the Lagrangian relaxation
problem of the problem (P ) with (ū, v̄, w̄). If x is feasible to the problem (P ) and satisfies
the complementarity condition :ūijk (1 − xij − xjk + xik ) = 0v̄ijk (1 − xij + xjk − xik ) = 0w̄ijk (1 + xij − xjk − xik ) = 0

then it is an optimal solution of the problem (P ).

for all (i, j, k) ∈ N<3 ,for all (i, j, k) ∈ N<3 ,for all (i, j, k) ∈ N<3 ,CLIQUE PARTITIONING PROBLEMS

5

3.2. Subgradient method. For the sake of notational simplicity we abbreviate ω(LR(u, v, w))to ω(u, v, w) within this section. The Lagrangian dual problem, denoted by (LD), is aproblem for finding the smallest value of ω(u, v, w), i.e., the best upper bound of ω(P ).
Namely,(LD) :$$ minimize ω(u, v, w)$$ subject to u, v, w ≥ 0.

The objective function ω(u, v, w) is convex, piecewise linear and not differentiable on theintersection of pieces, and one of the most widely used methods for this problem is the
subgradient method. See for example Fisher [6].
Definition 3.2. (g u , g v , g w ) is said to be a subgradient of ω at (ū, v̄, w̄) ≥ 0 whenω(ū, v̄, w̄) + *g u , u − ū+ + *g v , v − v̄+ + *g w , w − w̄+ ≤ ω(u, v, w)

holds for any (u, v, w) ≥ 0, where *·, ·+ means the inner product.
The following lemma is well known.

Lemma 3.3. Let x(u, v, w) denote an optimal solution of the Lagrangian relaxation problem (LR(u, v, w)). Then (g u , g v , g w ) defined asu
gijk:= 1 − x(u, v, w)ij − x(u, v, w)jk + x(u, v, w)ik

for (i, j, k) ∈ N<3 ,w
gijk:= 1 + x(u, v, w)ij − x(u, v, w)jk − x(u, v, w)ik

for (i, j, k) ∈ N<3v
gijk:= 1 − x(u, v, w)ij + x(u, v, w)jk − x(u, v, w)ik

is a subgradient of ω at (u, v, w).

for (i, j, k) ∈ N<3 ,The subgradient method uses the following rule to update the iterate of multiplier vector(u, v, w) to the next iterate (u+ , v + , w+ ).
%&ω(u, v, w) − ωlow u+(3.3)
uijk := max 0, uijk − µgfor (i, j, k) ∈ N<3 ,,(g u , g v , g w ),2 ijk%&
ω(u, v, w) − ωlow v+(3.4)
vijk := max 0, vijk − µgfor (i, j, k) ∈ N<3 ,,(g u , g v , g w ),2 ijk%&
ω(u, v, w) − ωlow w+(3.5)
wijk := max 0, wijk − µgfor (i, j, k) ∈ N<3 ,,(g u , g v , g w ),2 ijk
where µ is a step size control parameter and ωlow is a lower bound of ω(P ). It is known
that if ωlow is replaced by the optimal value ω(P ), the sequence generated will converge
to an optimal solution of the Lagrangian dual problem (LD), see e.g., Geoffrion [7] andLarsson et al. [10]. Therefore the quality of lower bound ωlow has a decisive influenceon the convergence of the method. The value ω(u, v, w) does not necessarily decrease
when the multiplier vector is updated. Often the parameter µ is initially set to 2.0 and ishalved whenever ω(u, v, w) fails to decrease in some fixed number of iterations. This rule
performs well empirically. The method we used is as follows.

Basic Subgradient Method (BSM)
Step 1 (Initialization)(a) l ← 0, µ ← 2.0, (u, v, w) ← (0, 0, 0).
(b) ωlow ← the objective value of a feasible solution obtained by some heuristic
algorithm (see Section 4). ωup ← ω(u, v, w).

6

SUKEGAWA, YAMAMOTO, AND ZHANG

Step 2 (Solving (LR(u, v, w))(a) Compute r(u, v, w)ij by (3.1) and set x(u, v, w)ij according to (3.2).
(b) ωup ← min{ωup , ω(u, v, w)}. If ωup is not improved, l ← l + 1. Otherwise,l ← 0.
Step 3 (Optimality check)
If x(u, v, w) satisfies the optimality condition in Corollary 3.4 with (u, v, w)(in this case, x(u, v, w) is optimal for the problem (P )), then terminate. If
ωup − ωlow < ε, then terminate where ε is a predetermined tolerance for the
duality gap.
Step 4 (Update of µ)
If µ ≤ 0.005, then terminate. If l reaches 30, µ ← µ/2.
Step 5 (Update of (u, v, w))
Update (u, v, w) according to (3.3), (3.4) and (3.5) and go to Step 2.
3.3. Modified Lagrangian relaxation. Since the Lagrangian dual problem (LD) has
to handle the multipliers as many as the constraints of the problem (P ), which amountto n(n − 1)(n − 2)/2, we need to reserve a huge capacity of memory to solve (LD) fora large CPP instance. To overcome this, we present a modified Lagrangian relaxation
problem, which temporarily discards some of the transitivity constraints and retrieves
when needed. Let U, V and W be subsets of N<3 , and (u, v, w) be a nonnegative multipliervector respectively. Then the problem we consider is(LR(U, V, W, u, v, w)) :$"$ maximize
cij xij$$
2(i,j)∈N<
$
$
$
$
$
$
$
$
$
$
$ subject to xij ∈ {0, 1}+"(i,j,k)∈U+"(i,j,k)∈V+"(i,j,k)∈W

uijk (1 − xij − xjk + xik )vijk (1 − xij + xjk − xik )wijk (1 + xij − xjk − xik )for all (i, j) ∈ N<2 .

This modified Lagrangian relaxation problem differs from the ordinary one in that themultipliers with indices not in (U, V, W ) are fixed to zero. This reduces the number ofmultipliers that we have to handle to |U | + |V | + |W |. The modified Lagrangian relaxation
problem could be viewed as an ordinary Lagrangian relaxation problem of the problem(P ) with some transitivity constraints being relaxed:(P (U, V, W )) :$$ maximize$$$$ subject to
$
$
$
$
$
$

"

2
(i,j)∈N<cij xijxij ∈ {0, 1}
xij + xjk − xik ≤ 1
xij − xjk + xik ≤ 1−xij + xjk + xik ≤ 1

for all (i, j) ∈ N<2 ,for all (i, j, k) ∈ U ,for all (i, j, k) ∈ V ,for all (i, j, k) ∈ W .

It is likely that some of the discarded inequality constraints would be satisfied. Thus, if
we can find a “good and small” triplet of (U, V, W ), the modified Lagrangian relaxation

CLIQUE PARTITIONING PROBLEMS

7

problem would be in performance comparable to and in computation better than theordinary one.
For the problem (LR(U, V, W, u, v, w)) the coefficient r(u, v, w)ij of variable xij in theobjective function is written as"""r(u, v, w)ij = cij −uijk −
ukij +
uikj
k:(i,j,k)∈U

(3.6)

−

"

k:(i,j,k)∈V+"

k:(i,j,k)∈W

k:(k,i,j)∈U"

vijk +

k:(k,i,j)∈V

wijk −"

k:(i,k,j)∈U

vkij −

k:(k,i,j)∈W

"

vikj

k:(i,k,j)∈V

wkij −"

wikj .

k:(i,k,j)∈W

and an optimal solution x(u, v, w) = (x(u, v, w)ij )(i,j)∈N<2 can be obtained by#
1 if r(u, v, w)ij > 0,
(3.7)
x(u, v, w)ij =0 if r(u, v, w)ij ≤ 0.
Similar to Theorem 3.1, the optimality condition is given as follows.

Corollary 3.4. If an optimal solution x(u, v, w) of the problem (LR(U, V, W, u, v, w))
is feasible to the problem (P ) and satisfies the restricted complementarity condition with(u, v, w):
uijk (1 − xij − xjk + xik ) = 0
vijk (1 − xij + xjk − xik ) = 0

wijk (1 + xij − xjk − xik ) = 0

then it is an optimal solution of the problem (P ).

for all (i, j, k) ∈ U ,for all (i, j, k) ∈ V ,for all (i, j, k) ∈ W ,Proof. We readily see that the Lagrangian relaxation problem (LR(U, V, W, u, v, w)) is an
ordinary Lagrangian relaxation problem for (ū, v̄, w̄) such that#uijk for (i, j, k) ∈ U ,ūijk =0for (i, j, k) ∈ N<3 \ U ,#vijk for (i, j, k) ∈ V ,v̄ijk =0for (i, j, k) ∈ N<3 \ V ,#wijk for (i, j, k) ∈ W ,w̄ijk =0for (i, j, k) ∈ N<3 \ W .

When x(u, v, w) meets the restricted complementarity condition with (u, v, w), it alsosatisfies the complementarity condition for all constraints with (ū, v̄, w̄). This together
with the feasibility of x(u, v, w) and Theorem 3.1 yields the desired result.
!
Let us abbreviate ω(LR(U, V, W, u, v, w)) to ω(U, V, W, u, v, w). Accordingly, we modify the Lagrangian dual problem as follows.
(LD(U, V, W )) :$$ minimize$$ subject to
$
$
$
$
$
$

ω(U, V, W, u, v, w)
u, v, w ≥ 0
uijk = 0 for all (i, j, k) ∈ N<3 \Uvijk = 0 for all (i, j, k) ∈ N<3 \Vwijk = 0 for all (i, j, k) ∈ N<3 \W .

8

SUKEGAWA, YAMAMOTO, AND ZHANG

In spite of the reduction in the number of multipliers, the problem of finding a “goodand small” triplet of (U, V, W ) remains unsolved. To solve the modified Lagrangian dual
problem, we propose the following modified subgradient method (MSM for short) including
the subgradient method and the updating scheme of the (U, V, W ). The MSM is composedof the inner and outer cycles. The inner cycle consisting of Step 2 to 6 generates a sequenceof Lagrangian multiplier vectors, and the outer cycle expands the constraint index set(U, V, W ).
Modified Subgradient Method (MSM)
Step 1 (Initialization)(a) k, l, q ← 0, µ ← 2.0, (u, v, w) ← (0, 0, 0).
(b) ωlow ← the objective value of a feasible solution obtained by some heuristic algorithm (see Section 4). ωup , ω̄up ← ω((LR(∅, ∅, ∅, u, v, w)) and (U, V, W ) ←the transitivity constraints violated by x(u, v, w).
Step 2 (Solving (LR(U, V, W, u, v, w))(a) Compute r(u, v, w)ij by (3.6) and set x(u, v, w)ij according to (3.7).
(b) ωup ← min{ωup , ω(LR(U, V, W, u, v, w))}. If ωup is not improved, l ← l + 1.
Otherwise, l ← 0.
Step 3 (Optimality check)
If x(u, v, w) satisfies the optimality condition in Theorem 3.4 with (u, v, w) (inthis case, x(u, v, w) is optimal for the problem (P )), then terminate.
Step 4 (Update of µ)(a) If µ ≤ 0.005, then µ ← 2.0 and go to Step 6 (we decide that there is nochance of improving the upper bound unless we expand (U, V, W )).
(b) If l reaches 5, µ ← µ/2.
Step 5 (Update of (u, v, w))
Update (u, v, w) according to (3.3), (3.4) and (3.5) and go to Step 2.
Step 6 (Lagrangian heuristics and local search)
Apply the Lagrangian heuristics and local search (see Section 4) to x(u, v, w) fora better solution x̃. ωlow ← max{ωlow , objective function value of x̃}.
Step 7 (Termination)(a) If the improvement of ωup compared to ω̄up is less than 0.1% or 1 in value
then q ← q + 1. ω̄up ← ωup(b) If q > 10 then terminate.
Step 8 (Update of (U, V, W ))(a) Find the transitivity constraints violated by x(u, v, w) and add them to U ,V and W .
(b) uijk , vijk , wijk ← 0 for newly added indices (i, j, k) and go to Step 2.
In Step 6, we apply the Lagrangian heuristics and the local search to make a good lower
bound from the optimal solution of the problem (LR(U, V, W, u, v, w)). The details will
be given in Section 4.

3.4. Theoretical aspect. We will show that the problem (LD) always has an optimal
solution having at most n(n − 1)/2 positive elements. This supports the adequacy of theMSM. We first focus on the optimal solution of the linear programming relaxation problemof (P ) given as follows.

CLIQUE PARTITIONING PROBLEMS(P̄ ) :$$ maximize$$$$ subject to
$
$
$
$
$
$

"

9

cij xij

2(i,j)∈N<0 ≤ xij ≤ 1
xij + xjk − xik ≤ 1
xij − xjk + xik ≤ 1−xij + xjk + xik ≤ 1

for all (i, j) ∈ N<2 ,for all (i, j, k) ∈ N<3 ,for all (i, j, k) ∈ N<3 ,for all (i, j, k) ∈ N<3 .

Its dual problem is(D̄) :""$$ maximize(uijk + vijk + wijk ) +
zij
$
$
3
2
(i,j,k)∈N<(i,j)∈N<$$ subject to u , v , w ≥ 0
ijk
ijk
ijk
$
$
zij ≥ 0$u + Dv + Dw + z ≥ c
$
Dij
ij
ij
ij
ij

where

"

u :=
Dij

uijk

+

3
(i,j,k)∈N<

v :=
Dij

"

"

ukij

3
(k,i,j)∈N<

vijk

3(i,j,k)∈N<

w := −
Dij"−

wijk +

3(i,j,k)∈N<

"

vkij

−

3
(k,i,j)∈N<"

uikj ,

3(i,k,j)∈N<"+

3(k,i,j)∈N<"

for all (i, j, k) ∈ N<3 ,for all (i, j) ∈ N<2 ,for all (i, j) ∈ N<2 ,vikj ,

3(i,k,j)∈N<

wkij

+

"

wikj .

3(i,k,j)∈N<

Since the Lagrangian relaxation problem (LR(u, v, w)) has the integrality property (wesay that a problem has the integrality property if the optimal value is unchanged when
integrality restriction is removed.), we readily seeω(P̄ ) = ω(D̄)= ω(LD)(by the duality theorem)(by the integrality property).

Further, if a multiplier vector (u∗ , v ∗ , w∗ , z ∗ ) is optimal for (D̄), then the multiplier vector(u∗ , v ∗ , w∗ ) is optimal for (LD) (see for instance Geoffrion [7]). Since the problem (D̄) isa linear programming problem, there exists a basic optimal solution (u∗ , v ∗ , w∗ , z ∗ ), andthe number of its positive elements is at most the number of inequality constraints, that is|N<2 | = n(n − 1)/2. Therefore the Lagrangian dual problem (LD) has an optimal solution
whose positive elements are at most n(n − 1)/2.
4. Heuristics
Finding a good lower bound ωlow is essential for the convergence of the subgradientmethod, and as we will see in Section 5, it also enhances the capabilities of the pegging
test. For this purpose we use the Forgotten Vertices Method with Variant D, FVD for
short, proposed by Charon and Hudry [4]. We repeatedly solve the Lagrangian relaxation
problem for different Lagrangian multiplier vectors in the MSM. The optimal solution
thus obtained may not be feasible to the problem (P ), however it provides us with useful
information to make a good feasible solution to (P ), hence a good lower bound. We firstdescribe the Lagrangian heuristics, then introduce the local search and FVD.

10

SUKEGAWA, YAMAMOTO, AND ZHANG

4.1. Lagrangian heuristics. Let K = (N, E) denote the complete undirected graphwith n vertices. For a given S ⊆ N of k vertices and A ⊆ E letIO vA (S) + vA (S) (k ≥ 2),
O (S)vA (S) :=vA(k = 1),+∞(k = 0),where|{ {i, j} ∈ E | {i, j} ∈/ A, i, j ∈ S }|,k(k − 1)|{ {i, j} ∈ E | {i, j} ∈ A, i ∈ S, j ∈/ S }|OvA(S) :=.
k(n − k)We readily see that if A is a clique partitioning of K = (N, E) induced by a partition
Γ = {W1 , W2 , . . . , Wk } of N , then vA (Wi ) = 0 for all i = 1, 2, . . . , k. For an optimal
solution x(u, v, w) of the problem (LR(U, V, W, u, v, w)), let Ā be a set of edges such that
IvA(S) :=

Ā := { {i, j} ∈ En | x(u, v, w)ij = 1 }.

It would be reasonable to think that a partition Γ = {W1 , W2 , . . . , Wk } of N minimizing
k
"

(4.1)

vĀ (Wi )i=1

is a feasible solution “nearest” to x(u, v, w). Here we use the following simple greedy
algorithm to minimize (4.1).
Basic Lagrangian Heuristics (BLH)
Step 1 (Initialization)(a) k := 1, Wk := {}, Γ := {Wk }(b) A := Ā, V := N , E := En and F :=N .
Step 2 (Local Search)(a) Compute vA (Wk ∪ {i}) for all i ∈ F .
(b) i∗ := argmin{ vA (Wk ∪ {i}) | i ∈ F }.
Step 3 (Update of Γ )(a) If vA (Wk ∪ {i∗ }) ≤ vA (Wk ) then Wk := Wk ∪ {i∗ }, F := F \ {i∗ } and go to
Step 4 if F = ∅ and go to Step 2 otherwise.
(b) Otherwise, V := V \ {Wk }, E := E(V ), A := A ∩ E, k := k + 1, Wk := {},Γ := Γ ∪ {Wk } and go to Step 4.
Step 4 (Termination) Output Γ and terminates.
The feasible solution obtained by BLH will be used as an initial feasible solution for the
local search to start with.
4.2. Local search. Given a partition Γ = {W1 , W2 , . . . , Wk }, we define the neighborhood
N (Γ ) of Γ as followsn!
N (Γ ) :=Ni (Γ )
where%
Ni (Γ ) := Γ #i=1

Γ # is obtained from Γ by moving i from its current set Ws
to another one Wt ∈ (Γ ∪ {∅}) \ {Ws }&.

This neighborhood is equivalent to that suggested by Charon and Hudry [4]. Since theincidence vector of a partition Γ is feasible for the problem (P ), in what follows, let f (Γ )

CLIQUE PARTITIONING PROBLEMS

11

denote the objective function value of Γ .
Definition 4.1. Given a partition Γ = {W1 , W2 , . . . , Wk } and a vertex i ∈ Vn , let Γi∗ ∈Ni (Γ ) denote the best partition obtained by moving i, i.e.,+,Γi∗ = argmax f (Γ # ) | Γ # ∈ Ni (Γ ) .
We say that i is in the best position if f (Γi∗ ) ≤ f (Γ ) holds.

To find a better feasible solution, we employ the following local search method basedon the neighborhood defined above.
Basic Local Search (BLS)
Step 1 (Initialization) Set Γ ← an initial feasible solution for the problem (P ).
Step 2 (Decide exploration order) Rank the vertices in a random order.
Step 3 (Exploration)(a) opt:= true.
(b) According to the order defined in Step 2, for all vertex i, if i is not in thebest position, opt:= false and Γ ← Γi∗ .
Step 4 (Termination)(a) If opt is false, then go to Step 2.
(b) Otherwise output Γ (in this case, Γ is a local optimal solution) and terminate.
4.3. Forgotten vertices method with variant D. The Lagrangian heuristics describedin the previous section seems to be successful when the corresponding multiplier vector
is near optimal. To obtain a good lower bound in an early stage of the computation, we
introduce a heuristic method called Forgotten vertices method with variant D (FVD for
short) by Charon and Hudry [4]. This method is mainly based on the above local searchbut sometimes forgets some vertices when calculating the objective function value. Morespecifically, when some vertices are forgotten, the perturbed objective function value fp ofa partition Γ = {W1 , W2 , . . . , Wk } is given by"
fp (Γ ) =ce
e={u,v}∈Au,v not forgotten

where A is a clique partitioning of Vn induced by Γ .
Definition 4.2. Given a partition Γ = {W1 , W2 , . . . , Wk } and a vertex i ∈ Vn , let Γ̃i∗ ∈Ni (Γ ) denote the perturbed best partition obtained by moving i, i.e.,+,Γ̃i∗ = argmax fp (Γ # ) | Γ # ∈ Ni (Γ ) .
We say that i is in the perturbed best position if fp (Γi∗ ) ≤ fp (Γ ) holds.

Forgetting some vertices can be considered as a “noise.” Charon and Hudry [4] proposed
18 variants of noise and reported that the following variant D was always the best andnumerically stable.
Forgotten vertices method with variant D (FVD)
Step 1 (Initialization)
Set Γ, Γ ∗ ← some partitioning which is randomly computed and r ← rmax .
Step 2 (Initialization for perturbed exploration) l ← 1.
Step 3 (Decide forgotten vertices and exploration order for perturbed exploration)

12

SUKEGAWA, YAMAMOTO, AND ZHANG(a) k ← the product of r by the number of vertices n. Select k vertices which
will be forgotten in Step 3.
(b) Rank the vertices in a random order.
Step 4 (Perturbed exploration)(a) According to the above order, for all vertex i, if i is not in the perturbed best
position, Γ ← Γ̃i∗ .
(b) If l ≥ 4 then go to Step 5. Otherwise l ← l + 1 and go to Step 3.
Step 5 (Decide exploration order for unperturbed exploration)
Rank the vertices in a random order.
Step 6 (Unperturbed exploration)(a) opt:= true.
(b) According to the above order, for all vertex i, if i is not in the best position,
opt:= false and Γ ← Γi∗ .
Step 7 (Termination for unperturbed exploration)(a) If opt is false, then go to Step 6.
(b) Otherwise, if the objective function value of Γ is larger than that of Γ ∗ , thenΓ ∗ ← Γ (Update the best solution).
Step 8 (Termination and update of r)(a) r ← r − d.
(b) If r > 0, then go to Step 2. Otherwise, output Γ ∗ and terminates.
In our algorithm, we set rmax ← 0.5 and d ← 0.05. See Charon and Hudry [4] for thedetails.
5. Pegging test
We start this section with the explanation of an ordinary pegging test, and then present
an idea for an improvement. For the sake of simplicity we abbreviate (LR(U, V, W, u, v, w))
to (LR(u, v, w)) within this section.
5.1. Ordinary pegging test. By the information obtained from the optimal solution of
the Lagrangian relaxation problem (LR(u, v, w)) we can see which variable takes one andwhich takes zero at the optimal solution of the problem (P ).
Let us choose (s, t) ∈ N<2 and suppose that the problem (P ) has an optimal solution
with xst = ξ, where ξ is either zero or one. Then the problem (P ) with an additional
constraint xst = ξ is equivalent to the problem (P ) in the sense that optimal values of thetwo problems coincide. Suppose further we have an incumbent value ωlow . Then clearlyω(P |xst = ξ) = ω(P ) ≥ ωlow .

Since the problem (LR(u, v, w)|xst = ξ) is a relaxation problem of the problem (P |xst = ξ)
we obtain
ω(LR(u, v, w)|xst = ξ) ≥ ω(P |xst = ξ),hence
ω(LR(u, v, w)|xst = ξ) ≥ ωlow .

Lemma 5.1. Let ξ be either zero or one. If ω(LR(u, v, w)|xst = ξ) < ωlow , then xst =1 − ξ for any optimal solution of the problem (P ).
Proof. The proof is clear from the above discussion.

!

Suppose that we have an optimal solution x(u, v, w) of the problem (LR(u, v, w)) andthat x(u, v, w)st = 0. By a simple calculation we see that
(5.1)

ω(LR(u, v, w)|xst = 1) = ω(LR(u, v, w)) + r(u, v, w)st .

CLIQUE PARTITIONING PROBLEMS

13

Note that x(u, v, w)st = 0 implies r(u, v, w)st ≤ 0. In the same way we see that
(5.2)

ω(LR(u, v, w)|xst = 0) = ω(LR(u, v, w)) − r(u, v, w)stwhen x(u, v, w)st = 1. Note also that r(u, v, w)st > 0 in this case.

Theorem 5.2 (ordinary pegging test). Let x(u, v, w) be an optimal solution of the Lagrangian relaxation problem LR(u, v, w). If
ω(LR(u, v, w)) − ωlow < |r(u, v, w)st |∗holds, then xst = x(u, v, w)st for any optimal solution x∗ of the problem (P ).
Proof. Substituting equation (5.1) or (5.2) for the condition in Lemma 5.1 will yield theassertion.
!
We say that the variable xst is pegged at x(u, v, w)st when the case holds in the theorem.
5.2. Improved pegging test. As the computation goes, we will have several variables
pegged. Let P0 and P1 denote the index sets of the variables that have been pegged atzero and one, respectively. Given Lagrangian multiplier vector (u, v, w), the problem(LR(u, v, w, P0 , P1 )) :$""$ maximizecx+uijk (1 − xij − xjk + xik )
ij
ij
$
$
2
(i,j,k)∈U(i,j)∈N<$"$$+
vijk (1 − xij + xjk − xik )$$(i,j,k)∈V"$$+
wijk (1 + xij − xjk − xik )$$(i,j,k)∈W$$ subject to xij ∈ {0, 1}
for all (i, j) ∈ N<2 ,$%$0
for all (i, j) ∈ P0 ,
$
xij =
$
1
for all (i, j) ∈ P1 .
is a relaxation problem of the problem (P ).
Let E(P1 ) be the set of edges {i, j} such that xij or xji has been pegged at one, i.e.,(5.3)

E(P1 ) := { {i, j} ∈ En | (i, j) ∈ P1 or (j, i) ∈ P1 }.

Definition 5.3. Given P1 and i, j ∈ N with i &= j, we say that i is connected to j when
there is a path from i to j on the edge set E(P1 ).
Definition 5.4. Given (s, t) ∈ N<2 \ (P0 ∪ P1 ) let

S= := {s} ∪ { i ∈ N | i is connected to s },T= := {t} ∪ { j ∈ N | j is connected to t },ST= := (S= × T= ) ∪ (T= × S= ).

Take a variable xst that has not yet been pegged, i.e., (s, t) ∈ N<2 \ (P0 ∪ P1 ), and fix
xst temporarily to one. Then every element connected to s should be connected to every
element connected to t by the transitivity. Namely, the variables must satisfy(5.4)

xij = 1for all (i, j) ∈ ST= ∩ N<2 .

When xst is fixed temporarily to zero, we have similarly(5.5)

xij = 0for all (i, j) ∈ ST= ∩ N<2 .

Now given nonnegative multiplier vectors u, v and w we define

14

SUKEGAWA, YAMAMOTO, AND ZHANG(LR(u, v, w, P0 , P1 )|xst = 1) :$"$ maximize
cij xij$$
2(i,j)∈N<
$
$
$
$
$
$
$
$
$
$
$ subject to xij ∈ {0, 1}
$
%
$
0
$
xij =
$
1

+

"

(i,j,k)∈U+"(i,j,k)∈V+"(i,j,k)∈W

uijk (1 − xij − xjk + xik )vijk (1 − xij + xjk − xik )wijk (1 + xij − xjk − xik )for all (i, j) ∈ N<2 ,for all (i, j) ∈ P0 ,for all (i, j) ∈ (ST= ∩ N<2 ) ∪ P1 .

This is a relaxation problem of the problem (P ) with a temporary constraint xst = 1
added.
Lemma 5.5. If (ST= ∩ N<2 ) ∩ P0 &= ∅, then x∗st = 0 for any optimal solution x∗ of the
problem (P ).
Proof. Suppose there is an element (i, j) in the set (ST= ∩ N<2 ) ∩ P0 . If (i, j) = (s, t),since (s, t) ∈ P0 , the assertion holds. If i &= s, by the transitivity constraint, we see thatthe variable corresponding to {i, s} must be one for any optimal solution. In a similarway, if j &= t, we see that the variable corresponding to {j, t} must be one for any optimal
solution. Thus, since xij ∈ P0 , to meet the transitivity constraint, xst must be zero forany optimal solution.
!
When xst is temporarily fixed to zero, we have the following problem and lemma.
(LR(u, v, w, P0 , P1 )|xst = 0) :$"$ maximize
cij xij$$
2(i,j)∈N<
$
$
$
$
$
$
$
$
$
$
$ subject to xij ∈ {0, 1}
$
%
$
0
$
x
=
ij
$
1

+

"

(i,j,k)∈U+"(i,j,k)∈V+"(i,j,k)∈W

uijk (1 − xij − xjk + xik )vijk (1 − xij + xjk − xik )wijk (1 + xij − xjk − xik )for all (i, j) ∈ N<2 ,for all (i, j) ∈ (ST= ∩ N<2 ) ∪ P0 ,for all (i, j) ∈ P1 .

Lemma 5.6. If (ST= ∩ N<2 ) ∩ P1 &= ∅, then x∗st = 1 for any optimal solution x∗ of the
problem (P ).
Since the problem (LR(u, v, wP0 , P1 )|xst = ξ) is a relaxation problem of the problem(P ) with a constraint xst = ξ added, we readily see the following lemma.
Lemma 5.7. Let ξ be either zero or one, and let xst be a variable that has not been pegged,i.e., (s, t) ∈ N<2 \ (P0 ∪ P1 ). If ω(LR(u, v, w, P0 , P1 )|xst = ξ) < ωlow , then xst = 1 − ξ forany optimal solution of the problem (P ).

CLIQUE PARTITIONING PROBLEMS

15

We have seen in (5.1) and (5.2) that
ω(LR(u, v, w)) − ω(LR(u, v, w)|xst = 1 − x(u, v, w)st ) = |r(u, v, w)st |holds. Namely, the objective function value deteriorates by |r(u, v, w)st | when the constraint xst = 1−x(u, v, w)st is added to the problem (LR(u, v, w)). In the similar manner
we see that"
ω(LR(u, v, w)) − ω(LR(u, v, w, P0 , P1 )|xst = 1 − x(u, v, w)st ) =|r(u, v, w)ij |,(i,j)∈R

where

R = ((ST= ∩ N<2 ) ∪ P1 ) ∩ {(i, j) | x(u, v, w)ij = 0}
(5.6)

when x(u, v, w)st = 0, andR = ((ST= ∩ N<2 ) ∪ P0 ) ∩ {(i, j) | x(u, v, w)ij = 1}
when x(u, v, w)st = 1.

The first subset of R of (5.6) corresponds to the variables that should be one but takes
zero at x(u, v, w), and the second subset to those that should be zero but takes one atx(u, v, w).
Theorem 5.8 (improved pegging test). Let ξ be either zero or one, and let xst be a variable
that has not been pegged, i.e., (s, t) ∈ N<2 \ (P0 ∪ P1 ). If"
ω(LR(u, v, w)) − ωlow <|r(u, v, w)ij |(i,j)∈R

holds, then xst = 1 − ξ for any optimal solution of the problem (P ).
Proof. The proof is clear from the above discussion.

!

Compared to the ordinary pegging test in Theorem 5.2, since the set R always contains(s, t), we see that"|r(u, v, w)st | ≤|r(u, v, w)ij |.
(i,j)∈R

Thus, given a multiplier vector (u, v, w), if a variable can be pegged by the ordinary pegging test then the variable can also be pegged- by the improved pegging test. In addition,
the difference between |r(u, v, w)st | and (i,j)∈R |r(u, v, w)ij | would be significant when
we have many variables pegged.
5.3. Post-processing of the pegging test. When xij and xjk have been pegged at
one by the pegging test, xik should be pegged at one to meet the transitivity constraint.
We do this post-processing by finding the connected components of the undirected graph(N, E(P1 )) through the depth-first search.
Step 1 Find the connected components {C1 , C2 , . . . , Ck } of the graph (N, E(P1 )) by thedepth-first search.
Step 2 For all pairs i and j with 1 ≤ i < j ≤ n, if they are in the same component in C
then we peg xij at one.
Further, suppose xij is pegged at zero and i and j are in different components Cp and Cq
respectively. From the transitivity constraints, we can peg xi" j " at zero for all pairs (i# , j # )in either Cp × Cq or Cq × Cp .
Step 3 For all pairs p and q with 1 ≤ p < q ≤ k, if there exists a variable xij which is
pegged at zero with (i, j) ∈ Cp × Cq or (i, j) ∈ Cq × Cp then we peg xi" j " at zerofor all pairs (i# , j # ) in either Cp × Cq or Cq × Cp .

16

by

SUKEGAWA, YAMAMOTO, AND ZHANG

Once we execute the above three steps, the two sets in Definition 5.4 are readily obtained

S= := {s} ∪ { i ∈ N | {s, i} ∈ Ē },T= := {t} ∪ { j ∈ N | {t, j} ∈ Ē }.
When some variables are newly pegged, we might peg some other variables by applying
the above three steps. Thus, we apply the three steps after every pegging test in ouralgorithm. See the next section for the details.
6. Our algorithm
The algorithm we propose in this paper is a combination of the MSM (described in
Section 3) and the pegging test (described in Section 5). We abbreviate it to MSMPT.
Modified Subgradient Method with Pegging Test (MSMPT)
Step 1 (Initialization)(a) k, l, q ← 0, µ ← 2.0, (u, v, w) ← (0, 0, 0).
(b) ωlow ← the objective value of a feasible solution obtained by FVD. ωup , ω̄up ←ω((LR(∅, ∅, ∅, u, v, w)) and (U, V, W ) ← the transitivity constraints violatedby x(u, v, w).
(c) P0 , P1 ← ∅.
Step 2 (Solving (LR(U, V, W, u, v, w, P0 , P1 ))(a) Compute r(u, v, w)ij by (3.6) and set x(u, v, w)ij according to (3.7).
(b) ωup ← min{ωup , ω(LR(U, V, W, u, v, w, P0 , P1 ))}. If ωup is not improved,l ← l + 1. Otherwise, l ← 0.
Step 3 (Optimality check)
If x(u, v, w) satisfies the optimality condition in Corollary 3.4 with (u, v, w) (inthis case, x(u, v, w) is optimal for the problem (P )), then terminate.
Step 4 (Update of µ)(a) If µ ≤ 0.005, then µ ← 2.0 and go to Step 6 (we decide that there is nochance of improving the upper bound unless we expand (U, V, W )).
(b) If l reaches 5, µ ← µ/2.
Step 5 (Update of (u, v, w))
Update (u, v, w) according to (3.3), (3.4) and (3.5) and go to Step 2.
Step 6 (Lagrangian heuristics and local search)
Apply the Lagrangian heuristics and local search to x(u, v, w) for a better solutionx̃. ωlow ← max{ωlow , objective function value of x̃}.
Step 7 (Pegging test)(a) If |P0 | = |P1 | = 0, then apply the ordinary pegging test (Theorem 5.2).
Otherwise, apply the improved pegging test (Theorem 5.8).
(b) Append the indices corresponding variables are newly pegged to P0 and P1
respectively.
(c) Apply the post-processing described in Section 5.
Step 8 (Termination)(a) If the improvement of ωup compared to ω̄up is less than 0.1% or 1 in value
then q ← q + 1. ω̄up ← ωup(b) If q > 10 then go to Step 10.
Step 9 (Update of (U, V, W ))(a) Find the transitivity constraints violated by x(u, v, w) and add them to U ,V and W .

CLIQUE PARTITIONING PROBLEMS

17(b) uijk , vijk , wijk ← 0 for newly added indices (i, j, k) and go to Step 2.
Since the pegged variables contribute to downscaling the feasible region of the Lagrangian
relaxation problem, possibly the more we peg the variables, the tighter the upper boundωup would be. In the following section, we will report numerical experiments of the
MSMPT.
7. Numerical experiments
In this section, we report numerical experiments with our algorithm described in Section 6. We coded the algorithm in Java, and run it on a PC with an Intel i3, 3.33 GHz
processor and 2 GB of memory.
Since all of the instances we solved have integral objective function, we stop the MSMPT
if the duality gap fell to below 1 in value (in this case an optimal solution is found).
7.1. Real-world instances. The real-world instances we solved are problems of aggregation of equivalence relations, which are available from Grötschel and Wakabayashi [8].
Suppose we have q different equivalence relations ∼1 , ∼2 , . . . , ∼q on N , and want toaggregate them into a single equivalence relation ∼X on N . Here, we define the “distance”
d(∼k , ∼X ) of two equivalence relations ∼k and ∼X as the total number of disagreementof ∼k with ∼X i.e.,d(∼k , ∼X ) := |{{i, j} | (i ∼k j ∧ i &∼X j) ∨ (i &∼k j ∧ i ∼X j) }|.

By representing an equivalence relation ∼1 , ∼2 , . . . , ∼q and ∼X as a binary vector i.e.,##
1 i ∼k j,
1 i ∼X j,krij =for all k ∈ {1, 2, . . . , q} and xij =0 otherwise,
0 otherwise,
we see that the quantity d(∼k , ∼X ) can be given by/2" .
kd(∼k , ∼X ) =
rij
− xij
2
(i,j)∈N<="2(i,j)∈N<="2(i,j)∈N<kkk , x ∈ {0, 1})rij− 2rijxij + xij (∵ rij
ij

.
/
k
k
rij
+ 1 − 2rij
xij .

The problem of aggregation of equivalence relations is to find an equivalence relation ∼X
minimizing the total dissimilarity


q
p
"
"
"
k
k

d(∼k , ∼X ) =rij+ (1 − 2rij
)xij 
k=1

=

k=1

2
(i,j)∈N<p""2
k=1 (i,j)∈N<

= c+

"

2
(i,j)∈N<

= c+

"

2
(i,j)∈N<

k
rij
+

4 p
"
k=1

p
"

"

2
k=1 (i,j)∈N<k(1 − 2rij
)xij

5

k
(1 − 2rij) xij

cij xij .

18

SUKEGAWA, YAMAMOTO, AND ZHANG

Table 1. Result for the instances described in Grötschel and Wakabayashi [8]instancen|N<2 ||N<3 | relative gap%!
%" time(s)
A1
36
630
21,4200 1.18
0
0.01
A2
30
435
12,1800 6.08 91.49
0.04
A3
34
561
17,952
1.E-03 13.27 91.44
1.96
A4
33
528
16,368
0 14.4000.16
A6
54 1,431
74,4120 3.140
0.03
A7.1
158 12,403 1,934,868
0 0.500
0.14
A7.2
158 12,403 1,934,868
0 0.920
0.43
To be an equivalence relation, (xij )(i,j)∈N<2 must satisfy all the inequality constraints of the
problem (P ). Since minimizing pk=1 d(∼k , ∼X ) is equivalent to maximizing − qk=1 d(∼k
, ∼X ), we see that the problem of aggregation of equivalence relations reduces to the
problem (P ).
Table 1 shows the result of the MSMPT for the real-world instances in Grötschel andWakabayashi [8] where• the relative gap is the ratio of the duality gap to the best lower bound (i.e.,(ωup − ωlow )/ωlow ),• %! shows the percentage of the transitivity constraints considered in the MSMPT(i.e., 100(|U | + |V | + |W |)/|N<3 |),• %" shows the percentage of the pegged variables in the MSMPT (i.e., 100(|P0 | +|P1 |)/|N<2 |).
Optimal solution was obtained within 0.5 seconds for all the instances except instance A3.
From the column %!, we observe that the MSM works well. In addition, we observed thatFVD finds an optimal solution in Step 1 for all the instances except instance A1. Also we
see that only for instance A1, our Lagrangian heuristics (BLH) contributes to finding an
optimal solution. Note that the instance A6 and A7 have too huge number of inequality
constraints to apply the BSM on an ordinary workstation.
Figure 1 shows the comparison of our algorithm MSMPT and the cutting plane algorithm described in Grötschel and Wakabayashi [8](GW for short). We cannot say about
which algorithm is more efficient only from the computational times since the computational environment has changed dramatically. However, at least we can judge from thefigure that the computation time of the MSMPT is less sensitive to the increase of the
problem size n than that of the cutting plane algorithm GW. It should be pointed outthat for instance A3, the cutting plane algorithm GW finds an optimal solution while the
MSMPT terminates with a narrow duality gap remained.
7.2. Randomly generated instances. In this section, we report numerical experiments
for randomly generated instances. The instances are generated by the following method.
Step 1 Decide the problem size n and the number of equivalence relations q.
Step 2 Make a table D = (dik ) of dimension n × p by choosing dik randomly from theuniform distribution on {0, 1, 2}.
Step 3 Make equivalence relations ∼1 , ∼2 , . . . , ∼q by setting#
1 dik = djk ,krij:=
for all k ∈ {1, 2, . . . , q}.
0 otherwise,
Step 4 Set cij ← −

q
"
k=1

k
(1 − 2rij) for all (i, j) ∈ N<2 .

CLIQUE PARTITIONING PROBLEMS

19

Figure 1. Comparison of the MSMPT with the cutting plane algorithm
described in Grötschel and Wakabayashi [8] (GW for short)
Table 2. Result for the random instances (n = 30, 50, 100)
instance relative gap
R30a
4.E-02R30b
2.E-02R30c
2.E-02
R30d
0
R30e
0
R50a
3.E-02R50b0
R50c
1.E-01R50d
3.E-02R50e
3.E-02R100a
2.E-01R100b
2.E-01R100c
2.E-01R100d
2.E-01R100e
2.E-01

MSMPT%!
%" time(s)1.158 57.701
0.25
1.026 84.368
0.11
0.591 82.989
0.10
0.435 96.552
0.02
1.190 0.000
0.00
0.777 32.163
0.51
0.855 60.327
0.08
2.214 0.000
0.49
0.968 32.082
0.52
1.776 52.653
1.02
2.648 0.000
4.35
2.616 0.000
4.47
2.362 0.000
5.83
2.806 0.000
4.52
2.653 0.000
4.82

Xpressratio time(s)
1.000
1.01
1.000
1.03
1.000
1.01
1.000
1.03
1.000
1.01
1.000
4.76
1.000
4.90
1.000
4.76
1.000
4.77
1.000
4.77
1.000
43.96
1.000
46.32
0.996
47.55
1.000
45.94
1.000
43.93

For each n = 30, 50, 100, we generated five problems, all of which the number of equivalence relations q is fixed to 10. The result is summarized in Table 2. As mentioned in
Section 3, since one of the possible competitors of the MSM is the dual of the LP relaxation, here we also presented the result of the dual of the LP relaxation. To solve it, weused Xpress Optimizer 21.01.06 to solve them and run on a PC with an Intel i7, 2.80GHz
processor and 6 GB of memory. The column “ratio” shows the ratio of the best upperbound obtained by the MSMPT and the optimal value of the dual of the LP relaxation.
We observe that the MSMPT fails to find an optimal solution for almost all of therandom instances. In addition, the pegging test does not work very well especially forthe instances with n = 100. However, we see that the MSMPT excels the dual of the LP
relaxation in computational times and upper bounds. This shows that the MSMPT has
the potential to replace the linear programming relaxation in exact methods such as thebranch-and-bound method.

20

SUKEGAWA, YAMAMOTO, AND ZHANG

Table 3. Similarity of the instances
A1
A2
A3
A4
A6
A71
A72
- 0.765 0.738 0.780 0.696 0.907 0.861
The real-world instances would be different from randomly generated instances. To
examine the difference quantitatively, we define the degree of similarity of given equivalence
relations ∼1 , ∼2 , . . . , ∼q . Let sim(s, t) denote the cosine similarity of ∼s and ∼t with1 ≤ s < t ≤ q, i.e.,
"
s t
rij
rij
sim(s, t) := 6 "2(i,j)∈N<2(i,j)∈N<

s )2
(rij

6 "

t )2
(rij

,

2
(i,j)∈N<and simall denote the geometric mean of the cosine similarities, i.e.,
1simall := (sim(1, 2) × sim(1, 3) × · · · × sim(q − 1, q)) q(q−1)/2 .

Note that, in this case, 0 ≤ sim(s, t) ≤ 1 holds for all 1 ≤ s < t ≤ q. If the value simall isclose to 1, then one could say that the given equivalence relations are similar. The value
simall of the random instances we solved is about 0.575 on average. Table 3 shows thevalue simall of the instances in Grötschel and Wakabayashi [8]. Similarity is not available
for A1 since it has some missing values.
Judging from the similarity defined above, the real-world instances are likely to be easierthan the random instances. Keeping this point in mind, we generate random instances
which imitates the real-world instances in the following way.
Step 1 Decide the problem size n and the number of equivalence relations q.
Step 2 Make an n-dimensional vector (di1 ) by choosing di1 randomly from the uniform
distribution on {0, 1, 2} (this vector will be used as a reference equivalence relation).
Step 3 For each k ∈ {2, 3, . . . , q} and for each i ∈ {1, 2, . . . , n}, choose p randomly fromthe uniform distribution in [0, 1], and if p ≥ p̄ then set dik ← di0 , otherwise set
randomly from the uniform distribution in {0, 1, 2}.
Step 4 Make equivalence relations ∼1 , ∼2 , . . . , ∼q by setting#
1 dik = djk ,krij:=
for all k ∈ {1, 2, . . . , q}.
0 otherwise,
Step 5 Set cij ← −

q
"
k=1

k
(1 − 2rij) for all (i, j) ∈ N<2 .

The parameter p̄ controls the value simall . Namely, the smaller p̄ we set, the greater
the value simall will be. We set p̄ to 0.4 and generated five problems for each n = 100,
200, 300. The result is summarized in Table 4. Note that the instances with n = 300
have 44,850 decision variables and 13,365,300 inequality constraints. We observe that the
MSMPT finds an optimal solution for almost all these instances. This result shows that
the MSMPT has a potential to solve larger real-world instances. Even when the MSMPT
fails to find an optimal solution, we have problems of reduced size at hand. Then we could
solve them by an appropriate solver. In the next section, we will discuss the performanceof the MSMPT as a preprocessing for an exact method.

CLIQUE PARTITIONING PROBLEMS

21

Table 4. Result for the random instances with similarity (n = 100, 200, 300)
instance simall relative gap
RS100a0.65
0RS100b0.68
0RS100c0.66
0RS100d0.67
0RS100e0.66
0RS200a0.67
0RS200b0.68
0RS200c0.67
0RS200d0.68
0RS200e0.66
6.E-05
RS300a0.67
0RS300b0.67
0RS300c0.66
0RS300d0.67
0RS300e0.67
3.E-05%!
%" time (s)
1.125
0.000
0.17
1.199
9.354
0.50
1.193
0.000
0.06
1.372
0.000
0.49
1.145 98.101
0.93
1.164
7.367
4.65
1.090
0.000
5.52
1.228
0.000
6.85
1.087 99.030
16.80
1.118 99.467
140.34
1.143
0.000
40.47
1.133 100.000
159.571.221 98.370
240.66
1.171
0.000
53.48
1.160 99.802
469.64

Table 5. Comparison of (n# , m# ) and (|N<2 |, |N<3 |)instancen|N<2 ||N<3 | n# /|N<2 | m# /|N<3 |
A3
34
561
17,952 9.E-02
5.E-02
RS200e 200 19,900 3,940,200 5.E-03
9.E-04
RS300e 300 44,850 13,365,300 2.E-03
5.E-05
Table 6. Result of Xpress for the reduced problems
MSMPT
instance
ωlow
ωup
A3
1042.00 1043.00
RS200e 17238.00 17239.00
RS300e 36714.00 36715.00

Xpresstotal
optimal value time(s) time(s)1042.00
0.04
2.00
17238.00
0.35 140.69
36714.00
0.73 470.37

7.3. MSMPT as a preprocessing. Once some variables are pegged, some of the transitivity constraints become redundant. Take the transitivity constraint xij + xjk − xik ≤ 1of type U for instance. If xik is pegged at one, this inequality is redundant. It is also
redundant if xij or xjk is pegged at zero. If xij is pegged at one and xik is pegged atzero, the third variable xjk must take zero to meet the transitivity constraint. Then wecan peg it at zero and remove the constraint. Repeating this procedure, we can reduce
the size of the problem. Let n# and m# denote the number of the decision variables andthe transitivity constraints of the reduced problem, respectively. Table 5 shows how the
MSMPT reduces the size of the instances which we reported the MSMPT failed to solve
to optimality.
We observe that the reduction ratio of the transitivity constraints is better than that ofthe decision variables. To solve the reduced problems exactly, we used Xpress Optimizer
21.01.06 and run on a PC with an Intel i7, 2.80GHz processor and 6 GB of memory.
Table 6 shows the result.
Although the MSMPT terminates with a small duality gap to close up for instance
RS200e and RS300e, the best solutions found by the MSMPT turned out to be optimal

22

SUKEGAWA, YAMAMOTO, AND ZHANG

solutions. We observe that the computation time of Xpress is very short due to the
problem reduction. We conclude that the MSMPT deserves a good preprocessing forexact methods.
8. Conclusion
In order to manage the high dimensionality of the Lagrangian multipliers arising fromthe O(n3 ) number of constraints, we proposed in this paper a modification of the conventional subgradient method. This enables us to apply the Lagrangian relaxation method tolarge instances that exceed the capability of the conventional method. Since the method
reached the optimal value of the Lagrangian dual problem in almost all instances that we
tested, this modification is practical as a tool for finding a good upper bound of the cliquepartitioning problem.
Taking advantage of the structural property of the clique partitioning problem, we madean improvement on the ordinary pegging test, which was experimentally confirmed to peg
more variables than the ordinary one. We showed that a slightly complicated computationof the improved pegging test is performed efficiently.
Our numerical experiments showed that the algorithm works well for moderate-sized
instances of the aggregation problem of equivalence relations, and for even larger instances
when the equivalence relations are correlated.
References[1] S. G. de Amorim, J-P. Barthélemy and C. C. Riberio, “Clustering and clique partitioning: simulated
annealing and tabu search approaches”, Journal of Classification 9 (1992) 17–41.
[2] E. Balas and M. C. Carrera, “A dynamic subgradient-based branch-and-bound procedure for set covering”, Operations Research 44 (1996) 875–890.
[3] M. J. Brusco and H. F. Köhn, “Clustering qualitative data based on binary equivalence relations:neighborhood search heuristics for the clique partitioning problem”, PSYCHOMETRIKA 74 (2009)
685–703.
[4] I. Charon and O. Hudry, “Noising methods for a clique partitioning problem”, Discrete Applied Mathematics 154 (2006) 754–769.
[5] U. Dorndorf, F. Jaehn and E. Pesch, “Modeling robust flight-gate scheduling as a clique partitioning
problem”, Transportation Science 42 (2008) 292–301.
[6] M.L. Fisher, “The Lagrangian relaxation method for solving integer programming problems”, Management Science 27 (1981) 1–18.
[7] A.M. Geoffrion, “Lagrangian relaxation for integer programming”, Mathematical Programming Study
2 (1974) 82–114.
[8] M. Grötschel and Y. Wakabayashi, “A cutting plane algorithm for a clustering problem”, Mathematical
Programming 45 (1989) 59–96.
[9] G. Kochenberger, F. Glover, B. Alidaee and H. Wang, “Clustering of microarray data via clique partitioning”, Journal of Combinatorial Optimization 10 (2005) 77–92.
[10] T. Larsson, M. Patriksson and A.-B. Strömberg, “Conditional subgradient optimization - Theory andapplications”, European Journal of Operational Research 88 (1996) 382–403.
[11] A. Mehrotra and M. A. Trick, “Cliques and clustering: A combinatorial approach”, Operations Research Letters 22 (1998) 1–12
[12] R. M. Nauss “An efficient algorithm for the 0-1 knapsack problem”, Management Science 23 (1976)
27–31.
[13] M. Oosten, J. H. G. C. Rutten and F. C. R. Spieksma, “The clique partitioning problem: Facets andpatching facets”, Networks 38 (2001) 209–226.
[14] B. You and T. Yamada “A pegging approach to the precedence-constrained knapsack problem”,
European Journal of Operational Research 183 (2007) 618–632
Faculty of Engineering, Information and Systems, University of Tsukuba, Tsukuba, Ibaraki
305-8573, Japan
E-mail address: sukegawa.n.aa@m.titech.ac.jp, yamamoto@sk.tsukuba.ac.jp, zhangliyuan120@gmail.com