Evaluating and Tuning n-fold Integer
Programming
Kateřina Altmanová
Department of Applied Mathematics, Charles University
Prague, Czech Republic
kacka@kam.mff.cuni.cz

Dušan Knop1

arXiv:1802.09007v1 [cs.DS] 25 Feb 2018

Department of Informatics, University of Bergen
Bergen, Norway
dusan.knop@uib.no
0000-0003-2588-5709

Martin Koutecký2
Faculty of Industrial Engineering and Management, Technion – Israel Institute of Technology
Haifa, Israel
koutecky@technion.ac.il
0000-0002-7846-0053

Abstract
In recent years, algorithmic breakthroughs in stringology, computational social choice, scheduling,
etc., were achieved by applying the theory of so-called n-fold integer programming. An n-fold
integer program (IP) has a highly uniform block structured constraint matrix. Hemmecke, Onn,
2
and Romanchuk [Math. Programming, 2013] showed an algorithm with runtime aO(rst+r s) n3 ,
where a is the largest coefficient, r, s, and t are dimensions of blocks of the constraint matrix
and n is the total dimension of the IP; thus, an algorithm efficient if the blocks are of small size
and with small coefficients. The algorithm works by iteratively improving a feasible solution with
augmenting steps, and n-fold IPs have the special property that augmenting steps are guaranteed
to exist in a not-too-large neighborhood. However, this algorithm has never been implemented
and evaluated.
We have implemented the algorithm and learned the following along the way. The original
algorithm is practically unusable, but we discover a series of improvements which make its evaluation possible. Crucially, we observe that a certain constant in the algorithm can be treated
as a tuning parameter, which yields an efficient heuristic (essentially searching in a smaller-thanguaranteed neighborhood). Furthermore, the algorithm uses an overly expensive strategy to find
a “best” step, while finding only an “approximatelly best” step is much cheaper, yet sufficient
for quick convergence. Using this insight, we improve the asymptotic dependence on n from n3
to n2 log n which yields the currently asymptotically fastest algorithm for n-fold IP.
Finally, we tested the behavior of the algorithm with various values of the tuning parameter
and different strategies of finding improving steps. First, we show that decreasing the tuning
parameter initially leads to an increased number of iterations needed for convergence and eventually to getting stuck in local optima, as expected. However, surprisingly small values of the
parameter already exhibit good behavior. Second, our new strategy for finding “approximatelly
best” steps wildly outperforms the original construction.
2012 ACM Subject Classification Software and its engineering → Software design engineering,
Theory of computation → Parameterized complexity and exact algorithms

1
2

Author supported supported by the project NFR MULTIVAL and the project P202/12/G061 of GA ČR.
Author supported by Technion postdoctoral fellowship.

23:2

Evaluating and Tuning n-fold Integer Programming

Keywords and phrases n-fold integer programming, integer programming, analysis of algorithms, primal heuristic, local search
Supplement Material https://github.com/katealtmanova/nfoldexperiment

1

Introduction

In this article we consider the general integer linear programming (ILP) problem in standard
form,
min {wx | Ax = b , l ≤ x ≤ u , x ∈ Zn } .

(ILP)

with A an integer m × n matrix, b ∈ Z , w ∈ Z , l, u ∈ (Z ∪ {±∞}) . It is well known to
be strongly NP-hard, but models many important problems in combinatorial optimization
such as planning [27], scheduling [12] and transportation [3] and thus powerful generic solvers
have been developed for it [24]. Still, theory is motivated to search for tractable special cases
One such special case is when the constraint matrix A has a so-called N -fold structure:


E1 E1 · · · E1
 E2 0 · · · 0 




A = E (N ) =  0 E2 · · · 0  .
 .
..
.. 
..
 ..
.
.
. 
m

0

0

···

n

n

E2

