Exploiting variable associations to configure efficient
local search algorithms in large-scale binary integer
programsI

arXiv:1604.08448v2 [] 12 May 2017

Shunji Umetani∗
Graduate School of Information Science and Technology, Osaka University,
1-5 Yamadaoka, Suita, Osaka, 565-0871, Japan

Abstract
We present a data mining approach for reducing the search space of local search
algorithms in a class of binary integer programs including the set covering and
partitioning problems. The quality of locally optimal solutions typically improves if a larger neighborhood is used, while the computation time of searching
the neighborhood increases exponentially. To overcome this, we extract variable
associations from the instance to be solved in order to identify promising pairs of
flipping variables in the neighborhood search. Based on this, we develop a 4-flip
neighborhood local search algorithm that incorporates an efficient incremental
evaluation of solutions and an adaptive control of penalty weights. Computational results show that the proposed method improves the performance of the
local search algorithm for large-scale set covering and partitioning problems.
Keywords: combinatorial optimization, heuristics, set covering problem, set
partitioning problem, local search

1. Introduction
The Set Covering Problem (SCP) and Set Partitioning Problem (SPP) are
representative combinatorial optimization problems that have many real-world
applications, such as crew scheduling (Barnhart et al., 1998; Hoffman & Padberg, 1993; Mingozzi et al., 1999), vehicle routing (Agarwal et al., 1989; Baldacci
et al., 2008; Bramel & Simchi-Levi, 1997; Hashimoto et al., 2009), facility location (Boros et al., 2005; Farahani et al., 2012) and logical analysis of data
(Boros et al., 2000; Hammer & Bonates, 2006). Real-world applications of SCP
and SPP are comprehensively reviewed in (Ceria et al., 1997) and (Balas &
Padberg, 1976), respectively.
IA

preliminary version of this paper was presented in (Umetani, 2015).
author. Tel.: +81 (6) 6879 7799.
Email address: umetani@ist.osaka-u.ac.jp (Shunji Umetani)

∗ Corresponding

Preprint submitted to Elsevier

January 3, 2018

Given a ground set of m elements i ∈ M = {1, . . . , m}, n subsets Sj ⊆ M
(|Sj | ≥ 1), and their costs cj ∈ R (R is the set Sof real values) for j ∈ N =
{1, . . . , n}, we say that X ⊆ N is a cover
S of M if j∈X Sj = M holds. We also
say that X ⊆ N is a partition of M if j∈X Sj = M and Sj1 ∩ Sj2 = ∅ for all
j1 , j2 ∈ X hold. The goals of SCP and SPP are to find a minimum cost cover
and partition X of M , respectively. In this paper, we consider the following
class of Binary Integer Programs (BIPs) including SCP and SPP:
X
minimize
cj xj
subject to

j∈N
X
j∈N
X
j∈N
X

aij xj ≤ bi ,

i ∈ ML ,

aij xj ≥ bi ,

i ∈ MG ,

aij xj = bi ,

i ∈ ME ,

(1)

j∈N

xj ∈ {0, 1},

j ∈ N,

where aij ∈ {0, 1} and bi ∈ Z+ (Z+ is the set of nonnegative integer values).
We note that aij = 1 if i ∈ Sj holds and aij = 0 otherwise, and xj = 1 if j ∈ X
holds and xj = 0 otherwise. That is, a column vector aj = (a1j , . . . , amj )> of
the matrix (aij ) represents the corresponding subset Sj = {i ∈ M | aij = 1},
and the vector x also represents the corresponding cover (or partition) X =
{j ∈ N | xj = 1}. For notational convenience, we denote M = ML ∪ MG ∪ ME .
For each i ∈ M , let Ni = {j ∈ N | aij =P1} be the index set of subsets Sj that
contain the elements i, and let si (x) = j∈N aij xj be the left-hand side of the
ith constraint.
Continuous development of mathematical programming has much improved
the performance of exact and heuristic algorithms and this has been accompanied by advances in computing machinery. Many efficient exact and heuristic algorithms for large-scale SCP and SPP instances have been developed (Atamtürk
et al., 1995; Barahona & Anbil, 2000; Bastert et al., 2010; Borndörfer, 1998;
Boschetti et al., 2008; Caprara et al., 1999, 2000; Caserta, 2007; Ceria et al.,
1998; Linderoth et al., 2001; Umetani & Yagiura, 2007; Wedelin, 1995; Yagiura
et al., 2006), many of which are based on a variant of the column generation
method called the pricing method that reduces the number of variables to be
considered in the search by using Linear Programming (LP) and/or Lagrangian
relaxation. However, many large-scale SCP and SPP instances still remain unsolved because there is little hope of closing the large gap between the lower
and upper bounds of the optimal values. In particular, the equality constraints
of SPP often make the pricing method less effective because they often prevent
solutions from containing highly evaluated variables together. In this paper, we
consider an alternative approach for extracting useful features from the instance
to be solved with the aim of reducing the search space of local search algorithms
for large-scale SCP and SPP instances.
In the design of local search algorithms for combinatorial optimization prob2

lems, the quality of locally optimal solutions typically improves if a larger neighborhood is used. However, the computation time of searching the neighborhood
also increases exponentially. To overcome this, extensive research has investigated ways to efficiently implement neighborhood search, which can be broadly
classified into three types: (i) reducing the number of candidates in the neighborhood (Pesant & Gendreau, 1999; Shaw et al., 2002; Yagiura et al., 1999, 2006),
(ii) evaluating solutions by incremental computation (Michel & Van Hentenryck,
2000; Yagiura et al., 1999; Van Hentenryck & Michel, 2005; Voudouris et al.,
2001), and (iii) reducing the number of variables to be considered by using LP
and/or Lagrangian relaxation (Caprara et al., 1999; Ceria et al., 1998; Umetani
et al., 2013; Yagiura et al., 2006).
To suggest an alternative, we develop a data mining approach for reducing
the search space of local search algorithms. That is, we construct a k-nearest
neighbor graph by extracting variable associations from the instance to be solved
in order to identify promising pairs of flipping variables in the neighborhood
search. We also develop a 4-flip neighborhood local search algorithm that flips
four variables alternately along 4-paths or 4-cycles in the k-nearest neighbor
graph. We incorporate an efficient incremental evaluation of solutions and an
adaptive control of penalty weights into the 4-flip neighborhood local search
algorithm.
2. 2-flip neighborhood local search
Local Search (LS) starts from an initial solution x and then iteratively replaces x with a better solution x0 in the neighborhood NB(x) until no better
solution is found in NB(x). For some positive integer r, let the r-flip neighborhood NBr (x) be the set of solutions obtainable by flipping at most r variables in
x. We first develop a 2-Flip Neighborhood Local Search (2-FNLS) algorithm as
a basic component of our algorithm. In order to improve efficiency, the 2-FNLS
first searches NB1 (x), and then searches NB2 (x) \ NB1 (x) only if x is locally
optimal with respect to NB1 (x).
The BIP is NP-hard, and the (supposedly) simpler problem of judging the
existence of a feasible solution is NP-complete. We accordingly consider the
following formulation of the BIP that allows violations of the constraints and
introduce the following penalized objective function with penalty weights wi+ ∈
R+ (R+ is the set of nonnegative real values) for i ∈ ML ∪ ME and wi− ∈ R+

3

for i ∈ MG ∪ ME .
minimize

z(x) =

subject to

X

X
j∈N

j∈N
X
j∈N
X

X

cj x j +

wi+ yi+ +

X

wi− yi−

i∈MG ∪ME

i∈ML ∪ME

aij xj − yi+ ≤ bi ,

i ∈ ML ,

aij xj + yi− ≥ bi ,

i ∈ MG ,

aij xj − yi+ + yi− = bi ,

i ∈ ME ,

(2)

j∈N

xj ∈ {0, 1},
yi+ ≥ 0,
yi− ≥ 0,

j ∈ N,
i ∈ ML ∪ M E ,
i ∈ MG ∪ ME .