Here, r, s, t, N ∈ N, u, l, w ∈ ZN t , b ∈ Zr+N s , E (N ) is an (r + N s) × N t-matrix, E1 ∈ Zr×t
is an r ×t-matrix and E2 ∈ Zs×t is an s × t-matrix. We call E (N ) the N -fold product of
(N )
1
E= E
is known as N -fold integer programming (N -fold
E2 . Problem (ILP) with A = E
IP). Hemmecke, Onn, and Romanchuk [15] prove the following.
I Proposition 1 ([15, Theorem 6.2]). There is an algorithm that solves3 (ILP) with A = E (N )
2
encoded with L bits in time ∆O(trs+t s) · n3 L, where ∆ = 1 + max{kE1 k∞ , kE2 k∞ }.
Recently, algorithmic breakthroughs in stringology [19], computational social choice [20],
scheduling [4, 17, 21], etc., were achieved by applying this algorithm and its subsequent
non-trivial improvements.
The algorithm belongs to the larger family of augmentation (primal) algorithms. It
starts with an initial feasible solution x0 ∈ ZN t and produces a sequence of increasingly
better solutions x1 , . . . , xs (better means wxs < wxs−1 < · · · < wx0 ). It is guaranteed that
the algorithm terminates, that xs is an optimal solution, and that the algorithm converges
quickly, i.e., s is polynomial in the length of the input. A key property of N -fold IPs is
that, if an augmenting step exists, then it can be decomposed into a bounded number of
elements of the so-called Graver basis of A, which in turn makes it possible to compute it
using dynamic programming [15, Lemma 3.1]. In a sense, this property makes the algorithm
a local search algorithm which is always guaranteed to find an improvement in a not-too-large
neighborhood. The bound on the number of elements or the size of the neighborhood which
needs to be searched is called the Graver complexity of E. This, in turn, implies that, if an
augmenting step exists, then there is always one with small `1 -norm; for a matrix A, we
denote this bound g1 (A) = maxg∈G(A) kgk1 [23, Theorem 4]. However, the algorithm has
never been implemented and evaluated.
3

Given an IP, we say that to solve it is to either (i) declare it infeasible or unbounded or (ii) find a
minimizer of it.

K. Altmanová, D. Knop, and M. Koutecký

1.1

23:3

Our Contributions

We have implemented the algorithm and tested it on two problems for which n-fold formulations were known: makespan minimization on uniformly related machines (Q||Cmax ) and
Closest String; we have used randomly generated instances.
In the course of implementing the algorithm we learn the following. The algorithm in its
initial form is practically unusable due to an a priori construction of the Graver basis G(E2 )
of size exponential in r, s, t and a related (even larger) set Z(E). However, we discover a
series of improvements (some building on recent insights [23]) which avoid the construction
of these two sets. Moreover, we adjust the algorithm to treat g1 (A) as a tuning parameter
g1 , which turns it into a heuristic.
We also study the augmentation strategy, which is the way the algorithm chooses an
augmenting step among all the possible options. The original algorithm uses an overly
expensive strategy to find a “best” step, which means that a large number of possible steps
is evaluated in each iteration. We show that finding only an “approximatelly best” step is
sufficient to obtain asymptotically equivalent convergence rate, and the work per iteration
decreases exponentially. Using this insight, we improve the asymptotic dependence on N from
N 3 to N 2 log N . Together with recent improvements, this yields the currently asymptotically
fastest algorithm for N -fold IP:
I Theorem 2. Problem (ILP) with A = E (N ) can be solved in time ∆r
where M = log(wx∗ − wx0 ) for some minimizer x∗ of wx.

2

s+rs2

(N t)2 log(N t)M ,

Finally, we evaluate the behavior of the algorithm. We ask how is the performance of the
algorithm (in terms of number of dynamic programming calls and quality of the returned
solution) influenced by
1. the choice of the tuning parameter 1 < g1 ≤ g1 (A)?
2. the choice of the augmentation strategy between “best step”, “approximate best step”
and “any step”?
As expected, as g1 moves from g1 (A) to 1, we first see an increase in the number of iterations
needed for convergence and eventually the algorithm gets stuck in local optima. However,
surprisingly small values (e.g. g1 = 5 when g1 (A) > 105 ) of the parameter already exhibit
close to optimal behavior. Second, our new strategy for finding “approximatelly best” steps
outperforms the original construction by orders of magnitude, while the naive “any step”
strategy behaves erratically.
We note that at this stage we are not (yet) interested in showing supremacy over existing
algorithms; we simply want to understand the practical behavior of an algorithm whose
theoretical importance was recently highlighted. For this reasons our experimental focus is
on the two aforementioned questions rather than simply measuring the time. Similarly, due
to the rigid format of E (N ) we are limited to few problems for which N -fold formulations are
known. For Closest String we use the same instances as Chimani et al. [5]; for Makespan
Minimization we generate our own data because standard benchmarks are not limited to
short jobs or few types of jobs.

1.2

Related Work

Our work mainly relates to primal heuristics [2] for MIPs which are used to help reach optimality faster and provide good feasible solutions early in the termination process. Specifically,
our algorithm is a neighborhood (or local) search algorithm. The standard paradigm is Large
Neighborhood Search (LNS) [26] with specializations such as for example Relaxation Induced

23:4

Evaluating and Tuning n-fold Integer Programming

Neighborhood Search (RINS) [6] and Feasibility Pump [1]. Using this paradigm, our proposed
algorithm searches in the neighborhood induced by the `1 -distance from the current feasible
solution and the search procedure is formulated as an ILP subproblem with the additional
constraint kxk1 ≤ g1 . In this sense the closest technique to ours is local branching [10] which
also searches in the `1 -neighborhood; however, it was mainly applied to binary programs
without any additional structure of the constraint matrix.
On the theoretical side, very recently Koutecký et al. [23] have studied parameterized
strongly polynomial algorithms for various block-structured ILPs, not just N -fold IP. Eisenbrand et al. [8] independently (and using slightly different techniques) arrive at the same
complexity of N -fold IP as our Theorem 2.

2

Preliminaries

For positive integers m, n we set [m, n] = {m, . . . , n} and [n] = [1, n]. We write vectors in
boldface (e.g., x, y) and their entries in normal font (e.g., the i-th entry of x is xi ). Given
the problem (ILP), we say that x is feasible for (ILP) if Ax = b and l ≤ x ≤ u.
Graver bases and augmentation. Let us now introduce Graver bases and discuss how they
can be used for optimization. We also recall N -fold IPs; for background, we refer to the
books of Onn [25] and De Loera et al. [7].
Let x, y be n-dimensional integer vectors. We call x, y sign-compatible if they lie in
the same orthant, that is, for each i ∈ [n] the sign of xi and yi is the same. We call
P i
i
i g a sign-compatible sum if all g are pair-wise sign-compatible. Moreover, we write
y v x if x and y are sign-compatible and |yi | ≤ |xi | for each i ∈ [n]. Clearly, v imposes
a partial order called “conformal order” on n-dimensional vectors. For an integer matrix
A ∈ Zm×n , its Graver basis G(A) is the set of v-minimal non-zero elements of the lattice of
A, kerZ (A) = {z ∈ Zn | Az = 0}. An important property of G(A) is the following.
I Proposition 3 ([25, Lemma 3.4]). Every integer vector x 6= 0 with Ax = 0 is a sign-comPt
patible sum x = i=1 αi gi , αi ∈ N, gi ∈ G(A) and t ≤ 2n − 2.
Let x be a feasible solution to (ILP). We call g an x-feasible step (or simply feasible
if x is clear) if x + g is feasible for (ILP). Further, we call a feasible step g augmenting
if w(x + g) < wx; note that g decreases the objective by wg. An augmenting step g
and a step length γ ∈ N form an x-feasible step pair with respect to a feasible solution x
if l ≤ x + γg ≤ u. A pair (γ, g) ∈ (N × G(A)) is a γ-Graver-best step pair and γg is a
γ-Graver-best step if it is feasible and for every feasible step pair (γ, g0 ), g0 ∈ G(A), we have
wγg ≤ wγg0 . An augmenting step g and a step length γ ∈ N form a Graver-best step pair if
it is γ-Graver-best and it minimizes wγ 0 g0 over all γ 0 ∈ N, where (γ 0 , g0 ) is a γ 0 -Graver-best
step pair. We say that γg is a Graver-best step if (γ, g) is a Graver-best step pair.
The Graver-best augmentation procedure for (ILP) with a given feasible solution x0 and
initial value i = 0 works as follows:
1. If there is no Graver-best step for xi , return it as optimal.
2. If a Graver-best step γg for xi exists, set xi+1 := xi + γg, i := i + 1, and go to 1.
I Proposition 4 (Convergence bound [25, Lemma 3.10]). Given a feasible solution x0 for (ILP),
the Graver-best augmentation procedure finds an optimum in at most (2n − 2) log M steps,
where M = w(x0 − x∗ ) and x∗ is any minimizer of wx.
By standard techniques (detecting unboundedness etc.) we can ensure that log M ≤ L.

K. Altmanová, D. Knop, and M. Koutecký

23:5

N -fold IP. The structure of E (N ) allows us to divide the N t variables of x into N bricks of
size t. We use subscripts to index within a brick and superscripts to denote the index of the
brick, i.e., xij is the j-th variable of the i-th brick with j ∈ [t] and i ∈ [N ].

3

Approximate Graver-best Steps

In this section we introduce the notion of a c-approximate Graver-best step (Definition 5),
show that such steps exhibit good convergence (Lemma 6), can be easily obtained (Lemma 7),
and result in a significant speed-up of the N -fold IP algorithm (Theorem 2).
I Definition 5 (c-approximate Graver-best step). Let c ∈ R, c ≥ 1. Given an instance of (ILP)
and a feasible solution x, we say that h is a c-approximate Graver-best step for x if, for every
x-feasible step pair (γ, g) ∈ (N × G(A)), we have wh ≤ 1c · γwg.
Recall the Graver-best augmentation procedure. We call its analogue where we replace
a Graver-best step with a c-approximate Graver-best step the c-approximate Graver-best
augmentation procedure.
I Lemma 6 (c-approximate convergence bound). Given a feasible solution x0 for (ILP), the
c-approximate Graver-best augmentation procedure finds an optimum of (ILP) in at most
c · (2n − 2) log M steps, where M = w(x0 − x∗ ) and x∗ is any minimizer of wx.
Proof. The proof is a straightforward adaptation of the proof of Proposition 4. Let x∗ be
P2n−2
a minimizer and let h = x∗ − x0 . Since Ah = 0, by Proposition 3, h = i=1 αi gi for
some αi ∈ N, gi ∈ G(A), i ∈ [2n − 2]. Thus, an x-feasible step pair (γ, g) such that γg is
1
M . In other words, any Graver-best step pair
a Graver-best step must satisfy wγg ≤ 2n−2
1
improves the objective function by at least a 2n−2
-fraction of the total optimality gap M ,
and thus (2n − 2) log M steps suffice to reach an optimum. It is straightforward to see that a
c
c-approximate Graver-best step satisfies wx − w(x + γg) ≤ 2n−2
M , and thus c(2n − 2) log M
steps suffice.
J
I Remark. Lemma 6 extends naturally to separable objectives; see the original proof [25,
Lemma 3.10].
I Lemma 7 (Powers of c step lengths). Let c ∈ N, x be a feasible solution of (ILP), and let

Γc-apx = ci | ∃g ∈ G(A) : l ≤ x + ci g ≤ u .
Let (γ, g) ∈ (Γc-apx × G(A)) be an x-feasible step pair such that γg ≤ γ 0 g0 for any x-feasible
step pair (γ 0 , g0 ) ∈ (Γc-apx × G(A)). Then γg is a c-approximate Graver-best step.
Proof. Let (γ, g) satisfy the assumptions, and let (γ̃, g̃) ∈ (N × G(A)) be a Graver-best step
pair. Let γ 0 be a nearest smaller power of c from γ̃, and observe that γ 0 g̃ is a c-approximate
Graver-best step because γ 0 ≥ γ̃c . On the other hand, since γg is a γ-Graver-best step, we
have γg ≤ γ 0 g̃ and thus γg is also a c-approximate Graver-best step.
J
I Theorem 2. Problem (ILP) with A = E (N ) can be solved in time ∆O(r
where M = log(wx∗ − wx0 ) for some minimizer x∗ of wx.

2

s+rs2 )

(N t)2 log(N t)M ,

Proof. Recall that ∆ = kAk∞ +1. Koutecký et al. [23, Theorem 2] show that a γ-Graver-best
2
2
step can be found in time ∆r s+rs N t. Moreover, Hemmecke et al. [14] prove a proximity

23:6

Evaluating and Tuning n-fold Integer Programming

theorem which allows the reduction of an instance of (ILP) to an equivalent instance with
new bounds l0 , u0 satisfying ku0 − l0 k∞ ≤ N tg∞ , with
g∞ = max kgk∞ ≤ max kgk1 ≤ (∆rs)O(rs) ,
g∈G(A)

g∈G(A)

where the last inequality can be found in the proof of [23, Theorem 4]. This bound implies
that Γ2-apx from in Lemma 7 satisfies |Γ2-apx | ≤ rs log(∆N trs). By Lemma 7, finding
a γ-Graver-best for each γ ∈ Γ2-apx and picking the minimum results in a 2-approximate
2
2
Graver-best step, and can be done in time ∆r s+rs (N t) log(N t). By Lemma 6, (4n−4) log M
steps suffice to reach the optimum.
J

4

Implementation

We first give an overview of the original algorithm, which is our starting point. Then we
discuss our specific improvements and mention a few details of the software implementation.

4.1

Overview of the Original Algorithm

The key property of the N -fold product E (N ) is that, for any N ∈ N, the number of nonzero
bricks of any g ∈ G(E (N ) ) is bounded by some constant g(E) called the Graver complexity
of E, and, moreover, that the sum of all non-zero bricks of g can be decomposed into at most
g(E) elements of G(E2 ) [15, Lemma 3.1]. This facilitates the following construction. Let
(
Z(E) =

z ∈ Zt | ∃g1 , . . . , gk ∈ G(E2 ), k ≤ g(E), z =

k
X

)
gi

.

i=1

Then, every prefix sum of the bricks of g ∈ G(E (N ) ) is contained in Z(E) and a γ-Graver-best
step, γ ∈ N, can be found using dynamic programming over the elements of Z(E).
To ensure that a Graver-best step is found, a set of step-lengths Γbest is constructed as

follows. Observe that any Graver-best (and thus feasible) step pair (γ, g) ∈ N × G(E (N ) ) ,
must satisfy that in at least one brick i ∈ [n] it is “tight”, that is, (γ, g) is x-feasible while
(γ + 1, g) is not specifically because li ≤ x + γg ≤ ui holds but li ≤ x + (γ + 1)g ≤ ui does
not. Thus, for each z ∈ Z(E) and each i ∈ [n], we find all the potentially “tight” step lengths
γ and add them to Γbest , which results in a bound of |Γbest | ≤ |Z(E)| · n.

4.2

Replacing Dynamic Programming with ILP

We have started off by implementing the algorithm exactly as it is described by Hemmecke et
al. [15]. The first obstacle is encountered almost immediately and is contained in the constant
g(E). This constant can be computed, but the computation is extremely difficult [9, 13].
Another possibility is to estimate it, in which case it is almost always larger than N and
thus is essentially meaningless. Finally, one can take the approach partially suggested in [15,
Section 7], where we consider g(E) in the construction of Z(E) to be a tuning parameter
and consider the approximate set Zgc (E), gc ∈ N, obtained by taking sums of at most gc
elements of G(E2 ). This makes the algorithm more practical, but turns it into a heuristic.
In spite of this sacrifice, already for small (r = 3, s = 1, t = 7, n = 10) instances and
extremely small value of gc = 3, the dynamic programming based on the Zgc (E) construction
was taking an unreasonably long time (over one minute). Admittedly this could be improved;
however, already for gc > 5, it becomes infeasible to compute Zgc (E), and for larger instances

K. Altmanová, D. Knop, and M. Koutecký

1
2
3
4
5
6
7
8
9

input : matrices E1 , E2 , positive integer N , and vectors b, l, u, w
output : optimal solution to (ILP) with A = E (N )
g = GraverComplexity(E1 , E2 );
x0 = FindFeasibleSolution(E, b, l, u, w), i = 0;
G(E1 ) = GraverBasis(E1 , g);
Z(E) = DynamicProgramStates(G(E1 ), g);
do
Γ = BuildGamma(xi );
i = i + 1;
foreach γ ∈ Γ do
gγ = gammaBestStep(γ, g);
xi = xi−1 + argmin{g

10
11
12

23:7

} wgγ ;

γ |γ∈Γ

while xi−1 6= xi ;
return xi ;
Algorithm 1: Pseudocode of the algorithm of Hemmecke, Onn, and Romanchuk.

(r > 5, t > 12) it becomes very difficult to compute even G(E2 ). For these reasons we sought
to completely replace the dynamic program involving Z(E).
Koutecký et al. [23] show that all instances of (ILP) with the property that the so-called
dual treedepth tdD (A) of A is bounded and the largest coefficient kAk∞ is bounded also have
the property that g1 (A) = maxg∈G(A) kgk1 is bounded, which implies that augmenting steps
can be found efficiently. This class of ILPs contains N -fold IP.
The interpretation of the above fact is that, in order to solve (ILP), it is sufficient to
repeatedly (for different x and γ) solve an auxiliary instance
min {wh | Ah = 0, l ≤ x + γh ≤ u, khk1 ≤ g1 (A)}

(AugILP)

in order to find good augmenting steps; we note that the constraint khk1 ≤ g1 (A) can be
linearized [23, Lemma 25]. The heuristic approach outlined above transfers easily: we replace
g1 (A) in (AugILP) with some integer g1 , 1 < g1 ≤ g1 (A); this makes (AugILP) easier to
solve at the cost of losing the guarantee that an augmenting step is found if one exists. In
theory, solving (AugILP) should be easier than solving the original instance (ILP) due to the
special structure of A [23, Lemma 25]. Our approach here is to simply invoke an industrial
MILP solver on (AugILP) in order to find a γ-Graver-best step.

4.3

Augmentation Strategy: Step Lengths

Logarithmic Γ. The majority of algorithms based on Graver basis augmentation rely on
the Graver-best augmentation procedure [4, 7, 15, 19, 21, 25] and thus require finding
(exact) Graver-best steps. In the aforementioned algorithms this is always done using the
construction of a set Γbest mentioned above, which is of size f (k) · n where k is the relevant
2
parameter (e.g., (ars)O(rst+st ) in the original algorithm for N -fold IP). We replace this
construction with Γ2-apx = {1, 2, 4, 8, . . . } which, combined with the proximity technique, is
only of size O(log n) (Theorem 2); in particular, independent of the function f (k).
Exhausting γ. Moreover, we have noticed that sometimes the algorithm finds a step g for
γ = 2k which is not tight in any brick, and then repeatedly applies it for shorter step-lengths
γ 0 < γ. In other words, the discovered direction g is not exhausted. Thus, for each γ ∈ N,

23:8

Evaluating and Tuning n-fold Integer Programming

1
2
3
4
5
6
7
8
9
10

input : matrices E1 , E2 , positive integers N , c, and g1 , and vectors b, l, u, w
output : a feasible solution to (ILP) with A = E (N )
x0 = FindFeasibleSolution(E, b, l, u, w), i = 0;
do
Γ = ∅; j = 0, i = i + 1;
do
γ = cj ;
gγ = min{wγg | Ag = 0, l ≤ x + γg ≤ u, kgk1 ≤ g1 , g ∈ ZN t };
γ 0 = ExhaustDirection(gγ );
Γ = Γ ∪ {γ 0 }, j = j + 1;
while gγ 6= 0;
xi = xi−1 + argmin{g |γ∈Γ} wgγ ;
γ

while xi−1 6= xi ;
12 return xi
Algorithm 2: Pseudocode of our new heuristic algorithm. The algorithm is exact if
g1 ≥ g1 (A) = maxg∈G(A) kgk1 .
11

upon finding the γ-Graver-best step g, we replace γ with the largest γ 0 ≥ γ for which (γ 0 , g)
is still x-feasible.
Early termination. Another observation is that in any given iteration of the algorithm, if
γ > 1 then some augmenting step has been found and if the computation is taking too long,
we might terminate it and simply apply the best step found so far.

4.4

Software and Hardware

We have implemented our solver in the SageMath computer algebra system [30]. This was a
convenient choice for several reasons. The SageMath system offers an interactive notebookstyle web-based interface, which allows rapid prototyping and debugging. Data types for
vectors and matrices and a unified interface for several MILP solvers are also readily available.
We have experimented with the open-source solvers GLPK [29], Coin-OR CBC [28], and
the commercial solver Gurobi [11] and have settled for using the latter since it performs the
best. The downside of SageMath is that an implementation of the original dynamic program
is likely much slower than a similar implementation in C; however this DP is impractical
anyway as explained in Section 4.2. For random instance generation and subsequent data
evaluation and graphing, we have used the Jupyter notebook environment [18] and Matplotlib
library [16]. The computations were performed on a computer with an Intel® Xeon® E5-2630
v3 (2.40GHz) CPU and 128 GB RAM.

5

Evaluation

We begin our evaluation with two main questions, specifically, how is the performance of
the algorithm (both in terms of the number of iterations and the quality of the returned
solution) influenced by:
1. the tuning parameter g1 and
2. the augmentation strategy?

K. Altmanová, D. Knop, and M. Koutecký

Regarding our first question, theoretically we should see either an increase in the number
of iterations, a decrease in the quality of the returned solution, or both. However, the range
of the tuning parameter g1 is quite large: any number between 2 and g1 (A) is a valid choice,
and in all our scenarios the true value of g1 (A) exceeds 300. Thus, we are interested in the
threshold values of g1 when the algorithm no longer finds the true optimum or when its
convergence rate drops significantly.
Regarding our second question, there are two main candidates for the set of step-lengths
Γ. We can either use the “best step” construction Γbest of the original algorithm, which
assures that we always make a Graver-best step before moving to the next iteration. Or,
we can use the “approximate best step” construction Γ2-apx of Theorem 2, which provides a
2-approximate Graver-best step. To make this comparison more interesting, we also consider
Γ5-apx and also the trivial “any step” strategy where we always make the 1-Graver-best step,
which corresponds to taking Γany = {1}. Recall that due to the trick of always exhausting
the discovered direction, this strategy actually has a chance at quick convergence, unlike if
we only made the step with γ = 1.

5.1

Instances

We choose two problems for which N -fold IP formulations were shown in the literature,
namely the Q||Cmax scheduling problem [21] and the Closest String problem [19].
Uniformly related machines makespan minimization (Q||Cmax )
Input: Set of m machines M , each with a speed si ∈ N. A set of n jobs J, each with a
processing time pj ∈ N.
Find:
Find an assignment of jobs to machines such that the time when last job finishes (the
makespan) is minimized; a job j scheduled on a machine i takes time pj /si to execute.
Closest String
Input: A set of k strings s1 , . . . , sk of length L over an alphabet Σ.
Find:
Find a string y ∈ ΣL minimizing maxki=1 dH (si , y), where dH is the Hamming distance.

For both problems we generate random instances as follows.
Scheduling. We view the problem as a decision problem where, given a number B, we
ask whether a schedule of makespan at most B exists. This is equivalent to the multi-sized
bin packing problem, where we have m bins of various capacities instead of m machines
of different speeds, and we adopt this view as it is more convenient. We also view it as a
high-multiplicity problem where the items are not given explicitely as a list of item sizes,
but succinctly by a vector of item multiplicities. Because the algorithm is primarily an
optimization algorithm, we follow the standard approach [15, Lemma 3.8] and turn the
feasibility problem into an auxiliary optimization instance where finding a starting feasible
solution is easy. Specifically for Q||Cmax this means introducing auxiliary slack variables for
“not-yet-scheduled” jobs and minimizing the total “not-yet-scheduled” length.
The input parameters of the instance generation are number of bins m, the smallest
and largest capacities S and L, respectively, item sizes p1 , . . . , pk and probability weights
Pk
w1 , . . . , wk , W = i=1 wi , and a slack ratio σ. The instance is then generated as follows.
First, we choose m capacities from [S, L] uniformly at random. This determines the total
capacity of the bins C. Our goal is to generate items whose total size is roughly σ · C.
We do this by repeatedly picking an item length from p1 , . . . , pk , where pj is selected with
probability wj /W , until the total size of items picked so far exceeds σ · C, when we terminate
and return the generated instance.

23:9

23:10

Evaluating and Tuning n-fold Integer Programming

Closest String. As before, we view the problem as a decision problem: given a number
d ∈ N, decide whether there is a string y with maxki=1 dH (si , y) ≤ d. The random instance is
generated exactly as done by Chimani et al. [5]: first, we generate a random “target” string
y ∈ ΣL and create k copies s1 , . . . , sk of it; then, we make α random changes in s1 , . . . , sk .
This way, we have an upper bound α on the optimum. The input parameters of the instace
generation are thus k, L, Σ and the distance ratio r such α = n/r. Again, we solve an auxiliary
instance where we essentially start with a string of “all blanks” and try to fill in all the
blanks while staying in the specified distance d; the objective is thus the remaining number
of blanks.

5.2

Results

Here we demonstrate the overall behavior of the algorithm on three selected instances; we encourage the reader to see the full data (incl. plots) at https://github.com/katealtmanova/
nfoldexperiment. We have determined that the sensible range of values of g1 is between
2 and roughly 100 and beyond that the behavior does not change significantly. For all
augmentation strategies there are values of g1 which take much longer to converge than any
other values; because of these outliers we clip our figures. To expose the behavior of the
algorithm we use two types of figures.
Outer loop. In the outer loop figure we focus on the behavior of the algorithm with respect
to the loop starting at line 2 of Algorithm 2, i.e., where each iteration corresponds to making
an augmenting step. There is a line plot for each tested value of g1 . The x axis shows the
iteration, the y axis shows the objective value attained in this iteration, i.e.,wxi for iteration
i. We indicate the expensiveness of computing one augmentation by the thickness of the line
in a given iteration – the thicker the line, the more times the (AugILP) has been solved in
this iteration. The legend indicates the exact number of times (AugILP) has been solved for
this value of g1 .
Inner loop. In the inner loop figure we focus on the loop starting at line 4 of Algorithm 2,
i.e., where iterations correspond to solutions of (AugILP). As before, each color corresponds
to a tested value of g1 . There is a line plot displaying the minimum over augmenting steps
found in each outer iteration; however now there is also a semiopaque region above this line,
indicating the values of all the augmenting steps (including the non-minimal ones) found in
this iteration.
We chose two instances among the tested once as representative of the overall behavior:
Our first instance is Q||Cmax with parameters m = 15, S = 2000, L = 10000, item sizes
(2, 3, 13, 35) (so that we have an instance with nontrivial kAk∞ ), and weights (6, 13, 2, 1)
and σ = 0.45. The theoretical upper bound on g1 (A) is (rskAk∞ + 1)O(rs) , and here we
have r = 4, s = 1 and kAk∞ = 35; thus, without computing g1 (A) exactly, we should
consider it to be at least (4 · 36)4 . See Figures 1 and 2.
The third instance is Closest String with parameters k = 5, |Σ| = 2, L = 10000 and
r = 1. The N -fold model has r = 5, s = 1 and kAk∞ = 1, thus, without computing g1 (A)
exactly, we should consider it to be at least (2 · 5)5 . See Figure 3.

5.3

Conclusions

Our main takeaway regarding Question #1 is that, while the theoretical upper bounds
for g1 (A) are huge, already small values of g1 (g1 > 5 for Closest String and g1 > 20

K. Altmanová, D. Knop, and M. Koutecký

35000
30000

objective

25000
20000
15000
10000

35000
30000
25000
20000
15000
10000

5000

5000
5

10

15

20
iteration

25

30

35

40

g1=2 AugILP: 1
g1=3 AugILP: 1
g1=4 AugILP: 5
g1=5 AugILP: 15
g1=6 AugILP: 15
g1=7 AugILP: 11
g1=8 AugILP: 16
g1=9 AugILP: 16
g1=10 AugILP: 12
g1=12 AugILP: 12
g1=14 AugILP: 5304
g1=16 AugILP: 111
g1=18 AugILP: 127
g1=20 AugILP: 114
g1=23 AugILP: 129
g1=26 AugILP: 107
g1=29 AugILP: 128
g1=33 AugILP: 106
g1=40 AugILP: 88
g1=45 AugILP: 81
g1=50 AugILP: 74
g1=60 AugILP: 58
g1=75 AugILP: 50
g1=90 AugILP: 48
g1=110 AugILP: 38

40000
35000
30000
25000
20000
15000
10000
5000

0

5

10

15

20
iteration

25

30

35

40

g1=2 AugILP: 1
g1=3 AugILP: 1
g1=4 AugILP: 2
g1=5 AugILP: 4
g1=6 AugILP: 4
g1=7 AugILP: 4
g1=8 AugILP: 4
g1=9 AugILP: 6
g1=10 AugILP: 4
g1=12 AugILP: 6
g1=14 AugILP: 5765
g1=16 AugILP: 40
g1=18 AugILP: 127
g1=20 AugILP: 257
g1=23 AugILP: 96
g1=26 AugILP: 119
g1=29 AugILP: 1864
g1=33 AugILP: 99
g1=40 AugILP: 36
g1=45 AugILP: 37
g1=50 AugILP: 30
g1=60 AugILP: 80
g1=75 AugILP: 30
g1=90 AugILP: 44
g1=110 AugILP: 52

40000
35000
30000
25000
objective

0

objective

g1=2 AugILP: 1
g1=3 AugILP: 1
g1=4 AugILP: 10
g1=5 AugILP: 25
g1=6 AugILP: 25
g1=7 AugILP: 25
g1=8 AugILP: 25
g1=9 AugILP: 24
g1=10 AugILP: 24
g1=12 AugILP: 24
g1=14 AugILP: 2379
g1=16 AugILP: 233
g1=18 AugILP: 228
g1=20 AugILP: 232
g1=23 AugILP: 237
g1=26 AugILP: 207
g1=29 AugILP: 193
g1=33 AugILP: 149
g1=40 AugILP: 139
g1=45 AugILP: 119
g1=50 AugILP: 127
g1=60 AugILP: 109
g1=75 AugILP: 94
g1=90 AugILP: 77
g1=110 AugILP: 69

40000

objective

g1=2 AugILP: 1
g1=3 AugILP: 1
g1=4 AugILP: 8
g1=5 AugILP: 119
g1=6 AugILP: 156
g1=7 AugILP: 195
g1=8 AugILP: 231
g1=9 AugILP: 273
g1=10 AugILP: 307
g1=12 AugILP: 386
g1=14 AugILP: 13161
g1=16 AugILP: 4449
g1=18 AugILP: 5167
g1=20 AugILP: 5580
g1=23 AugILP: 5937
g1=26 AugILP: 6013
g1=29 AugILP: 6290
g1=33 AugILP: 3747
g1=40 AugILP: 3616
g1=45 AugILP: 3210
g1=50 AugILP: 3093
g1=60 AugILP: 2383
g1=75 AugILP: 2186
g1=90 AugILP: 2122
g1=110 AugILP: 1786

40000

23:11

20000
15000
10000
5000

0

5

10

15

20
iteration

25

30

35

40

0

5

10

15

20
iteration

25

30

35

40

Figure 1 Outer loop results for Makespan Minimization when using (left to right) Γbest , Γ2-apx ,
Γ5-apx , and Γany ; clipped to 40 outer iterations.

for Makespan Minimization) are sufficient for convergence. We remark that, in the
case of Closest String, this hints at the possibility that the maximum value of any
feasible augmenting step g ∈ G(A) is bounded by k O(1) rather than k O(k) , which would
imply an algorithm with runtime k O(k) log L while the currently best algorithm runs in time
2
k O(k ) log L [19].
Regarding Question #2, we see that Γ2-apx provides essentially the same convergence
rate as Γbest but is orders of magnitude cheaper to compute. The “any step” augmentation
strategy Γany usually converges surprisingly quickly, but our results make it clear that its
behavior is erratic and unpredictable. The inner loop Figure 2 reveals that a good step
(close to a Graver-best step) is usually found for larger step-length γ; this motivates adding
the Γ5-apx augmentation strategy to the comparison, as it spends less time on short steplengths than Γ2-apx . Figure 1 shows that using 5-approximate Graver-best steps instead of
2-approximate does not affect the outer loop convergence much, and Figure 2 shows that in
terms of the total number of (AugILP) calls it performs better.
Furthermore, we observe that solving (AugILP) using a MILP solver such as Gurobi
typically takes essentially as much time as solving (ILP) itself; in other words, current MILP
solvers are (without tuning) unable to make any use neither of the extra structure of E (N ) ,
nor the fact that we are seeking a solution with small `1 norm and the right hand side
is 0. Moreover, with a growing number of bricks N , the time to solve (AugILP) using a
MILP solver grows superlinearly, suggesting that, for large enough N , a specialized dynamic
programming algorithm might be competitive with generic MILP solvers.

Evaluating and Tuning n-fold Integer Programming

35000

30000

objective

25000

20000

15000

g1=2 AugILP: 1
g1=3 AugILP: 1
g1=4 AugILP: 8
g1=5 AugILP: 119
g1=6 AugILP: 156
g1=7 AugILP: 195
g1=8 AugILP: 231
g1=9 AugILP: 273
g1=10 AugILP: 307
g1=12 AugILP: 386
g1=14 AugILP: 13161
g1=16 AugILP: 4449
g1=18 AugILP: 5167
g1=20 AugILP: 5580
g1=23 AugILP: 5937
g1=26 AugILP: 6013
g1=29 AugILP: 6290
g1=33 AugILP: 3747
g1=40 AugILP: 3616
g1=45 AugILP: 3210
g1=50 AugILP: 3093
g1=60 AugILP: 2383
g1=75 AugILP: 2186
g1=90 AugILP: 2122
g1=110 AugILP: 1786

35000

30000

25000

20000

15000

10000

10000

5000

5000

20

40

60

iteration

80

100
g1=2 AugILP: 1
g1=3 AugILP: 1
g1=4 AugILP: 5
g1=5 AugILP: 15
g1=6 AugILP: 15
g1=7 AugILP: 11
g1=8 AugILP: 16
g1=9 AugILP: 16
g1=10 AugILP: 12
g1=12 AugILP: 12
g1=14 AugILP: 5304
g1=16 AugILP: 111
g1=18 AugILP: 127
g1=20 AugILP: 114
g1=23 AugILP: 129
g1=26 AugILP: 107
g1=29 AugILP: 128
g1=33 AugILP: 106
g1=40 AugILP: 88
g1=45 AugILP: 81
g1=50 AugILP: 74
g1=60 AugILP: 58
g1=75 AugILP: 50
g1=90 AugILP: 48
g1=110 AugILP: 38