For a given x ∈ {0, 1}n , we can easily compute optimal yi+ = |si (x) − bi |+ and
yi− = |bi − si (x)|+ , where we denote |x|+ = max{x, 0}.
Because the region searched by a single application of LS is limited, LS
is usually applied many times. When a locally optimal solution is found, a
standard strategy is to update the penalty weights and to resume LS from the
obtained locally optimal solution. We accordingly evaluate solutions by using an
alternative function z̃(x), where the original penalty weight vectors w+ and w−
e + and w
e − , respectively, and these are adaptively controlled
are replaced by w
during the search (see Section 6 for details).
We first describe how 2-FNLS is used to search NB1 (x), which is called the
1-flip neighborhood. Let
(
∆z̃j↑ (x) j ∈ N \ X
∆z̃j (x) =
(3)
∆z̃j↓ (x) j ∈ X,
be the increase in z̃(x) due to flipping xj ← 1 − xj , where
X
X
∆z̃j↑ (x) = cj +
w
ei+ −
i∈Sj ∩(ML ∪ME )∩{l|sl (x)≥bl }

∆z̃j↓ (x)

= −cj −

X

i∈Sj ∩(MG ∪ME )∩{l|sl (x)<bl }

w
ei+

i∈Sj ∩(ML ∪ME )∩{l|sl (x)>bl }

w
ei− ,

+

X

w
ei− ,

i∈Sj ∩(MG ∪ME )∩{l|sl (x)≤bl }

(4)
are the increases in z̃(x) due to flipping xj = 0 → 1 and xj = 1 → 0, respectively. 2-FNLS first searches for an improved solution obtainable by flipping
xj ← 1 − xj for j ∈ N . If an improved solution exists, it chooses j with the
minimum value of ∆z̃j (x) and flips xj ← 1 − xj .
We next describe how 2-FNLS is used to search NB2 (x) \ NB1 (x), which is
called the 2-flip neighborhood. We derive conditions that reduce the number
of candidates in NB2 (x) \ NB1 (x) without sacrificing the solution quality by
expanding the results as shown in (Yagiura et al., 2006). Let ∆z̃j1 ,j2 (x) be the
increase in z̃(x) due to simultaneously flipping the values of xj1 and xj2 .

4

Lemma 1. If a solution x is locally optimal with respect to NB1 (x), then
∆z̃j1 ,j2 (x) < 0 holds only if Sj1 ∩ Sj2 6= ∅ and xj1 6= xj2 .
Proof. By the assumption of the lemma, ∆z̃j1 (x) ≥ 0 and ∆z̃j2 (x) ≥ 0 hold.
It is clear from (4) that ∆z̃j1 ,j2 = ∆z̃j1 (x) + ∆z̃j2 (x) ≥ 0 holds if Sj1 ∩ Sj2 = ∅.
We show that ∆z̃j1 ,j2 (x) ≥ 0 holds if Sj1 ∩ Sj2 6= ∅ and xj1 = xj2 . First,
we consider the case of xj1 = xj2 = 1. If si (x) = bi + 1 holds for i ∈ Sj1 ∩
Sj2 ∩ (ML ∪ ME ), then decrease of the violation yi+ partly cancels by flipping
xj1 = 1 → 0 and xj2 = 1 → 0 simultaneously. Similarly, if si (x) = bi + 1
holds for i ∈ Sj1 ∩ Sj2 ∩ (MG ∪ ME ), then a new violation yi− occurs by flipping
xj1 = 1 → 0 and xj2 = 1 → 0 simultaneously. We then have
∆z̃j1 ,j2 (x) = ∆z̃j↓1 (x) + ∆z̃j↓2 (x) +

X

w
ei+

i∈Sj1 ∩Sj2 ∩(ML ∪ME )∩{l|sl (x)=bl +1}

X

+

w
ei− ≥ 0. (5)

i∈Sj1 ∩Sj2 ∩(MG ∪ME )∩{l|sl (x)=bl +1}

Next, we consider the case of xj1 = xj2 = 0. If si (x) = bi − 1 holds
for i ∈ Sj1 ∩ Sj2 ∩ (ML ∪ ME ), then a new violation yi+ occurs by flipping
xj1 = 0 → 1 and xj2 = 0 → 1 simultaneously. Similarly, if si (x) = bi − 1 holds
for i ∈ Sj1 ∩ Sj2 ∩ (MG ∪ ME ), then decrease of the violation yi− partly cancels
by flipping xj1 = 0 → 1 and xj2 = 0 → 1 simultaneously. We then have
∆z̃j1 ,j2 (x) = ∆z̃j↑1 (x) + ∆z̃j↑2 (x) +

X

w
ei+

i∈Sj1 ∩Sj2 ∩(ML ∪ME )∩{l|sl (x)=bl −1}

X

+

w
ei− ≥ 0. (6)

i∈Sj1 ∩Sj2 ∩(MG ∪ME )∩{l|sl (x)=bl −1}


Based on this lemma, we consider only the case of xj1 = 1 and xj2 = 0. If
si (x) = bi holds for i ∈ Sj1 ∩ Sj2 ∩ (ML ∪ ME ), then increase of the violation yi+
by flipping xj2 = 0 → 1 cancels by flipping xj1 = 1 → 0 simultaneously, while
no decrease of the violation yi+ occurs by flipping xj1 = 1 → 0 independently.
Similarly, if si (x) = bi holds for i ∈ Sj1 ∩ Sj2 ∩ (MG ∪ ME ), then increase
of the violation yi− by flipping xj1 = 1 → 0 cancels by flipping xj2 = 0 → 1
simultaneously, while no decrease of the violation yi− occurs by flipping xj2 =
0 → 1 independently. We then have
X
X
w
ei+ −
w
ei− ,
∆z̃j1 ,j2 (x) = ∆z̃j↓1 (x) + ∆z̃j↑2 (x) −
i∈S̄(x)∩(ML ∪ME )

i∈S̄(x)∩(MG ∪ME )

(7)
where S̄(x) = {i ∈ Sj1 ∩ Sj2 | si (x) = bi }. From these results, the 2-flip
neighborhood can be restricted to the set of solutions satisfying xj1 6= xj2 and
S̄(x) 6= ∅. However, it might not be possible to search this set efficiently without
5

first extracting it. We thus construct a neighbor list that stores promising pairs
of variables xj1 and xj2 for efficiency (see Section 4 for details).
To increase the efficiency of 2-FNLS, we decompose the neighborhood NB2 (x)
(j )
into a number of sub-neighborhoods. Let NB2 1 (x) = {x0 ∈ NB2 (x) | xj1 =
1, x0j1 = 0} be the subset of NB2 (x) obtainable by flipping xj1 = 1 → 0. 2FNLS searches NB2 1 (x) for each j1 ∈ X in ascending order of ∆z̃j↓1 (x). If
an improved solution exists, it chooses the pair j1 and j2 with the minimum
(j )
value of ∆z̃j1 ,j2 (x) among NB2 1 (x) and flips xj1 = 1 → 0 and xj2 = 0 → 1.
2-FNLS skips to search the remaining sub-neighborhoods and immediately return to search NB1 (x) whenever an improved solution is obtained in a sub(j )
neighborhood NB2 1 (x). The first version of 2-FNLS is formally described as
Algorithm 1.
(j )

e+, w
e−)
Algorithm 1 2-FNLS(x, w
e + and w
e−.
Input: A solution x and penalty weight vectors w
Output: A solution x.
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:

START:
x0 ← x
for j ∈ N do
if z̃(x) + ∆z̃j (x) < z̃(x0 ) then
x0 ← x, x0j ← 1 − x0j
end if
end for
if z̃(x0 ) < z̃(x) then
x ← x0
goto START
end if
for j1 ∈ X in ascending order of ∆z̃j↓1 (x) do
x0 ← x
for j2 ∈ N \ X do
if z̃(x) + ∆z̃j1 ,j2 (x) < z̃(x0 ) then
x0 ← x, x0j1 ← 0, x0j2 ← 1
end if
end for
if z̃(x0 ) < z̃(x) then
x ← x0
goto START
end if
end for

{Search NB1 (x)}

{Search NB2 (x)}
(j )

{Search NB2 1 (x)}

3. Efficient incremental evaluation
The 2-FNLS requires O(σ) time to compute the value of the evaluation
function z̃(x) for the current solution x if implemented naively, where σ =

6

P

P

j∈N aij denote the number of nonzero elements in the constraint matrix (aij ). To overcome this, we first develop a standard incremental evaluation
of ∆z̃j↑ (x) and ∆z̃j↓ (x) in O(|Sj |) time by keeping the values of the left-hand
side of constraints si (x) for i ∈ M in memory. We further develop an improved
incremental evaluation of ∆z̃j↑ (x) and ∆z̃j↓ (x) in O(1) time by keeping additional auxiliary data in memory. By using this, 2-FNLS is also able to evaluate
∆z̃j1 ,j2 (x) in O(|Sj |) time by (7).
We first consider a standard incremental evaluation of ∆z̃j↑ (x) and ∆z̃j↓ (x)
in O(|Sj |) time using the following formulas:
i∈M

∆z̃j↑ (x)
∆p̃↑j (x)

= cj + ∆p̃↑j (x) + ∆q̃j↑ (x),
X
=
w
ei+ (|(si (x) + 1) − bi |+ − |si (x) − bi |+ ) ,
i∈Sj ∩(ML ∪ME )

∆q̃j↑ (x)

=

X

w
ei−

(8)

(|bi − (si (x) + 1)|+ − |bi − si (x)|+ ) ,

i∈Sj ∩(MG ∪ME )

∆z̃j↓ (x)
∆p̃↓j (x)

= −cj + ∆p̃↓j (x) + ∆q̃j↓ (x),
X
=
w
ei+ (|(si (x) − 1) − bi |+ − |si (x) − bi |+ ) ,
i∈Sj ∩(ML ∪ME )

∆q̃j↓ (x)

=

X

w
ei−

(9)

(|bi − (si (x) − 1)|+ − |bi − si (x)|+ ) ,

i∈Sj ∩(MG ∪ME )

where 2-FNLS keeps the values of the left-hand side of constraints si (x) for
i ∈ M in memory. 2-FNLS updates si (x) for i ∈ Sj in O(|Sj |) time by si (x0 ) ←
si (x) + 1 and si (x0 ) ← si (x) − 1 when the current solution x moves to x0 by
flipping xj = 0 → 1 and xj = 1 → 0, respectively.
We further develop an improved incremental evaluation of ∆z̃j↑ (x) and ∆z̃j↓ (x)
in O(1) time by directly keeping ∆p̃↑j (x), ∆q̃j↑ (x) for j ∈ N \ X and ∆p̃↓j (x),
∆q̃j↓ (x) for j ∈ X in memory. When the current solution x moves to x0 by
flipping xj = 0 → 1, 2-FNLS first updates si (x) for i ∈ Sj in O(|Sj |) time by
si (x0 ) ← si (x) + 1, and then updates ∆p̃↑k (x), ∆q̃k↑ (x) for k ∈ Ni \ X, i ∈ Sj
P
and ∆p̃↓k (x), ∆q̃k↓ (x) for k ∈ Ni ∩ X, i ∈ Sj in O( i∈Sj |Ni |) time using the
following formulas:
X

∆p̃↑k (x0 ) ← ∆p̃↑k (x) +
w
ei+ ∆yi+ (x0 ) − ∆yi+ (x) ,
i∈Sj ∩Sk ∩(ML ∪ME )

∆q̃k↑ (x0 )

← ∆q̃k↑ (x) +

X


w
ei− ∆yi− (x0 ) − ∆yi− (x) ,

i∈Sj ∩Sk ∩(MG ∪ME )

∆p̃↓k (x0 )

∆p̃↓k (x)

←

∆q̃k↓ (x0 )

← ∆q̃k↓ (x) +

+

X


w
ei+ ∆yi+ (x0 ) − ∆yi+ (x) ,

i∈Sj ∩Sk ∩(ML ∪ME )

X


w
ei− ∆yi− (x0 ) − ∆yi− (x) ,

i∈Sj ∩Sk ∩(MG ∪ME )

7

(10)





(aij ) = 



0
1
1
0
1
0

1
0
1
0
1
0

1
1
0
1
0
0

1
0
0
0
1
1

0
1
0
1
0
1

0
0
0
1
1
0

0
1
1
1
0
0

0
0
1
1
1
0

1
0
0
1
0
1

0
1
1
0
0
1





 −→



x1
x2
x3
x4
x5
x6
x7
x8

x9
x10

x2

x7

x8

x1

x4

x8

x5

x7

x9

x2

x9

x3

x7

x9

x10

x10

x8
x1

x3

x5

x8

x1

x2

x6

x7

x3

x4

x5

x1

x5

x7

x10

Figure 1: Example of the neighbor list

where

∆yi+ (x0 ) = |(si (x0 ) + 1) − bi |+ − |si (x0 ) − bi |+ ,
∆yi+ (x) = |si (x0 ) − bi |+ − |si (x) − bi |+ ,
∆yi− (x0 ) = |bi − (si (x0 ) + 1)|+ − |bi − si (x0 )|+ ,
∆yi− (x) = |bi − si (x0 )|+ − |bi − si (x)|+ .

(11)

Similarly, when the current solution x moves to x0 by flipping xj = 1 → 0, 2FNLS first updates si (x) for i ∈ Sj in O(|Sj |) time, and then updates ∆p̃↑k (x),
↓
↓
∆q̃k↑ (x)
P for k ∈ Ni \ X, i ∈ Sj and ∆p̃k (x), ∆q̃k (x) for k ∈ Ni ∩ X, i ∈ Sj
in O( i∈Sj |Ni |) time. We note that the computation time for updating the
auxiliary data has little effect on the total computation time of 2-FNLS because
the number of solutions actually visited is much less than the number of neighbor
solutions evaluated in most cases.
4. Exploiting variable associations
It is known that the quality of locally optimal solutions improves if a larger
neighborhood is used. However, the computation time to search the neighborhood NBr (x) also increases exponentially with r, since |NBr (x)| = O(nr )
generally holds. A large amount of computation time is thus needed in practice
in order to scan all candidates in NB2 (x) for large-scale instances with millions
of variables. To overcome this, we develop a data mining approach that identifies
promising pairs of flipping variables in NB2 (x) by extracting variable associations from the instance to be solved using only a small amount of computation
time.
From the results in Section 2, the 2-flip neighborhood can be restricted to
the set of solutions satisfying xj1 6= xj2 and S̄(x) 6= ∅. We further observe from
(7) that it is favorable to select pairs of flipping variables xj1 and xj2 which
gives a larger size |Sj1 ∩ Sj2 | in order to obtain ∆z̃j1 ,j2 (x) < 0. Based on this
observation, we keep a limited set of pairs of variables xj1 and xj2 for which
8

|Sj1 ∩ Sj2 | is large in memory, and we call this the neighbor list (Figure 1). We
note that |Sj1 ∩ Sj2 | represents a certain kind of similarity between the subsets
Sj1 and Sj2 (or column vectors aj1 and aj2 of the constraint matrix (aij )) and
we keep the k-nearest neighbors for each variable xj for j ∈ N in the neighbor
list. (In Figure 1, we keep the set of pairs of variables xj1 and xj2 having
|Sj1 ∩ Sj2 | ≥ 2 in the neighbor list, because small examples often have many
ties.)
For each variable xj1 for j1 ∈ N , we store the first k = min{|N (j1 ) |, α|M |}
variables xj2 (j2 6= j1 ) in descending order of |Sj1 ∩ Sj2 | in the j1 th row of the
neighbor list, where N (j1 ) = {j2 ∈ N | j2 6= j1 , Sj1 ∩Sj2 6= ∅} and α is a program
parameter that we set to five. Let L[j1 ] be the index set of variables xj2 stored
in the j1 th row of the neighbor list. We then reduce the number of candidates
in NB2 (x) by restricting the pairs of flipping variables xj1 and xj2 to pairs in
the neighbor list j1 ∈ X and j2 ∈ (N \ X) ∩ L[j1 ].
We note that it is still expensive to construct the whole neighbor list for
large-scale instances with millions of variables. To overcome this, we develop a
lazy construction algorithm for the neighbor list. That is, 2-FNLS starts from
an empty neighbor list and generates the j1 th row of the neighbor list L[j1 ] only
(j )
when 2-FNLS searches NB2 1 (x) for the first time. The improved version of
2-FNLS is formally described as Algorithm 2.
A similar approach has been developed in local search algorithms for the
Euclidean Traveling Salesman Problem (TSP) in which a sorted list containing
only the k-nearest neighbors is stored for each city by using a geometric data
structure called the k-dimensional tree (Johnson & McGeoch, 1997). However, it
is not suitable for finding the k-nearest neighbors efficiently in high-dimensional
spaces. We thus extend it for application to the high-dimensional column vectors
aj ∈ {0, 1}m for j ∈ N of BIPs by using a lazy construction algorithm for the
neighbor list.
5. 4-flip neighborhood local search
We can regard the neighbor-list in Section 4 as an adjacency-list representation of a directed graph, and represent associations between variables by a
corresponding directed graph called the k-nearest neighbor graph. Figure 2
illustrates an example of the k-nearest neighbor graph corresponding to the
neighbor-list in Figure 1.
Using the k-nearest neighbor graph, we extend 2-FNLS to search a set of
promising neighbor solutions in NB4 (x). For each variable xj1 for j1 ∈ X, we
keep the variables xj2 for j2 ∈ (N \ X) ∩ L[j1 ] having the minimum value of
∆z̃j1 ,j2 (x) in memory as j2 = π(j1 ). The extended 2-FNLS, called the 4-Flip
Neighborhood Local Search (4-FNLS) algorithm, then searches for an improved
solution by flipping xj1 = 1 → 0, xπ(j1 ) = 0 → 1, xj3 = 1 → 0 and xπ(j3 ) =
0 → 1 for j1 ∈ X and j3 ∈ X ∩ L[π(j1 )] satisfying j1 6= j3 and π(j1 ) 6= π(j3 ),
i.e., flipping the values of variables alternately along 4-paths or 4-cycles in the
k-nearest neighbor graph. Let ∆z̃j1 ,j2 ,j3 ,j4 (x) be the increase in z̃(x) due to

9

e+, w
e−)
Algorithm 2 2-FNLS(x, w
e + and w
e−.
Input: A solution x and penalty weight vectors w
Output: A solution x.
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:

for j1 ∈ N do
Set L[j1 ] ← ∅
end for
START:
x0 ← x
for j ∈ N do
if z̃(x) + ∆z̃j (x) < z̃(x0 ) then
x0 ← x, x0j ← 1 − x0j
end if
end for
if z̃(x0 ) < z̃(x) then
x ← x0
goto START
end if
for j1 ∈ X in ascending order of ∆z̃j↓1 (x) do
if L[j1 ] = ∅ then
Generate L[j1 ]
end if
x0 ← x
for j2 ∈ (N \ X) ∩ L[j1 ] do
if z̃(x) + ∆z̃j1 ,j2 (x) < z̃(x0 ) then
x0 ← x, x0j1 ← 0, x0j2 ← 1
end if
end for
if z̃(x0 ) < z̃(x) then
x ← x0
goto START
end if
end for

10

{Initialize L[j1 ]}

{Search NB1 (x)}

{Search NB2 (x)}
{Generate L[j1 ]}

(j )

{Search NB2 1 (x)}

x5
x7
x10

x3

x8

x6

x9

x2
x4
x1

Figure 2: Example of the k-nearest neighbor graph

simultaneously flipping xj1 = 1 → 0, xj2 = 0 → 1, xj3 = 1 → 0 and xj4 = 0 →
1. 4-FNLS computes ∆z̃j1 ,j2 ,j3 ,j4 (x) in O(|Sj |) time by applying the standard
incremental evaluation alternately (see Section 3 for details). 4-FNLS is formally
described with as Algorithm 3.
6. Adaptive control of penalty weights
In our algorithm, solutions are evaluated by the alternative evaluation function z̃(x) in which the fixed penalty weight vectors w+ and w− in the original
e + and w
e − , respectively, and
evaluation function z(x) has been replaced by w
+
−
the values of w
ei for i ∈ ML ∪ ME and w
ei for i ∈ MG ∪ ME are adaptively
controlled in the search.
It is often reported that a single application of LS tends to stop at a locally
optimal solution of insufficient quality when large penalty weights are used.
This is because it is often unavoidable to temporally increase the values of
some violations yi+ and yi− in order to reach an even better solution from a
good solution through a sequence of neighborhood operations, and large penalty
weights thus prevent LS from moving between such solutions. To overcome this,
we incorporate an adaptive adjustment mechanism for determining appropriate
values of penalty weights w
ei+ for i ∈ ML ∪ ME and w
ei− for i ∈ MG ∪ ME
(Nonobe & Ibaraki, 2001; Yagiura et al., 2006; Umetani et al., 2013). That is,
LS is applied iteratively while updating the values of the penalty weights w
ei+
−
for i ∈ ML ∪ ME and w
ei for i ∈ MG ∪ ME after each call to LS. We call this
sequence of calls to LS the Weighting Local Search (WLS) according to (Selman
& Kautz, 1993; Thornton, 2005). This strategy is also referred as the breakout
algorithm (Morris, 1993) and the dynamic local search (Hutter et al., 2002) in
the literature.
Let x be the solution at which the previous local search stops. WLS ree + and w
e − . We
sumes LS from x after updating the penalty weight vectors w
+
assume that the original penalty weights wi for i ∈ ML ∪ ME and wi− for

11

e+, w
e−)
Algorithm 3 4-FNLS(x, w
e + and w
e−.
Input: A solution x and penalty weight vectors w
Output: A solution x.
1:
2:
3:
4:
5:
6:
7:
8:
9:
10:
11:
12:
13:
14:
15:
16:
17:
18:
19:
20:
21:
22:
23:
24:
25:
26:
27:
28:
29:
30:
31:
32:
33:
34:
35:
36:
37:
38:
39:
40:
41:
42:
43:
44:

for j1 ∈ N do
{Initialize L[j1 ]}
Set L[j1 ] ← ∅
end for
START:
x0 ← x
for j ∈ N do
{Search NB1 (x)}
if z̃(x) + ∆z̃j (x) < z̃(x0 ) then
x0 ← x, x0j ← 1 − x0j
end if
end for
if z̃(x0 ) < z̃(x) then
x ← x0
goto START
end if
for j1 ∈ X in ascending order of ∆z̃j↓1 (x) do
{Search NB2 (x)}
if L[j1 ] = ∅ then
{Generate L[j1 ]}
Generate L[j1 ]
end if
x0 ← x, ∆z̃ ∗ ← ∞
(j )
for j2 ∈ (N \ X) ∩ L[j1 ] do
{Search NB2 1 (x)}
∗
if ∆z̃j1 ,j2 (x) < ∆z̃ then
{Update π(j1 )}
∆z̃ ∗ ← ∆z̃j1 ,j2 (x), π(j1 ) ← j2
end if
if z̃(x) + ∆z̃j1 ,j2 (x) < z̃(x0 ) then
x0 ← x, x0j1 ← 0, x0j2 ← 1
end if
end for
if z̃(x0 ) < z̃(x) then
x ← x0
goto START
end if
end for
for j1 ∈ X in ascending order of ∆z̃j1 ,π(j1 ) (x) do
{Search NB4 (x)}
x0 ← x
for j3 ∈ X ∩ L[π(j1 )] satisfying j3 6= j1 and π(j1 ) 6= π(j3 ) do
if z̃(x) + ∆z̃j1 ,π(j1 ),j3 ,π(j3 ) (x) < z̃(x0 ) then
x0 ← x, x0j1 ← 0, x0π(j1 ) ← 1, x0j3 ← 0, x0π(j3 ) ← 1
end if
end for
if z̃(x0 ) < z̃(x) then
x ← x0
goto START
end if
end for

12

i ∈ MG ∪ME are sufficiently large. Starting from the original penalty weight vece+, w
e − ) ← (w+ , w− ), the penalty weight vectors w
e + and w
e − are updated
tors (w
∗
as follows. Let x be the best feasible solution obtained so far for the original
formulation (1). We modify 4-FNLS to keep the best feasible solution x∗ , and
update it whenever an improved feasible solution is found. If z̃(x) ≥ z(x∗ ) holds,
e+, w
e − ) ← β(w
e+, w
e − ),
WLS uniformly decreases the penalty weights by (w
where 0 < β < 1 is a program parameter that is adaptively computed so that
the new value of ∆z̃j↓ (x) becomes negative for 10% of variables xj for j ∈ X.
Otherwise, WLS increases the penalty weights by
z(x∗ ) − z̃(x)

yi+ , i ∈ ML ∪ ME ,
+2
−2
(y
+
y
)
l∈M l
l
z(x∗ ) − z̃(x)
w
ei− ← w
ei− + P
yi− , i ∈ MG ∪ ME .
+2
−2
(y
+
y
)
l∈M l
l

w
ei+ ← w
ei+ + P

(12)

e + and w
e−
WLS iteratively applies LS, updating the penalty weight vectors w
after each call to LS until the time limit is reached. WLS is formally described
as Algorithm 4. Note that we set the initial solution to x = 0 in practice.
Algorithm 4 WLS(x)
Input: An initial solution x.
Output: The best feasible solution x∗ .
e+, w
e − ) ← (w+ , w− )
1: x̃ ← x, z ∗ ← ∞, (w
2: repeat
e+, w
e−)
3:
(x̃, x0 ) ← 4-FNLS(x̃, w
0
e+, w
e − )}
4:
{x is the best feasible solution obtained in 4-FNLS(x̃, w
5:
if z(x0 ) < z ∗ then
6:
x∗ ← x0 , z ∗ ← z(x0 )
7:
end if
8:
if z̃(x̃) ≥ z ∗ then
e+, w
e − ) ← β(w
e+, w
e−)
9:
Decrease the penalty weights by (w
10:
else
11:
Increase the penalty weights by (12)
12:
end if
13: until The time limit is reached

7. Computational results
We report computational results of our algorithm for the SCP instances from
(Beasley, 1990; Umetani et al., 2013) and the SPP instances from (Borndörfer,
1998; Chu & Beasley, 1998; Koch et al., 2011). Tables 1 and 2 summarize the
information about the SCP and SPP instances, respectively. The first column
shows the name of the group (or the instance), and the number in parentheses
shows the number of instances in the group. In the subsequent part of this paper,
we show the average value in each cell of the tables for the instances of the group.
13

Table 1: Benchmark instances for SCP
instance
?G.1–5 (5)
?H.1–5 (5)
?I.1–5 (5)
?J.1–5 (5)
?K.1–5 (5)
?L.1–5 (5)
?M.1–5 (5)
?N.1–5 (5)
RAIL507
RAIL516
RAIL582
RAIL2536
?RAIL2586
?RAIL4284
?RAIL4872

zLP
149.48
45.67
138.97
104.78
276.67
209.34
415.78
348.93
172.15
182.00
209.71
688.40
935.92
1054.05
1509.64

zbest
166.4
59.6
158.0
129.0
313.2
258.0
549.8
503.8
∗174
∗182
∗211
∗689
947
1064
1530

original
#cst.
#var.
1000.0
10000.0
1000.0
10000.0
1000.0
50000.0
1000.0
100000.0
2000.0
100000.0
2000.0
200000.0
5000.0
500000.0
5000.0
1000000.0
507
63009
516
47311
582
55515
2536
1081841
2586
920683
4284
1092610
4872
968672

presolved
#cst.
#var.
1000.0
10000.0
1000.0
10000.0
1000.0
49981.0
1000.0
99944.8
2000.0
99971.0
2000.0
199927.6
5000.0
499988.0
5000.0
999993.2
440
20700
403
37832
544
27427
2001
480597
2239
408724
3633
607884
4207
482500

time limit
600 s
600 s
1200 s
1200 s
1800 s
1800 s
3600 s
3600 s
600 s
600 s
600 s
3600 s
3600 s
3600 s
3600 s

Table 2: Benchmark instances for SPP
instance
aa01–06 (6)
us01–04 (4)
v0415–0421 (7)
v1616–1622 (7)
t0415–0421 (7)
?t1716–1722 (7)
?ds
?ds-big
?ivu06-big
?ivu59

zLP
40372.75
9749.44
2385764.17
1021288.76
5199083.74
121445.76
57.23
86.82
135.43
884.46

zbest
∗40588.83
∗9798.25
∗2393303.71
∗1025552.43
5453475.71
157516.29
187.47
731.69
166.02
1878.83

original
#cst.
#var.
675.3
7587.3
121.3
295085.0
1479.3
30341.6
1375.7
83986.7
1479.3
7304.3
475.7
58981.3
656
67732
1042
174997
1177
2277736
3436
2569996

presolved
#cst.
#var.
478.7
6092.7
65.5
85772.5
263.9
7277.0
1171.9
51136.7
820.7
2617.4
475.7
13193.6
656
67730
1042
173026
1177
2197774
3413
2565083

time limit
600 s
600 s
600 s
600 s
600 s
3600 s
3600 s
3600 s
3600 s
3600 s

The detailed computational results are in the online supplement. The second
column “zLP ” shows the optimal values of the LP relaxation problems. The third
column “zbest ” shows the best upper bounds among all algorithms in this paper.
The fourth and sixth columns “#cst.” show the number of constraints, and the
fifth and seventh columns “#vars.” show the number of variables. Since several
preprocessing techniques that often reduce the size of instances by removing
redundant rows and columns are known (Borndörfer, 1998), all algorithms are
tested on the presolved instances. The instances marked “?” are hard instances
that cannot be solved optimally within at least 1 h by the tested Mixed Integer
Programming (MIP) solvers (i.e., CPLEX12.6, Gurobi5.6.3 and SCIP3.1).
We compare the results of our algorithm with those of its variations. Our
algorithm was implemented in C language and tested on a MacBook Pro laptop
computer with a 2.7 GHz Intel Core i7 processor and 16 GB memory. All
variations of our algorithm were run on a single thread under MacOS10.12
operating system with time limits as shown in Tables 1 and 2.
best
Tables 3 and 4 show the relative gap z(x)−z
× 100 (%) of the best feasible
z(x)
solutions for the original formulation (1) achieved by the variations of our algorithm for the SCP and SPP instances, respectively. We note that all variations
of our algorithm found feasible solutions for all SCP instances. The second
column “no-list” shows the results of our algorithm without the neighbor list,
and the third column “no-inc” shows the results of our algorithm without the
improved incremental evaluation (i.e., only applying the standard incremental
evaluation in Section 3). The fourth column “2-FNLS” shows the results of our
algorithm without the 4-flip neighborhood search (i.e., only applying 2-FNLS in
Section 4). In Table 4, the number in parentheses shows the number of instances
14

Table 3: Computational results of variations of the proposed algorithm for SCP instances
instance
?G.1–5 (5)
?H.1–5 (5)
?I.1–5 (5)
?J.1–5 (5)
?K.1–5 (5)
?L.1–5 (5)
?M.1–5 (5)
?N.1–5 (5)
RAIL507
RAIL516
RAIL582
RAIL2536
?RAIL2586
?RAIL4284
?RAIL4872
avg. (all)
avg. (with stars)

no-list
0.00%
0.31%
1.24%
2.42%
2.12%
3.44%
10.97%
19.11%
0.00%
0.00%
0.47%
2.68%
2.57%
5.42%
4.43%
4.55%
4.89%

no-inc
0.12%
0.31%
0.86%
1.67%
1.69%
3.51%
8.33%
22.06%
0.57%
0.00%
0.47%
2.27%
2.87%
5.17%
3.47%
4.42%
4.75%

2-FNLS
0.00%
0.31%
0.50%
1.68%
1.32%
2.35%
2.79%
4.76%
0.00%
0.00%
0.47%
1.29%
2.27%
2.74%
2.36%
1.65%
1.77%

proposed
0.00%
0.00%
0.50%
1.53%
1.26%
2.05%
2.65%
5.47%
0.00%
0.00%
0.00%
0.72%
1.56%
2.12%
1.80%
1.56%
1.69%

Table 4: Computational results of variations of the proposed algorithm for SPP instances
instance
aa01–06 (6)
us01–04 (4)
v0415–0421 (7)
v1616–1622 (7)
t0415–0421 (7)
?t1716–1722 (7)
?ds
?ds-big
?ivu06-big
?ivu59
avg. (all)
avg. (with stars)

no-list
2.33%(6)
0.04%(4)
0.00%(7)
0.62%(7)
1.46%(5)
5.46%(7)
36.03%
29.11%
5.31%
15.75%
3.76%(40/42)
10.37%(12/13)

no-inc
2.26%(6)
1.16%(4)
0.00%(7)
0.17%(7)
1.30%(6)
4.33%(7)
33.80%
0.00%
3.83%
11.39%
2.60%(41/42)
6.61%(12/13)

2-FNLS
2.07%(6)
0.63%(4)
0.00%(7)
0.09%(7)
0.29%(7)
5.71%(7)
24.13%
40.75%
2.25%
16.01%
3.35%(42/42)
9.48%(13/13)

proposed
1.60%(6)
0.04%(4)
0.00%(7)
0.09%(7)
0.92%(6)
2.45%(7)
0.00%
0.00%
0.00%
0.00%
0.81%(41/42)
1.43%(12/13)

for which the algorithm obtained at least one feasible solution within the time
limit, and the relative gap shows the average value for them.
Tables 5 and 6 show the computational efficiency of variations of our algorithm with respect to the number of calls to 4-FNLS (and 2-FNLS in the fourth
column “2-FNLS”), where the bottom rows show average factors normalized so
that that of our algorithm is set to one.
From these results, we observe that the proposed method improves the computational efficiency of the local search algorithm for the SCP and SPP instances
in comparison with its variations, and our algorithm attains good performance
even when the size of the neighbor list is considerably small. We also observe
that the 4-flip neighborhood search substantially improves the performance of
our algorithm even though there are fewer calls to 4-FNLS compared to 2-FNLS.
rows
Tables 7 and 8 show the ratio #generated
× 100 (%) of generated rows
|N |
in the neighbor list. We observe that our algorithm achieves good performance
while generating only a small part of the neighbor list for the large-scale instances.
We compare the result of our algorithm with those of general purpose solvers,
i.e., the latest MIP solvers called CPLEX12.6 (2014), Gurobi5.6.3 (2014) and
SCIP3.1 (Achterberg, 2009) and a local search solver for BIP (including nonlinear constraints and objective functions) called LocalSolver3.1 (Benoist et al.,
2011). LocalSolver3.1 is a simulated annealing based on ejection chain moves
specialized for maintaining the feasibility of Boolean constraints and an efficient
incremental evaluation using a directed acyclic graph. LocalSolver3.1 is not the

15

Table 5: The number of calls to 4-FNLS of variations of the proposed algorithm for SCP
instances
instance
?G.1–5 (5)
?H.1–5 (5)
?I.1–5 (5)
?J.1–5 (5)
?K.1–5 (5)
?L.1–5 (5)
?M.1–5 (5)
?N.1–5 (5)
RAIL507
RAIL516
RAIL582
RAIL2536
?RAIL2586
?RAIL4284
?RAIL4872
avg. factor

no-list
1528.2
1040.8
845.2
494.0
364.0
262.4
103.0
78.8
4498
2004
3079
429
416
182
197
0.37

no-inc
1094.2
434.4
750.6
454.6
444.8
273.8
117.2
66.0
6482
3345
4037
359
399
195
229
0.30

2-FNLS
5755.8
2191.6
5072.0
2580.4
2828.8
1862.6
1412.8
929.4
38195
16915
22123
4045
5398
2665
4214
2.21

proposed
3001.4
1248.2
2189.8
1057.8
1445.8
809.0
535.0
290.8
23692
12350
15140
3206
3028
1683
2374
1.00

Table 6: The number of calls to 4-FNLS of variations of the proposed algorithm for SPP
instances
instance
aa01–06 (6)
us01–04 (4)
v0415–0421 (7)
v1616–1622 (7)
t0415–0421 (7)
?t1716–1722 (7)
?ds
?ds-big
?ivu06-big
?ivu59
avg. factor

no-list
10152.2
35418.5
1516414.7
1197.7
9412.6
32689.3
3109
852
236
141
0.20

no-inc
20544.8
36491.75
1783742.3
4029.1
19885.4
60992.9
3463
1019
305
208
0.38

2-FNLS
93694.3
92397.5
4274501.6
17626.9
145916.4
281283.4
14142
4370
1293
620
1.70

proposed
56584.8
86031.0
2749825.9
10048.7
53751.7
174405.1
11414
3294
1060
438
1.00

latest version, but it performs better than the latest version (LocalSolver4.5)
for the SCP and SPP instances.
We note the following issues in comparing the performance of our algorithm
with that of the latest MIP solvers. Due to the lack of pruning mechanism, it
is inherently difficult to find optimal solutions by local search algorithms even
for instances having a small gap between the lower and upper bounds of the
optimal values, while the latest MIP solvers find optimal solutions quickly by
the branch-and-cut procedure. Of course, local search algorithms often obtain
good upper bounds close to the optimal values for the instances. On the other
hand, the latest MIP solvers often prefer running primal heuristics rather than
the branch-and-cut procedure for instances having a large gap between the lower
and upper bounds of the optimal values. Indeed, the latest MIP solvers include
dozens of primal heuristics and spend much computation time on finding good
feasible solutions (Lodi, 2013), e.g., SCIP3.1 reported that it spent 1644.73 s
out of 3600 s running primal heuristics for solving the “ds” instance, while it
spent 837.47 s for solving LP relaxation problems.
We also compare our algorithm with a 3-flip local search algorithm specially
tailored for SCP developed by Yagiura et al. (2006) (denoted by YKI). In order
to handle large-scale SCP instances, many heuristic algorithms have introduced
the pricing method that reduces the number of variables to be considered by
using LP and/or Lagrangian relaxation. YKI introduced a pricing method based
on Lagrangian relaxation that substantially reduces the number of variables to
be considered to 1.05% from the original SCP instances on average. Hence, in

16

Table 7: The ratio of generated rows in the neighbor list for SCP instances
instance
?G.1–5 (5)
?H.1–5 (5)
?I.1–5 (5)
?J.1–5 (5)
?K.1–5 (5)
?L.1–5 (5)
?M.1–5 (5)
?N.1–5 (5)
RAIL507
RAIL516
RAIL582
RAIL2536
?RAIL2586
?RAIL4284
?RAIL4872

1 min
3.58%
2.12%
1.39%
0.82%
1.20%
0.57%
0.20%
0.01%
13.41%
4.73%
8.19%
0.30%
0.44%
0.20%
0.37%

10 min
3.95%
2.43%
1.63%
1.10%
1.49%
0.98%
0.58%
0.21%
22.02%
8.94%
10.66%
1.36%
1.67%
1.01%
1.56%

20 min

30 min

1 h

1.71%
1.16%
1.57%
1.06%
0.73%
0.29%

1.61%
1.10%
0.81%
0.35%

0.93%
0.49%

1.87%
2.19%
1.43%
2.11%

2.17%
2.54%
1.73%
2.46%

2.73%
3.28%
2.29%
3.10%

Table 8: The ratio of generated rows in the neighbor list for SPP instances
instance
aa01–06 (6)
us01–04 (4)
v0415–0421 (7)
v1616–1622 (7)
t0415–0421 (7)
?t1716–1722 (7)
?ds
?ds-big
?ivu06-big
?ivu59

1 min
40.47%
3.93%
31.44%
6.55%
83.47%
61.00%
2.29%
0.21%
0.01%
0.01%

10 min
49.17%
5.17%
31.64%
7.42%
90.00%
94.38%
12.63%
2.06%
0.07%
0.05%

30 min

1 h

97.12%
27.11%
4.96%
0.23%
0.11%

97.98%
40.05%
8.10%
0.45%
0.16%

addition to the original SCP instances, we also tested all algorithms for reduced
SCP instances by another pricing method based on LP relaxation (Umetani
& Yagiura, 2007). Table 9 summarize the reduced SCP instances, where we
applied preprocessing to the reduced SCP instances as well as the original SCP
instances. The eighth column “#free var.” shows the number of variables to
be considered in YKI for the original SCP instances. We note that it turned
out that many of reduced SPP instances by the pricing method (Umetani &
Yagiura, 2007) were infeasible because the equality constraints of SPP often
prevent solutions from containing highly evaluated variables together and make
the pricing method less effective.
Tables 10 and 11 show the relative gap of the best feasible solutions for
the original formulation (1) achieved by the algorithms for the original and
reduced SCP instances, respectively. We note that all algorithms found feasible
solutions for all SCP instances. All algorithms were tested on a MacBook Pro
laptop computer with a 2.7 GHz Intel Core i7 processor and 16 GB memory,
and were run on a single thread under MacOS10.12 operating system with time
limits as shown in Tables 1 and 9. That is, we tested all algorithms under the
same amount of available computational resources for fair comparison.
We observe that our algorithm achieves better upper bounds than general
purpose solvers (CPLEX12.6.3, Gurobi5.6.3, SCIP3.1 and LocalSolver3.1) for
the original and reduced SCP instances. We also observe that our algorithm
achieves good upper bounds close to those of YKI for the reduced SCP instances.
It shows that our algorithm attains good performance close to that of the specially tailored algorithm (YKI) for SCP by combining a pricing method. In view
of these, our algorithm achieves sufficiently good upper bounds compared to the
other algorithms for the SCP instances.
17

Table 9: Reduced benchmark instances for SCP
instance
?G.1–5 (5)
?H.1–5 (5)
?I.1–5 (5)
?J.1–5 (5)
?K.1–5 (5)
?L.1–5 (5)
?M.1–5 (5)
?N.1–5 (5)
RAIL507
RAIL516
RAIL582
RAIL2536
?RAIL2586
?RAIL4284
?RAIL4872

zLP
149.48
45.67
138.97
104.78
276.67
209.34
415.78
348.93
172.15
183.00
209.71
688.68
935.92
1054.05
1510.87

zbest
166.4
59.8
158.4
129.4
313.8
259.0
551.6
505.0
∗174
∗183
∗211
∗691
948
1066
1532

reduced
#cst.
#var.
1000.0
441.8
1000.0
236.2
1000.0
721.6
1000.0
703.6
2000.0
1434.2
2000.0
1421.0
5000.0
3245.2
5000.0
4471.4
507
2649
516
3788
582
2091
2536
13746
2586
13349
4284
21728
4872
21329

presolved
#cst.
#var.
1000.0
441.8
1000.0
236.2
1000.0
721.6
1000.0
703.6
2000.0
1434.2
2000.0
1421.0
5000.0
3245.2
5000.0
4471.4
402
1019
350
3088
491
1493
1391
5782
2083
7377
3189
14565
3577
11404

Yagiura et al.
#free var.
339.0
174.6
479.4
429.6
959.4
856.4
1836.4
1660.6
394
492
513
1598
2089
2639
3650

time limit
600 s
600 s
1200 s
1200 s
1800 s
1800 s
3600 s
3600 s
600 s
600 s
600 s
3600 s
3600 s
3600 s
3600 s

Table 10: Computational results of the latest solvers and the proposed algorithm for SCP
instances
instance
?G.1–5 (5)
?H.1–5 (5)
?I.1–5 (5)
?J.1–5 (5)
?K.1–5 (5)
?L.1–5 (5)
?M.1–5 (5)
?N.1–5 (5)
RAIL507
RAIL516
RAIL582
RAIL2536
?RAIL2586
?RAIL4284
?RAIL4872
avg. (all)
avg. (with stars)

CPLEX12.6
0.37%
1.92%
2.81%
8.37%
4.77%
9.57%
18.43%
33.13%
0.00%
0.00%
0.00%
0.00%
2.27%
5.34%
1.73%
8.64%
9.45%

Gurobi5.6.3
0.49%
2.28%
2.72%
4.30%
4.38%
8.44%
10.10%
12.49%
0.00%
0.00%
0.00%
0.00%
2.17%
1.57%
1.73%
4.92%
5.38%

SCIP3.1
0.24%
1.93%
1.85%
3.59%
2.55%
3.52%
30.71%
42.32%
0.00%
0.00%
0.00%
0.86%
2.27%
30.55%
2.67%
10.00%
10.91%

LocalSolver3.1
45.80%
61.54%
41.38%
58.40%
51.22%
57.79%
71.08%
75.63%
5.43%
3.19%
5.80%
3.50%
5.39%
6.50%
5.61%
49.99%
54.22%

Yagiura et al.
0.00%
0.00%
0.00%
0.00%
0.00%
0.00%
0.00%
0.00%
0.00%
0.00%
0.00%
0.29%
0.00%
0.00%
0.00%
0.01%
0.00%

proposed
0.00%
0.00%
0.50%
1.53%
1.26%
2.05%
2.65%
5.47%
0.00%
0.00%
0.00%
0.72%
1.56%
2.12%
1.80%
1.56%
1.69%

Table 12 shows the relative gap of the best feasible solutions achieved by the
algorithms for the SPP instances. All algorithms were tested on a MacBook Pro
laptop computer with a 2.7 GHz Intel Core i7 processor and 16 GB memory,
and were run on a single thread under MacOS 10.12 operating system with time
limits as shown in Table 2. In Table 12, the number in parentheses shows that of
instances for which the algorithm obtained at least one feasible solution within
the time limit, and the relative gap shows the average value for them.
We also compare the computational results of our algorithm with those of a
Lagrangian heuristic algorithm for BIP called the Wedelin’s heuristic (Bastert
et al., 2010; Wedelin, 1995) and a branch-and-cut algorithm specially tailored
for SPP developed by Borndörfer (1998). The computational results for the
Wedelin’s heuristic are taken from (Bastert et al., 2010), where it was tested on
a 1.3 GHz Sun Ultra Sparc-IIIi and run with a time limit of 600 s. The computational results for the branch-and-cut algorithm are taken from (Borndörfer,
1998), where it was tested on a Sun Ultra Sparc 1 Model 170E and run with a
time limit of 7200 s.
We first observe that our algorithm achieves better upper bounds than the
latest MIP solvers (CPLEX12.6, Gurobi5.6.3 and SCIP3.1) for hard SPP instances and good upper bounds close to the optimal values for the others. We
next observe that our algorithm achieves better upper bounds than the general purpose heuristic solvers (LocalSolver3.1 and the Wedelin’s heuristic) and

18

Table 11: Computational results of the latest solvers and the proposed algorithm for reduced
SCP instances
instance
?G.1–5 (5)
?H.1–5 (5)
?I.1–5 (5)
?J.1–5 (5)
?K.1–5 (5)
?L.1–5 (5)
?M.1–5 (5)
?N.1–5 (5)
RAIL507
RAIL516
RAIL582
RAIL2536
?RAIL2586
?RAIL4284
?RAIL4872
avg. (all)
avg. (with stars)

CPLEX12.6
0.60%
1.92%
2.59%
4.43%
2.84%
4.77%
10.82%
15.95%
0.00%
0.00%
0.00%
0.00%
0.84%
0.47%
0.84%
4.72%
5.16%

Gurobi5.6.3
0.24%
1.62%
1.64%
3.99%
2.48%
4.85%
4.56%
11.36%
0.00%
0.00%
0.00%
0.00%
0.52%
0.93%
1.23%
3.33%
3.64%

SCIP3.1
0.47%
1.62%
2.10%
3.42%
2.66%
3.13%
31.15%
34.17%
0.00%
0.00%
0.00%
0.00%
0.73%
0.74%
1.10%
8.43%
9.21%

LocalSolver3.1
3.78%
1.90%
1.74%
2.99%
2.18%
2.41%
4.89%
6.47%
0.57%
1.61%
0.47%
0.29%
1.66%
1.66%
1.73%
2.98%
3.18%

Yagiura et al.
0.00%
0.00%
0.00%
0.15%
0.00%
0.00%
0.00%
0.00%
0.00%
0.00%
0.00%
0.00%
0.00%
0.00%
0.00%
0.02%
0.02%

proposed
0.00%
0.00%
0.00%
0.31%
0.63%
0.77%
1.01%
1.06%
0.00%
0.00%
0.00%
0.00%
0.84%
0.84%
0.84%
0.45%
0.50%

Table 12: Computational results of the latest solvers and the proposed algorithm for SPP
instances
instance
CPLEX12.6
Gurobi5.6.3
aa01–06 (6)
0.00%(6)
0.00%(6)
us01–04 (4)
0.00%(4)
0.00%(4)
v0415–0421 (7)
0.00%(7)
0.00%(7)
v1616–1622 (7)
0.00%(7)
0.00%(7)
t0415–0421 (7)
0.66%(7)
0.60%(7)
?t1716–1722 (7)
8.34%(7)
16.58%(7)
?ds
8.86%
55.61%
?ds-big
62.16%
24.03%
?ivu06-big
20.86%
0.68%
?ivu59
28.50%
4.36%
avg. (all)
4.37%(42/42)
4.88%(42/42)
avg. (with stars)
14.10%(13/13)
15.66%(13/13)
† 600 s on a 1.3 GHz Sun Ultra Sparc-IIIi
‡ 7200 s on a Sun Ultra Sparc 1 Model 170E

SCIP3.1
0.00%(6)
0.00%(3)
0.00%(7)
0.00%(7)
1.61%(6)
3.51%(7)
40.53%
72.01%
17.90%
37.84%
5.06%(40/42)
15.07%(13/13)

LocalSolver3.1
13.89%(1)
11.26%(2)
0.05%(7)
4.60%(7)
— (0)
37.08%(1)
85.17%
92.69%
52.54%
48.95%
17.52%(22/42)
63.29%(5/13)

Bastert et al.†
—
—
0.71%(6)
6.64%(3)
1.30%(5)
12.55%(7)
8.82%
—
—
—
5.79%(22/29)
10.79%(9/10)

Borndörfer‡
0.00%(6)
0.00%(4)
0.00%(7)
0.01%(7)
1.83%(7)
1.63%(7)
—
—
—
—
0.64%(38/38)
1.77%(9/9)

good upper bounds comparable to those of the branch-and-bound algorithm
(Borndörfer, 1998) for many SPP instances. In view of these, our algorithm
achieves sufficiently good upper bounds compared to the other algorithms for
the SPP instances, particularly for hard SPP instances.
8. Conclusion
We present a data mining approach for reducing the search space of local
search algorithms for a class of BIPs including SCP and SPP. In this approach,
we construct a k-nearest neighbor graph by extracting variable associations from
the instance to be solved in order to identify promising pairs of flipping variables
in the 2-flip neighborhood. We also develop a 4-flip neighborhood local search
algorithm that flips four variables alternately along 4-paths or 4-cycles in the
k-nearest neighbor graph. We incorporate an efficient incremental evaluation of
solutions and an adaptive control of penalty weights into the 4-flip neighborhood
local search algorithm. Computational results show that the proposed method
improves the performance of the local search algorithm for large-scale SCP and
SPP instances.
We expect that data mining approaches could also be beneficial for efficiently
solving other large-scale combinatorial optimization problems, particularly for

19

proposed
1.60%(6)
0.04%(4)
0.00%(7)
0.09%(7)
0.92%(6)
2.45%(7)
0.00%
0.00%
0.00%
0.00%
0.81%(41/42)
1.43%(12/13)

hard instances having large gaps between the lower and upper bounds of the
optimal values.
Acknowledgment
This work was supported by the Grants-in-Aid for Scientific Research (JP26282085).
References
Achterberg, T. (2009). SCIP: Solving constraint integer programs. Mathematical
Programming Computation, 1, 1–41.
Agarwal, Y., Mathur, K., & Salkin, H. M. (1989). A set-partitioning-based exact
algorithm for the vehicle routing problem. Networks, 19, 731–749.
Atamtürk, A., Nemhauser, G. L., & Savelsbergh, M. W. P. (1995). A combined
Lagrangian, linear programming, and implication heuristic for large-scale set
partitioning problems. Journal of Heuristics, 1, 247–259.
Balas, E., & Padberg, M. W. (1976). Set partitioning: A survey. SIAM Review,
18, 710–760.
Baldacci, R., Christofides, N., & Mingozzi, A. (2008). An exact algorithm for
the vehicle routing problem based on the set partitioning formulation with
additional cuts. Mathematical Programming, 115, 351–385.
Barahona, F. & Anbil, R. (2000). The volume algorithm: Producing primal
solutions with a subgradient method. Mathematical Programming, 87, 385–
399.
Barnhart, C., Johnson, E. L., Nemhauser, G. L., Savelsbergh, M. W. P., &
Vance, P. H. (1998). Branch-and-price: Column generation for solving huge
integer programs. Operations Research, 46, 316–329.
Bastert, O., Hummel, B. & de Vries, S. (2010). A generalized Wedelin heuristic
for integer programming. INFORMS Journal on Computing, 22, 93–107.
Beasley, J. E. (1990). OR-Library: Distributing test problems by electronic mail.
Journal of the Operational Research Society, 41, 1069–1072.
Benoist, T., Estellon, B., Gardi, F., Megel, R., & Nouioua, K. (2011). LocalSolver 1.x: A black-box local-search solver for 0-1 programming. 4OR — A
Quarterly Journal of Operations Research, 9, 299–316.
Borndörfer, R. (1998). Aspects of set packing, partitioning and covering. Ph. D.
Dissertation, Berlin: Technischen Universität.
Boros, E., Hammer, P. L., Ibaraki, T., Kogan, A., Mayoraz, E., & Muchnik, I.
(2000). An implementation of logical analysis of data. IEEE Transactions on
Knowledge and Data Engineering, 12, 292–306.
20

Boros, E., Ibaraki, T., Ichikawa, H., Nonobe, K., Uno, T., & Yagiura, M.
(2005). Heuristic approaches to the capacitated square covering problem. Pacific Journal of Optimization, 1, 465–490.
Boschetti, M. A., Mingozzi, A. & Ricciardelli, S. (2008). A dual ascent procedure
for the set partitioning problem. Discrete Optimization, 5, 735–747.
Bramel, J., & Simchi-Levi, D. (1997). On the effectiveness of set covering formulations for the vehicle routing problem with time windows. Operations
Research, 45, 295–301.
Caprara, A., Fischetti, M., & Toth, P. (1999). A heuristic method or the set
covering problem. Operations Research, 47, 730–743.
Caprara, A., Toth, P., & Fischetti, M. (2000). Algorithms for the set covering
problem. Annals of Operations Research, 98, 353–371.
Caserta, M. (2007). Tabu search-based metaheuristic algorithm for large-scale
set covering problems. In W. J. Gutjahr, R. F. Hartl, & M. Reimann (eds.),
Metaheuristics: Progress in Complex Systems Optimization (pp. 43–63).
Berlin: Springer.
Ceria, S., Nobili, P., Sassano, A. (1997). Set covering problem. In M. Dell’Amico,
F. Maffioli & S. Martello (eds.): Annotated Bibliographies in Combinatorial
Optimization, (pp. 415–428). New Jersey: John Wiley & Sons.
Ceria, S., Nobili, P., & Sassano, A. (1998). A Lagrangian-based heuristic for
large-scale set covering problems. Mathematical Programming, 81, 215–288.
Chu P. C., & Beasley, J. E. (1998). Constraint handling in genetic algorithms:
The set partitioning problem. Journal of Heuristics, 11, 323–357.
CPLEX Optimizer. (2014). http://www-01.ibm.com/software/commerce/
optimization/cplex-optimizer/index.html Accessed 2017.01.23.
Farahani, R. Z., Asgari, N., Heidari, N., Hosseininia, M., & Goh, M. (2012).
Covering problems in facility location: A review. Computers & Industrial
Engineering, 62, 368–407.
Gurobi Optimizer. (2014). http://www.gurobi.com/ Accessed 2017.01.23.
Hammer, P. L., & Bonates, T. O. (2006). Logical analysis of data — An
overview: From combinatorial optimization to medical applications. Annals
of Operations Research, 148, 203–225.
Hashimoto, H., Ezaki, Y., Yagiura, M., Nonobe, K., Ibaraki, T., & Løkketangen,
A. (2009). A set covering approach for the pickup and delivery problem with
general constraints on each route. Pacific Journal of Optimization, 5, 183–200.
Hoffman, K. L., & Padberg, A. (1993). Solving airline crew scheduling problems
by branch-and-cut. Management Science, 39, 657–682.
21

Hutter, F., Tompkins, D. A. D., & Hoos, H. H. (2002). Scaling and probabilistic smoothing: Efficient dynamic local search for SAT. Proceedings of International Conference on Principles and Practice of Constraint Programming
(CP), 233–248.
Johnson, D. S., & McGeoch, L. A. (1997). The traveling salesman problem: A
case study. In E. Aarts, & K. Lenstra (eds.), Local Search in Combinatorial
Optimization (pp. 215–310). New Jersey: Princeton University Press.
Koch, T., Achterberg, T., Andersen, E., Bastert, O., Berthold, T., Bixby, R. E.,
Danna, E., Gamrath, G., Gleixner, A. M., Heinz, S, Lodi, A., Mittelmann, H.,
Ralphs, T., Salvagnin, D., Steffy, D. E., & Wolter, K. (2011). MIPLIB2010:
Mixed integer programming library version 5. Mathematical Programming
Computation, 3, 103–163.
Linderoth, J. T., Lee, E. K., & Savelbergh, M. W. P. (2001). A parallel, linear
programming-based heuristic for large-scale set partitioning problems. INFORMS Journal on Computing, 13, 191–209.
Lodi, A. (2013) The heuristic (dark) side of MIP solvers. In E. -G. Talbi (ed.),
Hybrid Metaheuristics (pp. 273–284). Berlin: Springer.
Michel, L. & Van Hentenryck, P. (2000). Localizer. Constraints: An International Journal, 5, 43–84.
Mingozzi, A., Boschetti, M. A., Ricciardelli, S., & Bianco, L. (1999). A set
partitioning approach to the crew scheduling problem. Operations Research,
47, 873–888.
Morris, P. (1993). The breakout method for escaping from local minima. Proceedings of National Conference on Artificial Intelligence (AAAI), 40–45.
Nonobe, K., & Ibaraki, T. (2001). An improved tabu search method for the
weighted constraint satisfaction problem. INFOR, 39, 131–151.
Pesant, G., & Gendreau, M. (1999). A constraint programming framework for
local search methods. Journal of Heuristics, 5, 255–279.
Selman, B., & Kautz, H. (1993). Domain-independent extensions to GSAT:
Solving large structured satisfiability problems. Proceedings of International
Conference on Artificial Intelligence (IJCAI), 290–295.
Shaw, P., Backer, B. D., & Furnon, V. (2002). Improved local search for CP
toolkits. Annals of Operations Research, 115, 31–50.
Thornton, J. (2005). Clause weighting local search for SAT. Journal of Automated Reasoning, 35, 97–142.
Umetani, S., & Yagiura, M. (2007). Relaxation heuristics for the set covering
problem. Journal of the Operations Research Society of Japan, 50, 350–375.

22

Umetani, S., Arakawa, M., & Yagiura, M. (2013). A heuristic algorithm for the
set multicover problem with generalized upper bound constraints. Proceedings
of Learning and Intelligent Optimization Conference (LION), 75–80.
Umetani, S. (2015). Exploiting variable associations to configure efficient local
search in large-scale set partitioning problems. Proceedings of AAAI Conference on Artificial Intelligence (AAAI), 1226–1232.
Van Hentenryck, P., & Michel, L. (2005). Constraint-Based Local Search, Cambridge: The MIT Press.
Voudouris, C., Dorne, R., Lesaint, D., & Liret, A. (2001). iOpt: A software
toolkit for heuristic search methods. Proceedings of Principles and Practice
of Constraint Programming (CP), 716–729.
Wedelin, D. (1995). An algorithm for large-scale 0-1 integer programming with
application to airline crew scheduling. Annals of Operations Research, 57,
283–301.
Yagiura, M. & Ibaraki, T. (1999). Analysis on the 2 and 3-flip neighborhoods
for the MAX SAT. Journal of Combinatorial Optimization, 3, 95–114.
Yagiura, M., Kishida, M., & Ibaraki, T. (2006). A 3-flip neighborhood local
search for the set covering problem. European Journal of Operational Research, 172, 472–499.

23