40000

35000

30000

25000

20000

15000

0

20

40

iteration

60

80

100
g1=2 AugILP: 1
g1=3 AugILP: 1
g1=4 AugILP: 2
g1=5 AugILP: 4
g1=6 AugILP: 4
g1=7 AugILP: 4
g1=8 AugILP: 4
g1=9 AugILP: 6
g1=10 AugILP: 4
g1=12 AugILP: 6
g1=14 AugILP: 5765
g1=16 AugILP: 40
g1=18 AugILP: 127
g1=20 AugILP: 257
g1=23 AugILP: 96
g1=26 AugILP: 119
g1=29 AugILP: 1864
g1=33 AugILP: 99
g1=40 AugILP: 36
g1=45 AugILP: 37
g1=50 AugILP: 30
g1=60 AugILP: 80
g1=75 AugILP: 30
g1=90 AugILP: 44
g1=110 AugILP: 52

40000

35000

30000

25000
objective

0

objective

g1=2 AugILP: 1
g1=3 AugILP: 1
g1=4 AugILP: 10
g1=5 AugILP: 25
g1=6 AugILP: 25
g1=7 AugILP: 25
g1=8 AugILP: 25
g1=9 AugILP: 24
g1=10 AugILP: 24
g1=12 AugILP: 24
g1=14 AugILP: 2379
g1=16 AugILP: 233
g1=18 AugILP: 228
g1=20 AugILP: 232
g1=23 AugILP: 237
g1=26 AugILP: 207
g1=29 AugILP: 193
g1=33 AugILP: 149
g1=40 AugILP: 139
g1=45 AugILP: 119
g1=50 AugILP: 127
g1=60 AugILP: 109
g1=75 AugILP: 94
g1=90 AugILP: 77
g1=110 AugILP: 69

40000

objective

40000

20000

15000

10000

10000

5000

5000

0

20

40

60

iteration

80

100

0

20

40

iteration

60

80

100

Figure 2 Inner loop results for Makespan Minimization when using (left to right) Γbest , Γ2-apx ,
Γ5-apx , and Γany ; clipped to 100 inner iterations.

1500

2000

1500

2000

1500

1000

1000

1000

500

500

500

0

0

5

10

15
iteration

20

25

30

0

0

5

10

15
iteration

20

25

30

g1=2 AugILP: 1
g1=3 AugILP: 11
g1=4 AugILP: 31
g1=5 AugILP: 31
g1=6 AugILP: 39
g1=7 AugILP: 37
g1=9 AugILP: 45
g1=10 AugILP: 41
g1=12 AugILP: 50
g1=14 AugILP: 44
g1=16 AugILP: 51
g1=18 AugILP: 50
g1=20 AugILP: 47
g1=23 AugILP: 41
g1=26 AugILP: 48
g1=29 AugILP: 49
g1=33 AugILP: 53
g1=40 AugILP: 51

2500

objective

2000

g1=2 AugILP: 1
g1=3 AugILP: 81
g1=4 AugILP: 176
g1=5 AugILP: 176
g1=6 AugILP: 151
g1=7 AugILP: 147
g1=9 AugILP: 115
g1=10 AugILP: 103
g1=12 AugILP: 97
g1=14 AugILP: 90
g1=16 AugILP: 81
g1=18 AugILP: 71
g1=20 AugILP: 65
g1=23 AugILP: 69
g1=26 AugILP: 60
g1=29 AugILP: 55
g1=33 AugILP: 43
g1=40 AugILP: 39

2500

objective

g1=2 AugILP: 1
g1=3 AugILP: 378
g1=4 AugILP: 758
g1=5 AugILP: 860
g1=6 AugILP: 634
g1=7 AugILP: 661
g1=9 AugILP: 499
g1=10 AugILP: 515
g1=12 AugILP: 396
g1=14 AugILP: 448
g1=16 AugILP: 378
g1=18 AugILP: 387
g1=20 AugILP: 398
g1=23 AugILP: 410
g1=26 AugILP: 398
g1=29 AugILP: 428
g1=33 AugILP: 278
g1=40 AugILP: 274

2500

objective

23:12

0

0

10

20

iteration

30

40

50

Figure 3 Outer loop results for the Closest String instance when using (left to right) Γbest ,
Γ2-apx , and Γany .

6

Outlook

We have initiated an experimental investigation of certain subclasses of ILP with block
structured constraint matrices. Our results show that, as theory suggests, for such ILPs
a primal algorithm always augmenting with steps of small `1 norm converges quickly. We
close with a few interesting research directions. First, is there a way to tune generic MILP
solvers to solve (AugILP) significantly faster than (ILP)? Second, what is the behavior of
our algorithm on instances other than N -fold IP? For example, how large must be g1 in
order to attain the optimum quickly for standard benchmark instances, e.g. MIPLIB [22]?
Third, the approach of Koutecký et al. [23] suggests that a key property for the efficient
solvability of (AugILP) is a certain “sparsity” and “shallowness” (formally captured by the
graph parameter tree-depth) of graphs related to A; what are “natural” instances with small
tree-depth?

K. Altmanová, D. Knop, and M. Koutecký

References
1

2
3
4

5

6
7

8
9
10
11
12

13
14

15
16
17

18

19

Livio Bertacco, Matteo Fischetti, and Andrea Lodi. A feasibility pump heuristic for general mixed-integer problems. Discrete Optimization, 4(1):63 – 76, 2007. Mixed Integer
Programming. doi:https://doi.org/10.1016/j.disopt.2006.10.001.
Timo Berthold. Measuring the impact of primal heuristics. Operations Research Letters,
41(6):611–614, 2013.
Ralf Borndörfer, Martin Grötschel, and Ulrich Jäger. Planning problems in public transit.
In Production Factor Mathematics, pages 95–121. Springer, 2010.
Lin Chen and Daniel Marx. Covering a tree with rooted subtrees–parameterized and approximation algorithms. In Proceedings of the Twenty-Ninth Annual ACM-SIAM Symposium
on Discrete Algorithms, pages 2801–2820. SIAM, 2018.
Markus Chimani, Matthias Woste, and Sebastian Böcker. A closer look at the closest
string and closest substring problem. In 2011 Proceedings of the Thirteenth Workshop on
Algorithm Engineering and Experiments (ALENEX), pages 13–24. SIAM, 2011.
Emilie Danna, Edward Rothberg, and Claude Le Pape. Exploring relaxation induced
neighborhoods to improve mip solutions. Mathematical Programming, 102(1):71–90, 2005.
Jesus A. De Loera, Raymond Hemmecke, and Matthias Köppe. Algebraic and Geometric
Ideas in the Theory of Discrete Optimization, volume 14 of MOS-SIAM Series on Optimization. SIAM, 2013.
Friedrich Eisenbrand, Christoph Hunkenschröder, and Kim-Manuel Klein. Faster algorithms for integer programs with block structure. arXiv preprint arXiv:1802.06289, 2018.
Elisabeth Finhold and Raymond Hemmecke. Lower bounds on the graver complexity of
m-fold matrices. Annals of Combinatorics, 20(1):73–85, 2016.
Matteo Fischetti and Andrea Lodi. Local branching. Mathematical programming, 98(13):23–47, 2003.
Gurobi Optimization, Inc. Gurobi optimizer reference manual, 2016. URL: http://www.
gurobi.com.
Stefan Heinz, Wen-Yang Ku, and Christopher J. Beck. Recent improvements using constraint integer programming for resource allocation and scheduling. In International Conference on AI and OR Techniques in Constriant Programming for Combinatorial Optimization
Problems, pages 12–27. Springer, 2013.
Raymond Hemmecke. Exploiting symmetries in the computation of graver bases. arXiv
preprint math/0410334, 2004.
Raymond Hemmecke, Matthias Köppe, and Robert Weismantel. Graver basis and proximity
techniques for block-structured separable convex integer minimization problems. Math.
Program., 145(1-2, Ser. A):1–18, 2014.
Raymond Hemmecke, Shmuel Onn, and Lyubov Romanchuk. n-fold integer programming
in cubic time. Math. Program., 137(1-2, Ser. A):325–341, 2013.
John D. Hunter. Matplotlib: A 2d graphics environment. Computing In Science & Engineering, 9(3):90–95, 2007. doi:10.1109/MCSE.2007.55.
Klaus Jansen, Kim-Manuel Klein, Marten Maack, and Malin Rau. Empowering the
configuration-ip-new ptas results for scheduling with setups times. arXiv preprint
arXiv:1801.06460, 2018.
Thomas Kluyver, Benjamin Ragan-Kelley, Fernando Pérez, Brian E Granger, Matthias
Bussonnier, Jonathan Frederic, Kyle Kelley, Jessica B Hamrick, Jason Grout, Sylvain Corlay, et al. Jupyter notebooks-a publishing format for reproducible computational workflows.
In ELPUB, pages 87–90, 2016.
Dušan Knop, Martin Koutecký, and Matthias Mnich. Combinatorial n-fold integer programming and applications. In Proc. ESA 2017, volume 87 of Leibniz Int. Proc. Informatics,
pages 54:1–54:14, 2017.

23:13

23:14

Evaluating and Tuning n-fold Integer Programming

20

21
22

23
24
25
26
27
28
29
30

Dušan Knop, Martin Koutecký, and Matthias Mnich. Voting and bribing in singleexponential time. In Proc. STACS 2017, volume 66 of Leibniz Int. Proc. Informatics,
pages 46:1–46:14, 2017.
Dušan Knop and Martin Koutecký. Scheduling meets n-fold integer programming. Journal
of Scheduling, Nov 2017. doi:10.1007/s10951-017-0550-0.
Thorsten Koch, Tobias Achterberg, Erling Andersen, Oliver Bastert, Timo Berthold,
Robert E Bixby, Emilie Danna, Gerald Gamrath, Ambros M Gleixner, Stefan Heinz, et al.
Miplib 2010. Mathematical Programming Computation, 3(2):103, 2011.
Martin Koutecký, Asaf Levin, and Shmuel Onn. A parameterized strongly polynomial
algorithm for block structured integer programs. arXiv preprint arXiv:1802.05859, 2018.
Andrea Lodi. Mixed integer programming computation. In 50 Years of Integer Programming 1958-2008, pages 619–645. Springer, 2010.
Shmuel Onn. Nonlinear discrete optimization. Zurich Lectures in Advanced Mathematics,
European Mathematical Society, 2010.
David Pisinger and Stefan Ropke. Large neighborhood search. In Handbook of metaheuristics, pages 399–419. Springer, 2010.
Yves Pochet and Laurence A Wolsey. Production planning by mixed integer programming.
Springer Science & Business Media, 2006.
Matthew J. Saltzman. Coin-or: an open-source library for optimization. In Programming
languages and systems in computational economics and finance, pages 3–32. Springer, 2002.
Tommi Sottinen. Operations research with gnu linear programming kit. ORMS, 1020:200,
2009.
The Sage Developers. SageMath, the Sage Mathematics Software System (Version 7.6),
2017. http://www.sagemath.org.

