Math. Program., Ser. A
DOI 10.1007/s10107-014-0761-5
FULL LENGTH PAPER

Automatic Dantzig–Wolfe reformulation of mixed
integer programs
Martin Bergner · Alberto Caprara · Alberto Ceselli ·Fabio Furini · Marco E. Lübbecke · Enrico Malaguti ·Emiliano Traversi

Received: 21 September 2012 / Accepted: 1 February 2014
© Springer-Verlag Berlin Heidelberg and Mathematical Optimization Society 2014

Abstract Dantzig–Wolfe decomposition (or reformulation) is well-known to provide
strong dual bounds for specially structured mixed integer programs (MIPs). However,A preliminary version of this paper appeared in [3].
Martin Bergner was supported by the German Research Foundation (DFG) as part of the Priority Program“Algorithm Engineering” under Grants No. LU770/4-1 and 4-2.
Alberto Caprara: Deceased.
Electronic supplementary material The online version of this article(doi:10.1007/s10107-014-0761-5) contains supplementary material, which is available to authorized users.
M. Bergner · M. E. Lübbecke (B)Operations Research, RWTH Aachen University, Kackertstraße 7, 52072 Aachen, Germany
e-mail: marco.luebbecke@rwth-aachen.deM. Bergner
e-mail: martin.bergner@rwth-aachen.de
A. Ceselli
Dipartimento di Informatica, Università degli Studi di Milano, Via Bramante 65, 26013 Crema, Italye-mail: alberto.ceselli@unimi.it
F. Furini
LAMSADE, Université Paris-Dauphine, Place du Maréchal de Lattre de Tassigny, 75775 Paris, France
e-mail: fabio.furini@dauphine.fr
A. Caprara · E. Malaguti
DEI, Università di Bologna, Viale Risorgimento 2, 40136 Bologna, Italye-mail: enrico.malaguti@unibo.it
E. Traversi
LIPN, Équipe AOC, Université Paris 13, 99 Avenue Jean-Baptiste Clément, 93430 Villetaneuse, France
e-mail: emiliano.traversi@lipn.univ-paris13.fr

123

M. Bergner et al.

the method is not implemented in any state-of-the-art MIP solver as it is considered
to require structural problem knowledge and tailoring to this structure. We provide
a computational proof-of-concept that the reformulation can be automated. That is,
we perform a rigorous experimental study, which results in identifying a score toestimate the quality of a decomposition: after building a set of potentially good candidates, we exploit such a score to detect which decomposition might be useful for
Dantzig–Wolfe reformulation of a MIP. We experiment with general instances fromMIPLIB2003 and MIPLIB2010 for which a decomposition method would not be thefirst choice, and demonstrate that strong dual bounds can be obtained from the automatically reformulated model using column generation. Our findings support the ideathat Dantzig–Wolfe reformulation may hold more promise as a general-purpose tool
than previously acknowledged by the research community.
Keywords Dantzig–Wolfe decomposition · Column generation · Block-diagonal
matrix · Matrix re-ordering · Automatic reformulation · Hypergraph partitioning
Mathematics Subject Classification (2000)

90C11 · 49M27 · 65K05

1 Introduction
Dantzig–Wolfe reformulation (DWR) of mixed integer programs (MIPs) is a computationally very successful approach to produce high-quality solutions for well-structured
discrete optimization problems like vehicle routing, crew scheduling, cutting stock,
p-median, generalized assignment, and many others [28]. The common structure ofthese problems is the bordered block-diagonal form of the constraint matrix of the MIP
formulation, which reflects that otherwise independent subproblems are coupled only
by some linking constraints. This structure usually gives rise to a column generation
based solution approach.
It is generally agreed that DWR needs tailoring to the application at hand and isquite far from being a general-purpose tool: It is the user who does not only know that
there is an exploitable structure present but also what it looks like, and how to exploit it
algorithmically. In particular, in view of the automatic use of general-purpose cutting
planes in all modern MIP solvers, this is unsatisfactory.
1.1 Our contribution
In this paper we give a computational proof-of-concept that the DWR process canbe automated and applied to a general MIP even when the latter seemingly does not
expose the matrix structure for which DWR is classically applied. We perform a suiteof experiments, the results of which can be taken as advice on what empirically constitutes a good (or bad) DWR. Remarkably, a key ingredient—re-arranging a matrix intobordered block-angular form—has been available for a long time. Also automatically
applying DWR to a given structure is a concept that is implemented in several frameworks. However, these two components have not been combined in a MIP context
before. In this paper, we provide the missing link by proposing how to re-arrange a123Automatic Dantzig–Wolfe reformulation

matrix into a structure that is experimentally well-suited for DWR, in the sense thata subsequent column generation approach consistently closes a significant portion of
the integrality gap.
Our main contributions are summarized as follows:– We reveal that the constraint matrix of a general MIP can be re-arranged for DWR
in many different ways, necessitating a way to a priori evaluate the quality of a
re-arrangement;– we perform a rigorous experimental study which results in a proxy measure for thisquality, which we call the relative border area;– besides the classical bordered block-diagonal matrix structure, also a doublebordered block-diagonal (also called arrowhead) matrix is amenable to DWR when
applying a variant of Lagrangian decomposition. The re-arrangement of matrices
into both forms having a good border area can be accomplished via a hypergraphpartitioning algorithm;– for a set of medium sized instances from MIPLIB2003 and MIPLIB2010 our reformulations on average lead to comparable or stronger root node dual bounds w.r.t.
a state-of-the-art MIP solver with default parameter settings (i.e., cutting planes,
preprocessing, etc. enabled);– on a non-negligible fraction of these instances, we could automatically identify
reformulations yielding computational performances that are globally better thanthose of state-of-the-art MIP solvers;– our computational success is based on the observation that the constraint matrix ofa MIP may not originally contain a Dantzig–Wolfe decomposable form, but can be“forced” into such a form in almost all cases;– performance variability has raised great attention in the computational MIP literature lately; we show experimentally that this phenomenon is very pronounced in
the context of DWR.
As we re-arrange matrices with the goal of applying a DWR, we will use the notionsof re-arrangement of a matrix and (Dantzig–Wolfe) decomposition interchangeably,as the former immediately leads to the latter in our context.
1.2 Related literature
For a general background on DWR of MIPs, column generation, and branch-and-price
we refer to the recent survey [28] and the primer [9] in the book [8] that also devotes
several chapters to the most important applications.
There are several frameworks which perform DWR of a general MIP, and handle
the resulting column generation subproblems in a generic way such as BaPCod [27],
DIP [23], G12 [22], and GCG [16]. In all cases, the bordered block-diagonal matrix
structure needs to be known and given to the algorithm by the user. In [16] it isshown that such a generic reformulation algorithm performs well on bin packing,graph coloring, and p-median problems. Tebboth, in his thesis [26], derives somedecomposable matrix structures of a linear program (LP) when it is formulated ina specific modeling language. A similar approach of indicating the matrix structure
via key words in the modeling language is taken in [7,13,22,25] among others. All

123

M. Bergner et al.

proposals have in common that the user, in one way or another, needs to make available
her knowledge about the decomposition to be applied.
Specially structured matrices, like bordered block-diagonal forms, play an important role in several fields, e.g., in numerical linear algebra. Therefore, several proposals
exist to re-arrange the rows and columns of a matrix in order to reveal such forms.
A typical motivation is to prepare a matrix for parallel computation, like for solving
linear equation systems, see, for instance, [2] and the many references therein. The
goal usually is to identify a given number of independent blocks (of almost equalsize) with as few constraints in the border as possible (see below for more formal
statements). Some works like [30] mention the possible use of such re-arrangements
in DWR of LPs, but we know of no actual experiment with MIPs. The method in [12],for speeding up interior point methods, is based on graph partitioning as is ours. An
exact branch-and-cut algorithm for detecting a bordered block-angular structure wasproposed in [4].
Attempts to evaluate the re-arranged matrices in terms of suitability for DWR weredone rudimentarily in the LP case only [26]. We are not aware of any attempts to evaluate the quality of a decomposition in terms of suitability of DWR for the MIP case(which has an undoubtedly larger potential). In fact, it is not even known what characterizes a “good” decomposition in this context, and our paper gives computationally
supported guidance in this respect.
2 Partial convexification and Dantzig–Wolfe reformulations
A sketch of DWR applied to a MIP is as follows (see e.g., [9] for details). Consider ageneral MIP(1)
max ct x : Ax ≤ b, Dx ≤ e, x ∈ Zn−r × Qr .
Let P := {x ∈ Qn : Dx ≤ e}. The polyhedron PI P := conv{P ∩ Zn−r × Qr } iscalled the integer hull w.r.t. constraints Dx ≤ e. In a DWR based on a single blockof constraints Dx ≤ e we express x ∈ PI P as a convex combination of the (finitely
many) extreme points Q and (finitely many) extreme rays R of PIP ; this leads tomax ct x(2)s.t. Ax ≤ b
λq q +
μr r
x=


q∈Q

r ∈R

λq = 1

q∈Q
|Q|

|R|

λ ∈ Q+ , μ ∈ Q+
x ∈ Zn−r × Qr

(3)

where each q ∈ Q and r ∈ R represent vectors encoding extreme points and rays,
respectively, and variables λ and μ correspond to weights in their combination.

123

Automatic Dantzig–Wolfe reformulation

(a)

(b)

(c)

Fig. 1 Matrix forms

It is well-known that the resulting LP relaxation is potentially stronger than thatof (1) when PI P  P [17], in which case the dual bound one obtains from (2)–(3) is stronger than the one from (1). This is a main motivation for performing thereformulation in the first place. This partial convexification w.r.t. the constraints Dx ≤e corresponds to implicitly replacing P with PI P in (1). This can be done, in principle,
by explicitly adding all valid inequalities for PI P to (1). When this is impracticable,
the implicit description is in a sense the best one can hope for.
The reformulated MIP (2)–(3) contains the master constraints Ax ≤ b, the convexity constraint, and the constraints linking the original x variables to the extended λ andμ variables. In general, MIP (2)–(3) has an exponential number of λ and μ variables,
so its LP relaxation needs to be solved by column generation. The pricing subproblem
to check whether there are variables with positive reduced cost to be added to thecurrent master LP calls for optimizing a linear objective function over PI P , so it isagain a MIP. The effectiveness of the overall approach hinges crucially on our ability
to solve this subproblem, either by a general-purpose solver or by a tailored algorithm
to exploit its specific structure, if known.
DWR is usually applied when matrix D has a block-diagonal form as depicted
in Fig. 1a where D i ∈ Qm i ×n i for i = 1, . . . , k (and D is 0 elsewhere). In thiscase, k disjoint sets of constraints are identified: Dx ≤ e partitions into D i x i ≤ ei ,i = 1, . . . , k, where x = (x 1 , x 2 , . . . , x k ), with x i being an n i -vector for i = 1, . . . , k.
Every D i x i ≤ ei individually is then convexified in the above spirit. We denote by
k the number of blocks of the reformulation. A matrix of the form as in Fig. 1b with
Ai ∈ Qm  ×n i , i = 1, . . . , k is called (single-)bordered block-diagonal. For general
MIPs, often enough, constraints are not separable by variable sets as above, and adouble-bordered block-diagonal (or arrowhead) form is the most specific structure we
can hope for, i.e., the constraint matrix of (1) looks like Fig. 1c with F i ∈ Qm i ×n  , i =1, . . . , k, and G ∈ Qm  ×n  . The m  constraints associated with the rows of Ai are called
the coupling (or linking) constraints and the n  variables associated with the columnsof F i are called the linking variables, denoted by x  . We intuitively speak of the kblocks D i , i = 1, . . . , k, and the remaining border. Each of the k groups of constraints
D i x i + F i x  ≤ ei can be convexified. Let P i := {x ∈ Qn i +n  : D i x i + F i x  ≤ ei }.
Denote by PIi P the associated integer hull and, respectively, by Q i and R i the set of
extreme points and rays of PIi P . Then, the resulting DWR reads
max ct x
s.t.

k


(4)
Ai x i + Gx  ≤ b

i=1

123

M. Bergner et al.

(x i , x  ) =

λiq q +



μri r

(∀i)

(5)

r ∈R iq∈Q iλiq = 1

(∀i)

q∈Q i
|Q i ||R i |λi ∈ Q+ , μi ∈ Q+
x ∈Z

n−r

×Q

r

(∀i)
(6)

where each q ∈ Q and r ∈ R is augmented with x  components. This formulation
generalizes the explicit master format [21] handling both, the presence of linking
variables and a generic number of blocks. It is a variant of Lagrangian decomposition ,similarly used in [12,26]. Instead of replacing each linking variable with a copy for eachblock it appears in, and adding constraints that ensure consistency of these copies, thiscoordination is taken care of by the original x variables. Furthermore, by keeping all x ivariables in the master, one can enable several features of general purpose MIP solvers,
like the separation of generic cuts, advanced branching techniques, preprocessing, etc.,
without any additional implementation issues. We use the above formulation in our
experiments. We can further obtain a master problem formulation by projecting out theoriginal variables, and introducing for each linking variable a set of constraints which
ensures consistency in different blocks. However, we remark that the techniques andfindings described in the rest of the paper are not strictly dependent on the particularmaster problem format.
3 Automatic detection of an arrowhead form
Potentially, every MIP model is amenable to DWR, even if its structure is not known
in advance (from the modeler or from other sources). In the latter case, we need todetect a structure algorithmically. We need hence to decide: i) which constraints of
the MIP (if any) to keep in the master problem; ii) the number of blocks k and iii) how
to assign the remaining constraints to the different blocks. In other words, we need topartition the set of the original constraints into one subset representing the master andseveral subsets representing the blocks.
To this end, we follow the ideas proposed in [2] and [12] which build on wellknown connections between matrices and hypergraphs. Precisely, once the number kof blocks is fixed, we propose a procedure for obtaining a decomposition consistingof two main ingredients: i) a systematic way to produce a hypergraph starting fromthe constraint matrix, and ii) an algorithm that partitions its vertex set into k subsets(blocks). The partition obtained corresponds to a re-arrangement of the constraint
matrix of the original problem in such a way that it presents an arrowhead form readyto be used as input for the model (4)–(6).
3.1 Hypergraphs for structure detection
We use two different kinds of hypergraphs defined in the following which are the inputof the partitioning algorithm.

123

Automatic Dantzig–Wolfe reformulation

Row-Net Hypergraph. Given the matrix A, we construct a hypergraph H = (V, R) as
follows. Each column j of matrix A defines a vertex v j ∈ V , and each row i of matrix
A defines a hyperedge ri ∈ R linking those vertices v j ∈ V whose corresponding
variable j has non-zero entry ai j = 0 in the row i.
Row-Column-Net Hypergraph. Given the matrix A, we construct a hypergraph H =(V, R ∪ C) as follows. Each entry ai j = 0 of matrix A defines a vertex vi j ∈ V . For
every row i of A, we introduce a hyperedge ri ∈ R which spans all vertices vi j ∈ Vwith ai j = 0; analogously we introduce a hyperedge c j ∈ C for every column jspanning all vertices for which ai j = 0.
3.2 Hypergraph partitioning
In order to obtain a decomposition we heuristically solve a minimum weight balanced
k-partition problem on one of the two hypergraphs. This problem is to partition thevertices of a hypergraph V into k components Vi (i = 1, . . . , k) such that the sumof weights on hyperedges containing vertices in more than one Vi is minimized. The
partition is said to be balanced if |Vi | ≈ n/k.
For a given k-partition of a hypergraph, the blocks are obtained by grouping together
constraints corresponding to hyperedges ri ∈ R consisting of vertices v belonging to
the same component; in the case of the row-net hypergraph, hyperedges r spanning different components correspond to linking constraints and are hence kept explicitly in themaster problem. Similar, for a row-column-net hypergraph, the hyperedges c spanning
different components correspond to linking variables. We recall that by construction,a partition on a row-net hypergraph does not allow partitions with linking variables.
The idea behind both constructions is to obtain “homogeneous” blocks, i.e., blocks
consisting of constraints on similar sets of variables and limiting at the same time
the number of linking constraints (in both cases) and linking variables (when a rowcolumn-net hypergraph is used). To take into account slightly un-balanced partitions,dummy nodes that are not incident to any hyperedge are included in each graph. The
possibility of assigning weights to hyperedges enables us to penalize different components of the border differently; for example linking variables may be less desirable
when they are discrete, or linking variables may be more or less difficult to handle inthe subsequent column generation process than linking constraints, etc.
3.3 A very special case: the temporal knapsack problem
The arrowhead structure of a matrix is so general that it cannot happen that there isno such structure to detect. Therefore, “detection” can also be understood as “forcing” a matrix into a particular structure. We illustrate this on the temporal knapsackproblem [5], also known as unsplittable flow on a line or resource allocation.
The problem is defined as follows. There are n items, the i-th of which has size wi ,a profit pi , and is active only during a time interval [si , ti ). A subset of items has tobe packed into a knapsack of capacity C such that the total profit is maximized andthe knapsack capacity is not exceeded at any point in time. It suffices to impose that123

M. Bergner et al.

(a)

(b)

Fig. 2 Constraint matrix of an instance of the temporal knapsack problem a original and b “forced” intoarrowhead form with k = 10 blocks

the capacity constraint is satisfied at the discrete points in time si , i = 1, . . . , n. Let
S j := {i : si ≤ j and j < ti } denote the set of active tasks at time j, and xi a binary
variable equal to 1 if task i is selected, a binary program for the problem reads:
max

⎧
n
⎨
⎩

i=1

pi xi :


i∈S j

⎫
⎬

wi xi ≤ C, ∀ j ∈ {s1 , . . . , sn }, xi ∈ {0, 1}, ∀i ∈ {1, . . . , n} .
⎭(7)

Figure 2a shows that the coefficient matrix of this formulation does not contain any
non-trivial blocks (i.e., the entire matrix is one block), although all non-zero entries
in a column are associated with consecutive rows.
Applying the procedure explained in the previous subsection based on a rowcolumn-net hypergraph produces a decomposition with only a few linking variables,
see Fig. 2b, and “reveals” a Dantzig–Wolfe decomposable form. This form immediately suggests how to convexify groups of constraints in the spirit of (2)–(3) to (much)better describe the integer hull implicitly. A DWR of the second matrix does not only
produce a strong relaxation but also allows to quickly solve instances to optimality
by branch-and-price, significantly outperforming a state-of-the-art MIP solver on thestandard formulation (7) [5].
4 Experimenting with a generic decomposition framework
An important lesson we learned from preliminary experiments [3] is that we are notlooking for “the” decomposition of a matrix. There is a much larger degree of freedom than originally expected: an arrowhead form is not unique; changing the input
parameters of our re-arrangement algorithm can very much vary the results, and mostimportantly, seemingly small changes in the resulting decomposition may lead to very
different behavior in the subsequent column generation process (both, in terms of dualbound quality and running time). Such sensitivity is a well-documented phenomenon
in MIP solving in general [20]. Figures 3 and 4 show for instance the influence ofthe number k of blocks and the choice of the weights on hyperedges, respectively, on
the visual shape of the resulting arrowhead form. Not least, this non-uniqueness of
decompositions immediately calls for a priori evaluation criteria for a given decomposition w.r.t. its usefulness for DWR. In this study, we concentrate on approaching
this question experimentally.

123Automatic Dantzig–Wolfe reformulation

(a)

(b)

(c)

(d)

Fig. 3 Influence of the number k of blocks on the shape of the arrowhead form (fiber); a smaller k tends
to give a smaller border

(a)

(c)

(b)

(d)

Fig. 4 Visual differences for different decompositions (resulting from different settings in the partitioning
procedure). Top row (a, b) shows instance 10teams with k = 3 blocks; bottom row (c, d) shows instancerout with k = 4 blocks

In the following we describe the benchmark instances used in our computational
experiments, the specific setting used in our experimental framework (i.e., model (4)–(6) of Sect. 2) and the sets of decompositions achieved using the procedure describedin Sect. 3 using three different settings.

123

M. Bergner et al.

4.1 Benchmark instances
We did not experiment on instances with known structure: in these cases the effectiveness of DWR has already been proven, and indeed, the purpose of our study
is to prove that DWR can work for general MIPs. Instead, in order to confirm thegenerality of the proposed method we tested our algorithm on MIPLIB2003 [1] andMIPLIB2010 [20] instances. We selected a subset of instances, for which (a) the density is between 0.05 and 5 %, (b) the number of non-zeros is not larger than 20,000,
and (c) the fraction of discrete variables is at least 20 %. The rationale behind thischoice is the following: (a) if the instance is not sparse enough, no useful arrowhead
form can be detected, and therefore it is easy to tell a priori that DWR is not a promising option; (b) for large instances the partitioning heuristics may fail in finding goodpartitions, and therefore the potential of our approach cannot be estimated; and (c)without a sufficient degree of integrality no gains in the dual bound are possible from a
decomposition.

4.2 Experimental framework
All experiments to assess the performances of model (4)–(6) are carried out on anIntel CoreTM i7-870 PC (2.93 GHz, 8MB cache, 8GB memory) running Linux 2.6.34(single thread). In all our experiments the CPU time limit is set to 1 hour.
For solving the column generation MIP subproblems we use the branch-and-cut
algorithm implemented in CPLEX 12.2, with single thread and default parametersettings. In order to obtain a generic multiple pricing strategy, each time a new incumbent solution having negative reduced cost is found during the optimization of a MIP
subproblem, the corresponding column is inserted into the master. Whenever a MIP
subproblem is found to be unbounded, its integrality conditions are relaxed, and an
extreme ray is generated and added to the master by solving the resulting LP problem
to optimality.
For solving the master problem LPs we use the simplex algorithm of CPLEX 12.2,
with default parameter settings, again single thread. We implemented the dual variablestabilization method described in [11]. We keep as a stability center the dual solution
giving the best Lagrangian bound so far, and we enforce a penalty factor  whenevereach dual variable takes values outside an interval of width δ around the stability
center. At each column generation iteration we change , independently for each dual
variable, to a value randomly chosen between 0 and 10−4 . Whenever the stability
center changes we set δ := 0.00005 · |z − v|, uniformly for all dual variables, wherez is the current master problem LP value, and v is the best Lagrangian bound foundso far. An upper bound of 50.0 is imposed on δ in any case. The stabilization is activated on problems with more than 400 constraints when, during the column generation
process, the gap |z −v| is between 1 and 80 % of the current |z| value. We experimented
with different parameter settings, but these values gave best performance; we experimentally observed that, on our datasets, this stabilization mechanism is enough toovercome potential convergence problems. As we did not go for efficiency, no further
performance improvement method was implemented.

123

Automatic Dantzig–Wolfe reformulation

4.3 Sets of decompositions
Different decompositions can be obtained using the hypergraphs described in Sect. 3.1
with different sets of input parameters. For instance, completely different decompositions can be derived by changing the number k of blocks, the number of dummy
vertices, and the weights of the hyperedges given as input parameters to the partitioning algorithm.
We do not optimally solve the NP-hard minimum weight balanced k-partition
problem. Instead, we use the heuristic multilevel hypergraph partitioning algorithmin [18], of which the package hMETIS [19] is an implementation. In particular, thehMETIS heuristics follow a multilevel recursive bisection paradigm, working in threephases: coarsening, partitioning, and uncoarsening. Coarsening aims at constructing
a sequence of successively smaller hypergraphs by contracting hyperedges with aso-called “FirstChoice” scheme; then, balanced partitioning is performed on the contracted hypergraph by bisection; finally, during uncoarsening, the contracted hypergraph is expanded by successive projections, following backward the coarsening steps,
and running after each projection a local refinement procedure (FM), that is based ontabu search. See [18] for details. hMETIS can solve the partitioning problem neededin our procedure in a few seconds. All hMETIS parameters were kept at their default
values, except the random seed that was set to 1.
Decompositions without Linking Variables. As first set, we generate decompositions
leading to the standard single bordered structure with no linking variables. This shapeis important because the majority of formulations used in practice for DWR present
this shape. We create a row-net hypergraph imposing a weight of 1 on hyperedges
corresponding to constraints. We add a number of dummy vertices equal to 20 % of
the number of non-zero entries in the constraint matrix. We set k = 2, . . . , 20 for
the number of blocks, obtaining 19 decompositions for each instance. An example
of a matrix in a MIP model, and the corresponding detected structure are reported inFig. 5a and d, respectively.
Balanced Decompositions. In the most general setting we search for decompositions
in which both constraints and variables can be part of the border. We create a rowcolumn-net hypergraph with a weight of 1 on hyperedges corresponding to continuous
variables, a weight of 2 on hyperedges corresponding to integer variables, and a weightof 5 on hyperedges corresponding to constraints. We add a number of dummy vertices
equal to 20 % of the number of non-zero entries in the constraint matrix. We set
k = 2, . . . , 10 for the number of blocks, obtaining 9 decompositions for each instance.
Figure 5b and e show an example matrix as given in a MIP model, and the structure
detected using this approach, respectively.
Decompositions with few Linking Constraints. We created a row-column-net hypergraph imposing a weight of 1 on hyperedges corresponding to continuous variables, aweight of 2 on hyperedges corresponding to integer variables and a weight of 105 on
hyperedges corresponding to constraints. We added a number of dummy vertices equal
to 20 % of the number of non-zero entries in the constraint matrix. We set k = 2, . . . , 10
for the number of blocks, obtaining 9 decompositions for each instance. See Fig. 5c123M. Bergner et al.

(a)

(d)

(b)

(e)

(c)

(f)

Fig. 5 Patterns of non-zero entries (black dots) in the coefficient matrix of a MIPs from the MIPLIB2003 [1];the areas with grey background in (d)–(f) emphasize the embedded structure; a–c show that matrix structure
directly from the MPS file, with a re-ordering of rows and columns according to a decomposition foundwith (d) no linking variables, e balanced settings, f no linking constraints

and f for an example of a MIP model matrix and the structure detected in this way,respectively.
Altogether, we obtain a benchmark sample of around 1500 decompositions of 39
instances, which cover a wide range of different possible shapes. These decomposition
are used for testing our computational framework in the following section.
5 Linking input and output measures
In the previous sections, we explained how the set of decompositions has been createdand solved via DWR. We now offer a complete analysis of the results obtained byconsidering the link between the following output and input measures:
Output Measures. The definition of “usefulness” of a given decomposition is alreadyan open issue. To obtain as much methodological insight as possible, we primarily
measure such a usefulness in terms of the root node dual bound of the obtained master relaxation, being a reliable indicator on the effectiveness of the corresponding
convexification process.
Efficiency is not the primary aim of our investigation. Yet, in order to evaluate adecomposition’s computational potential, we consider as a second usefulness measure
the computing time needed to complete the column generation process at the root
node (including the time needed to obtain the decomposition, which is negligible). We
intuitively call this the computation time of a decomposition.
Input Measures. In an effort to understand which parameters affect most the output
measures, we take into account the following four direct indicators: number of blocks,
percentage of linking variables, percentage of linking constraints, and the percentage
of border area. The procedure explained in Sect. 3 can produce decompositions withdifferent values of the direct indicators, this can be done by explicitly fixing the num-123Automatic Dantzig–Wolfe reformulation

ber of blocks a priori or by properly changing the objective function of the graphpartitioning problem. As second set of indicators we consider the average density ofthe blocks and the average integrality gap of the subproblems. These are called indirect indicators because it is not possible with our procedure to obtain decompositions
minimizing or maximizing them.
Methodology of the Analysis. In each of the following subsections we consider eachof these input measures independently, we indicate the expected behavior from theory
and we present our experimental results, highlighting an eventual correlation to output
measures.
In particular, each subsection synthesizes the results of our experiments in two
figures, having the following structure. Each dot corresponds to one of the decompositions generated by the partitioning algorithm. The values on the horizontal (vertical)
axis correspond to a particular input (output) measure.
The integrality gap closed by a decomposition is computed by 1 − |OPT −DWR|/|OPT − LP|, where OPT is the value of a best integer solution found by
CPLEX in one hour of computing time, DWR is the dual bound found by column
generation using a particular decomposition, and L P is the LP relaxation value ofthe corresponding instance. On a few decompositions column generation could not
be completed within the given time limit; in these cases DWR was set to the best
Langrangian relaxation value found during column generation. We normalized thecomputing time t of every decomposition of each instance h as (t − μh )/σh with μh
being the mean running time and σh being the standard deviation of running time ofall decompositions created for that specific instance. For each correlated input-output
pair, we fit a LOESS curve [6] with a span size of s = 0.75.
We further plot the different types of decompositions with different colors, where
the decompositions with no linking variables are plotted with a black dot (•), balanced
decompositions with a dark grey dot (•) and decompositions with few linking constraints with a light grey dot (•). For every figure in the sequel, we provide a figure
disaggregated by instance in the appendix.
In the end of the analysis we propose a proxy measure, suitable for guiding an apriori choice among several potential decompositions. This proxy measure is used in
Sect. 6, where more computationally related benchmarks are discussed and where wecompare the performances of our decompositions with CPLEX.

5.1 Direct indicators
As sketched above, direct indicators are those that can be either chosen as input
parameters or optimized in the hypergraph partitioning algorithm.

5.1.1 Number of blocks
We first observe the influence of the number k of blocks, that is chosen as a parameterin our partitioning algorithms, on the output measures.

123

M. Bergner et al.

(a)

(b)

Fig. 6 Influence of the number k of blocks on the output measures

Expectation 1 As the number of blocks increases, both the integrality gap closed at
the root node and the computing time decrease.
In fact, on one hand a higher number of blocks means more (smaller) subproblems,that can be optimized faster than a few large ones. On the other hand, more blocks
means a more fragmented convexification, that might lead to looser bounds. We also
remark that, even if we expect the master relaxation to improve as k decreases, there
is no guarantee that a simple reduction in k is useful in terms of bound quality. In
fact, different values of k give different decompositions, in which different regions areconvexified, and no a priori dominance can be established between them. The resultsof our experiments are synthesized in Fig. 6.
Observation 1 The average computing time decreases as the number of blocks
increases.
Observation 2 Decompositions without linking variables usually have smallest computation times among all decompositions per instance.
Observation 3 When the number of blocks is small, it is more likely to close a highintegrality gap.
This can be observed by looking at the distribution of the dots in Fig. 6, and becomes
more evident by looking at the instance-by-instance plots reported in the appendix.
Keeping the number of blocks relatively small is a key point if we want to keep theadvantage of a strong dual bound. From the experiments it is hence also clear that weshould focus on decompositions with no more than 5 or 6 blocks. On the other hand, it
is important to notice that decompositions including 20 blocks exist, in which 100 %of the integrality gap is closed at the root node.
This partially confirms our theoritical expectations: few blocks yield tight bounds,
but as their number increases both good and bad decompositions exist, thus motivating
the search for measures ranking good ones.
5.1.2 Percentage of linking variables
Second, we consider the fraction of linking variables nn in each decomposition (seeFig. 7).

123

Automatic Dantzig–Wolfe reformulation

(a)

(b)

Fig. 7 Influence of the percentage of linking variables on the output measures

Expectation 2 As the fraction of variables in the border increases, the computing
time increases.
In fact, including variables in the border is not reducing the number of variablesin the subproblem, but is only allowing different copies of the same variable to takedifferent values in the subproblems. In a column generation setting, then, different
subproblems may blindly generate solutions which are found to be incompatible at amaster stage, thus slowing down the convergence of the process.
Expectation 3 As the fraction of variables in the border increases, the integrality gap
closed decreases.
As recalled in Sect. 2, including variables in the border is similar to including
constraints: due to the analogy with Lagrangian decomposition ,more variables in theborder lower the potential gain by the reformulation.
Observation 4 The fraction of variables in the border and the integrality gap closed
are linked by weak negative correlation.
In fact, as the fraction of variables in the border increases, the computing timeincreases, and therefore many instances run into time limit, so the integrality gap closed
decreases. Hence, our expectations are confirmed. An opposite trend can be observedfor very low fractions of variables in the border. Indeed, good decompositions exist
for some instances, that rely on including linking variables into the border. For these
instances, the single border decompositions generated by our partitioning algorithms
may be substantially worse than those with a few linking variables, thereby worsening
the average gap closed.
Observation 5 Computing times tend to get worse as the fraction of variables in theborder increases.
The expected unstable behavior is observed also in practice.
5.1.3 Percentage of linking constraints
Third, in Fig. 8 we take into account the fraction of linking constraints mm in eachdecomposition.

123

M. Bergner et al.

(a)

(b)

Fig. 8 Influence of the percentage of linking constraints on the output measures

Expectation 4 As the fraction of constraints in the border increases, the integrality
gap closed at the root node decreases.
The rationale is the following: when a constraint is included in a block, its corresponding region of feasibility is convexified, possibly yielding tighter bounds.
At the same time, in general subproblems are MIPs, and therefore it is hard to tella priori the effect on computing times obtained by adding constraints.
Observation 6 As the fraction of constraints in the border increases, the integrality
gap closed at the root node decreases.
This matches our theoretical expectation.
Observation 7 Computing time is higher for very low percentages of constraints inthe border.
Such an observation is not fully matching a modeler intuition, as when the numberof linking constraints is zero, we may have a pure block-diagonal form, which performs
perfectly in a decomposition setting. In an attempt to explain the experimental phenomenon, we consider two factors. First, the decompositions with very low percentagesof linking constraints include those with largest percentage of linking variables, that
according to Observation 5 yield higher computing time. Second, allowing a relevant
percentage of constraints to stay in the border means having more freedom in thecomposition of blocks. In other words, it facilitates the partitioning algorithm to fit theremaining constraints in the best block. This explanation meets also with practice: inseveral partitioning problems, for instance, a natural decomposition is obtained only
by assigning a relevant fraction of the constraints to the border.
In contrast to our results in Sect. 5.1.2, where having very few linking variables
does not guarantee to provide tight bounds, having very few linking constraints does.
That is, adjusting the fraction of linking constraints experimentally offers more control
on the integrality gap closed with respect to adjusting the fraction of linking variables.

123Automatic Dantzig–Wolfe reformulation

(a)

(b)

Fig. 9 Influence of the relative border area on the output measures

5.1.4 Percentage of border area
Fourth, we investigate the effect of adjusting linking variables and linking constraints
simultaneously, by considering the following relative border areaβ=m  · n + m · n − m  · n
.
m·n

(8)

This is the ratio between the border “area” and the entire matrix “area,” and is a wayof encoding the intuition on visually appealing decompositions sketched in Sect. 4:
having a small border area is what a human expert usually aims to, while trying to makea block diagonal structure appear by selectively choosing border elements. Trivially,
for a fixed number of blocks, a MIP whose constraint matrix is in pure block diagonal
form can be decomposed with border area zero.
Expectation 5 Decompositions having a reduced set of constraints and variables inthe border tend to yield tighter dual bounds.
The rationale comes directly from Expectations 3 and 4: each constraint that appears
in the border corresponds to a constraint that is not convexified; in the same way, eachcolumn that appears in the border requires duplicated variables and linking constraints
to be created and relaxed in a Lagrangian decomposition fashion.
From the theoretical point of view, it is hard to predict the relation between relative
border area and computing times. However, following a modeler intuition, two extreme
cases can be considered: when all the constraints and variables are in the border, theproblem becomes an LP, that can be optimized efficiently; on the opposite, when no
constraints nor variables are in the border, we are facing the original MIP. Therefore,such an observation suggests to investigate on the following with the results of ourexperiments summarized in Fig. 9:
Expectation 6 Decompositions having a reduced set of constraints and variables inthe border tend to yield higher computing times.
Observation 8 A lower relative border area correlates with a better dual bound.

123

M. Bergner et al.

As can be observed, the trend in Fig. 9 is more clear with respect to that of resultson both percentage of linking constraints (Fig. 8) and linking variables (Fig. 7).
Observation 9 Tests whose computing time is high, often correspond to decompositions having relative border area either considerably low or considerably high.
The higher percentage of slow runs with low relative border area meets Expectation 6, and that of slow runs with considerably high relative border area meets Expectation 2 and matches Observation 5. Furthermore, fewer linking constraints implies
more constraints in the subproblems, and hence potentially more extreme subproblem
solutions to generate via column generation, thus yielding higher computing time.
5.2 Indirect indicators
We also consider indirect indicators, that is, measures that cannot be reflected in theobjective function of our heuristic minimum k-way hypergraph partitioning approach.
Therefore, they are not optimized directly during the construction of a decomposition,but only assessed a posteriori.
5.2.1 Average density of the blocks
First, we investigate the influence of the density of subproblems, computed as
δi =

nz i
ni · m i

1
δi ,
k
k

and δ =

i=1

where nz i denotes the number of non-zero coefficients in D i . Each decomposition isevaluated according to the average value δ over its blocks.
Expectation 7 Dense subproblems are more difficult to solve, thus yielding higher
computing time.
This is often observed in the behaviour of MIP solvers. Our computational results
are summarized in Fig. 10.

(a)

(b)

Fig. 10 Influence of the average subproblem density on the output measures

123Automatic Dantzig–Wolfe reformulation
Fig. 11 Correlation between
relative border area and average
subproblem density

Observation 10 As the average subproblem density increases, the integrality gap
closed tends to decrease.
Observation 11 For low values of average subproblem density, as the average subproblem density increases, the computing time tends to decrease, while no particulartrend is observed as the average subproblem density further increases.
In order to explain this phenomenon we compared the average subproblem densityof each decomposition in our experiments with the corresponding relative border area,
see Fig. 11. We observe a strong correlation, and in fact Observations 8 and 10 are
coherent; moreover, the correlation in case of Observations 9 and 11 is reversed, i.e.,lower density decompositions close a large fraction of the integrality gap but takelonger to compute whereas decompositions with a lower relative border area bothtend to have good dual bounds and a lower computation time.
5.2.2 Average integrality gap of subproblems
Finally, we analyze the relation between average integrality gap of the subproblems
and output measures. In order to obtain such an analysis, we considered the MIP
subproblems of each block of each decomposition independently, we set its objectivefunction coefficients to those in the original objective function of the problem, we
measured its LP relaxation value (LP), and then let CPLEX run until proven optimality(that is, in the vast majority of our experiments) or until a time limit of 60 seconds
was reached (value ILP). The integrality gap was estimated as |(ILP − LP)/ILP|, andfor each decomposition the average among its subproblems was computed, excluding
those subproblems for which CPLEX could not find any feasible solution within thetime limit, and those having zero ILP value.
We stress that such a measure, besides being indirect in the sense discussed at the
beginning of the section, requires a potentially non-trivial MIP computation for each
subproblem, and is therefore interesting for a methodological insight only.
Theory suggests that
Expectation 8 A larger integrality gap in the subproblems may imply a larger potential for a stronger master dual bound.
This expectation has also been mentioned in the context of Lagrangean relaxation [17]. Trivially, a subproblem integrality gap of zero for whatever objective func-123

M. Bergner et al.

(a)

(b)

Fig. 12 Influence of the average integrality gap at the root node of the subproblems on the output measurestion indicates that no improvement of the dual bound is possible compared to theoriginal LP relaxation. For the computing time, we expect the following
Expectation 9 The computing time increases as the average integrality gap in the
subproblems increases.
This is due to the fact that closing the integrality gap is the main task of any MIP
solver in tackling each subproblem. Our results are reported in Fig. 12.
Decompositions corresponding to points on the horizontal axis are actually thosetests hitting time limits. In our experiments we found the following.
Observation 12 Both integrality gap closed and computing time correlate weakly toaverage subproblem integrality gap.
Future experiments may reveal the usefulness of a measure related to the subproblem integrality gap, perhaps by considering a more suitable definition of the objectivefunction.
5.3 Evaluating and exploiting a proxy measure
We remark that, when evaluating a decomposition, a strong dual bound cannot constitute the ultimate quality measure, as this can be trivially maximized by reformulating
the whole original problem into a single subproblem (which amounts to optimallysolving the original MIP in the subproblem in a single shot).
From the analysis reported in the previous sections we drew the conclusion that
relative border area is the most appropriate proxy measure to estimate the quality of adecomposition. In particular, it is the only measure exhibiting three appealing featuresat the same time. First, the measure ranges from zero (best) to one (worst), being
zero for block-angular matrices. Second, it correlates negatively with integrality gap
closed. Third, it is eventually positively correlated also with the time needed to computethe bound of the decompositions. In order to be more detailed on this third feature,
we present a further analysis on our results: we rank the decompositions by relative
border area, and we consider four classes of decompositions, corresponding to the four123Automatic Dantzig–Wolfe reformulation
Table 1 Average output
measures for relative borderarea quartiles

Quartile

Range

Norm. int. gap
closedNorm. cpu time1

(0.0031, 0.0546]

71.97

0.10

2

(0.0546, 0.1380]

47.46

−0.07

3

(0.1380, 0.3290]

29.42

−0.15

4

(0.3290, 0.9360]

12.80

0.11

quartiles in such a ranking; in Table 1 we report the average output measures for eachquartile: the table is composed of four rows, one of each quartile, and four columns,
reporting quartile, range of values corresponding to the quartile, average normalized
integrality gap closed, and average normalized computing time over decompositions
in each quartile.
We come to the following
Observation 13 A better tradeoff between integrality gap closed and computing time
can be obtained by looking for decompositions with low relative border area.
In fact, as reported in Table 1, decompositions having relative border area in the
fourth quartile are on average worse than those in all other quartiles in terms of bothcomputing time and integrality gap closed.
For the remaining quartiles, decompositions having low relative border area tend
to provide higher integrality gap closed at the price of higher computing time. Indeed,this matches Expectations 5 and 6, and supports the general intuition that additional
computational effort is needed to obtain better bounds.
In Table 2 we report, for each MIPLIB instance, the details of the DWR withminimum relative border area among all those generated by our graph partitioning
algorithms, as described in Sect. 4.3. Listed are: the instance features (name, numberof rows and columns), number of blocks, linking variables, and linking constraints,
maximum number of variables and constraints in a block, number of LP problems
solved during column generation and corresponding CPU time, number of MIP subproblems solved during the column generation process and corresponding CPU time,number of columns and extreme rays generated during the optimization process. The
table is split in two horizontal blocks, corresponding to MIPLIB 2003 and MIPLIB
2010 instances, respectively. As can be seen, no obvious correlation appears between
the prediction given by the proxy measure and the single features of the chosen decomposition; therefore our proxy measure proves to be unbiased with respect to the otherdecomposition parameters.
Observation 14 Among those generated by our graph partitioning algorithms, adecomposition corresponding to the best measure has either linking constraints or
linking variables, but not both, in the vast majority (36 out of 39) of the instances.
We remark that this phenomenon is not simply induced by Definition (8), as a low
relative border area can also be achieved by including in the border very few variables
and constraints at the same time. That is, our proxy measure seems to match an
intuitive behavior reported in the literature: for each instance, that is the representative

123

123

5,325

mkc

3,411

6,480

128

769

2,756

240

240

556

712

397

675

noswotopt1217

p2756

pp08a

pp08aCUTS

rout

set1chtimtab1

timtab2

294

171

492

291

246

136

755

64

182

3,321

manna81

112

291

2,993

harp2

396

1248

1,392

478

479

422

322

glass4

modglob

1,224

1,224

878

gesa2-o

363

1,298

fiber

fixnet6

gesa2

1,442

842

2,728

aflow40b

2

2

20

5

8

8

3

4

5

2

2

2

10

2

2

2

2

2

10

2

9

aflow30a

230

2,025

10teams

17

13

0

0

0

0

0

0

0

0

0

77

0

12

0

26

0

0

0

0

0

0

0

12

16

8

8

16

16

9

8

31

0

39

0

26

0

14

22

39

28

95

348

206

35

111

30

30

918

192

26

211

2,911

1,854

369

178

611

625

436

713

348

430

225

Var

Var

Cons

MaximalLinking

BlocksCols

Instance

Rows

Best decomposition features

Instance features

148

86

24

55

32

16

257

12

35

144

1,815

3,546

9

212

611

696

235

178

179

230

15

Cons

359

548

53

10

22

14

20

34

31

405

31

5

1,385

201

896

86

20

462

75

899

1,168

#

LP

0.03

0.03

0.02

0

0

0.01

32.86

0.01

0

2,727.36

21.43

300.23

0.24

0.01

7.12

0.53

2,664.52

11.7

55.59

40.82

0.53

Time (s)
468

20

42

260

95

168

160

1,212

120

20

2,768

400

1,790

850

37

922

148

1,796

2,334

3,580

1,094

#

MIP

3,271.48

244.56

0.59

5.9

0.58

0.46

202.090.06

1.26

876.1784.78

234.591.34

2,657.16

106

17.14

173.92528.77

14.89

52.66

1.69

Time (s)

365

228

320

364

319

224

254

5,048

120

31

17,595

1,467

2,962

793

256

3,191

860

3,721

4,658

3,490

2,097

Col

Generated

Table 2 Details of the DWR of minimum proxy measure (8), among those generated by our graph partitioning algorithms, for the corresponding instance

0

0

0

0

12

10

0

0

0

0

0

0

0

1

0

0

0

0

0

0

0

Rays

M. Bergner et al.

3,117rococoC10-001000

670

reblock67

1,096

512

ran16x16

rmine6

1,059

pw-myciel4

3,660

neos-686190

490

900

mine-90-10

pigeon-10

830

mine-166-5

1,747

mcsched

500

2,260

macrophagem100n500k4r1

169

338

1,205

enlight13

gmu-35-40

1,293

7078

2,523

288

8,164

931

3,664

6,270

8,429

2,107

3,164

100

424

351

1,758

csched010

1,750

2,500

beasleyC3

234

8

2

2

16

2

2

3

2

2

5

2

2

2

2

2

2

2

2

378

vpm2

750

1,080

tr12-30

0

0

0

0

0

96

59

0

0

0

7

351

0

25

1

0

7

12

82

354

69

16

165

40

118

92

596

32

0

0

9

0

55

26

0

0

447

557

335

32

530

274

1,324

450

416

349

1,203

440

432

196

823

1,252

196

552

Var

Var

Cons

MaximalLinking

BlocksCols

Instance

Rows

Best decomposition features

Instance features

Table 2 continued

407

3,439

1,228

17

4,009

477

1,265

3,090

3,928

415

1,702

54

217

92

171

862

117

375

Cons

#
2
19

160

35

21

25

59

183

82

23

214

45

8

542

18

35

102

309

LP

3.36

0.21

0.06

0.01

0.12

18.26

75.33

0.03

4.97

0.06

0.03

38.48

0.04

0.05

0.81

2,358.21

0.01

0Time (s)
3

26

68

202

616

36

1,272

68

40

384

114

364

244

44

426

220

16

1,082

#

MIP

10.37

836.68529.66

0.92

27.38

25.85

3,330.84

2,983.29

1,421.13

84.24

3,547.05

48.62

0.08

1,068.67

1,040.41

1,242.35

1.04

0.43

Time (s)

32
140

875

331

299

475

146

1,769

1,548

316

2,278

780

92

2,393

30

190

953

5,219

Col

Generated

0

0

0

0

1

0

0

0

0

0

0

0

4

0

0

0

0

0

Rays

Automatic Dantzig–Wolfe reformulation123

M. Bergner et al.

of a potentially difficult combinatorial optimization problem, Lagrangian relaxation
may be more promising than Lagrangian decomposition, or vice versa. Still, in a few
cases, combining the two techniques may be promising.
Observation 15 Among those generated by our graph partitioning algorithms, adecomposition corresponding to the best measure has no linking variables in the
majority (26 out of 39) of the instances.
This confirms another important modeler’s intuition: DWR unfolds its potential
when relaxing a few linking constraints is enough to decompose the problem in more
tractable subproblems.
6 Comparison to a state-of-the-art MIP solver
Finally, we sketch an answer to a very fundamental computational question: are there
generic MIP instances for which it is worthwhile to perform DWR with an automatically generated decomposition instead of applying traditional cutting-planes methods?
The overall analysis presented in Sect. 5 motivates us to propose the following two“automatic” algorithms.
DWR auto: for each instance, consider all the decompositions generated by the graphpartitioning algorithms as described in Sect. 4.3, pick one minimizing the relative
border area (8), keep the original variables in the master, including the possible linking
ones (as explained in Sect. 2), and perform a DWR of the resulting blocks; then solve
the resulting master problem by column generation.
DWR best: for each instance, consider all the decompositions generated by the graphpartitioning algorithms as described in Sect. 4.3, perform column generation on eachof the resulting master problems, and keep a decomposition giving the best dual bound.
That is, the algorithm “DWR auto” provides a benchmark for assessing the potentialof using the border area for finding decompositions yielding tight dual bounds, whereas“DWR best” aims at assessing the potential of the overall automatic DWR mechanism.
Note that “DWR best” is extremely time consuming and only meant as an assessmentof the method’s potential.
As a benchmark MIP solver we use CPLEX 12.2 with full default settings, inparticular with default preprocessing and cutting planes enabled. As for DWR, the
CPLEX integrality gap closed is computed by 1−|OPT−CPLEX|/|OPT−LP|, where
OPT is the value of the best integer solution found by CPLEX in one hour of computing
time, CPLEX is the dual bound found by CPLEX in the root node, after preprocessing
and cut generation, and L P is the LP relaxation value of the corresponding instance.
6.1 Quality of the dual bound
The first set of experiments aims at assessing the quality of the bound that can beobtained using the DWR approach. Our goal is to show that by using such an algorithm,
it is often possible to obtain a better description of the convex hull of the full problem
than that obtained by CPLEX through a cutting plane process.

123Automatic Dantzig–Wolfe reformulation

Fig. 13 Distribution of fraction of integrality gap closed for all decompositions (MIPLIB2003)

Fig. 14 Distribution of fraction of integrality gap closed for all decompositions (MIPLIB2010)

Figures 13 and 14 depict the distribution of fraction of integrality gap closed for
all the about 1500 decompositions (see Subsect. 4.3) we tried for MIPLIB2003 and
MIPLIB2010 instances, respectively. Every dot represents one decomposition. On
some instances (glass4, modglob, timtab2, tr12-30, beasleyC3, macrophage, neos686190), column generation could not be completed within the time limit, and therefore
the gap of the best Lagrangian bound (potentially worse than the LP relaxation bound),obtained during the column generation iterations, is reported [17]. One can see that adecomposition with the respective best proxy measure (marked by a ×) most of thetime leads to a very good choice among all decompositions we tried with respect tothe root node integrality gap closed, also compared to the bound obtained by CPLEX
after the cutting plane process at the root node (marked by a ).
In Table 3 we detail, for each instance in the dataset, our results for both the “DWRauto” and “DWR best” methods, compared to the results obtained by CPLEX when
stopped after the cutting plane process at the root node as well as the time CPLEX
needs to solve the original problem to optimality.
Listed are: the instance name, the value of an optimal solution or best solution
computed by CPLEX in 1 hour (opt*), the value of the LP relaxation (LP) and thecorresponding duality gap (LP gap). Then, three vertical blocks follow, corresponding

123

123

1005.66

33.33
0.61

25779900.00

1200010000.00 800002400.00 33.33

gesa2-o1077.5654537.80

764772.00

routtimtab15480.61pp08aCUTS 7350.00

set1ch2748.35

7350.00

pp08a−20.0228694.00

32007.73

981.862688.75−16.00−41.00

noswot3124.00−43.0020740500.00

modglob

opt1217−611.8520430947.62−563.85

mkc

p2756

1.01

−13297.00

96.25

41.31

8.88

25.43

62.61

13.93

25.13

4.88

1.49

8.51

0.61

−73899800.00 −74353341.50−13164.00

harp2

1.18

manna8125476489.68

0.00

16.40 82.96

0.10 99.76

0.68 92.34

2.58 89.86

2.50 96.01

0.27 98.05

25.13

0.49 90.00

0.79 47.30

0.16 98.13

0.08 92.05

0.00

0.00

0.78 33.53

0.03 97.03

18.89 72.96

glass4

1.18

69.85

3983.00

25476489.68

0.00

1.07 98.27

13.90

25779900.00

1200.88

13.90

61.55

2.59
1.07

0.49

0.12

0.16

0.04

0.61

8.33

0.05

0.03

0.07

0.68

2.27

2.27

0.27

244.59 14.60

0.61

5.90

0.58

0.47

234.980.07 25.13

1.26

3600.00

106.25535.491.59

3600.00

113.1417.67

2838.48 18.83

540.5370.52 13.90

93.51 10.23

1.46

84.84

99.84

92.34

91.09

96.38

98.07

0.00

90.00

0.00

0.21

0.10

6.06

1.89

6.47

5.35

106.25

4.41

0.33

3.78

0.92

8.86

3.82

2.52

7.47

0.00

4.88

0.14

326.55 39.05

0.88

5.90

9.72

5.74

111.69

0.02

1.26

0.37
0.00

0

0

23,1100020,646

0

58

0

29

0

0

59.43 0.15

97.77 0.07

0.26 0.06

84.97 0.08

95.97 0.08

46.41 0.07

100.00 0.07

38451

511

8921

224

4

199

0

0.00 0.01 69,02,434

90.48 0.03

55.63 0.32

100.00 0.11

39.60 0.05

0.00 0.02

82.38 0.08

91.51 0.06

91.32 0.47

96.92 0.04

53.44 4.15

64.55 1.93

100.00 0.47

Ratio

50.04 4.89

0.34 1.79

8.74 0.68

0.23 2.52

0.09 5.22

0.20 1174.90

0.01 7.00

731.62 0.00

0.01 360356.00

43.86 2.42

0.11 4868.09

0.03 53.00

10.37 346.86

0.03 3771.33

0.13 135.92

0.02 141924.00

0.08 6756.63

0.04 1763.00

0.04 2337.75

0.21 10.57

Time

CPLEX default

Gap (%) Closed Time Nodes(%)

CPLEX root

3.79 33.33

183.6426.75

292.51

540.53550.81

978.9891.74 3,501.33

98.13

96.23

0.00

75.00

96.16

97.82

73.05

98.27

0.00

32.27

100.00Gap (%) Closed Time(%)

DWR best gap

2.22 0.00

Closed Time(%)

0.05 92.86
14.71

fixnet6

156082.52

0.76

15.10

Gap
(%)

LP gap DWR auto gap
(%)

gesa2

1168.00

405935.00

aflow40b

fiber

983.17917.00924.001158.00

10teamsLP

Opt*aflow30a

Instance

Table 3 Comparison of the root node dual bounds provided by DWR and CPLEX: in particular we report the percentage of instances for which DWR was faster than CPLEX(<1), within one order of magnitude (<10), and within two orders of magnitude (<102 )M. Bergner et al.

13.75

vpm2

10.00

3116.43−39339910.92−462.313823.00−34627142.08−457.19

ran16x16

reblock67

rmine60.00−10000.00−9000.00

pigeon-10

pw-myciel4

11.11

1.12

13.61

18.48

100.000.00

0.42 62.23

0.13 99.02

10.77 41.71

60.00 40.00

11.11

0.00

5134.816730.00

neos-686190

23.70

0.49 96.24−784302337.63 −887165318.51 13.12

mine-90-10
23.70

0.06 99.29
7.33 83.74

8.56

211913.00−566395707.87 −821763677.67 45.09

193774.75

0.00

0.00

0.37 99.63

4.17

0.01

100.004.17

0.01

7.91 57.28
57.55 42.45

mcsched0.00

374.00macrophage

18.52

100.00mine-166-5−25.00m100n500k4r1 −24.00−2406943.56−2406677.75

gmu-35-40

332.420.00

408.0071.00

csched010

enlight13754.00<102
beasleyC3
15.38 83.74

11.29 62.53

5.59

1.71

2.83

7.91
0.01

0.49

7.33

0.06

0.05

0.00

836.89529.730.42

0.13

0.93 10.77

27.50

44.13 11.11

3600.00 17.60

2983.32

1426.1284.30

3600.00

87.10 4.17

0.12

1068.72 57.55

1041.233600.00 15.38

1.05

3600.00

3600.00 25.26
1.05

8.24

7.75

6.44

0.68

0.02

0.09

4.17

0.01

84.30

8.54

1.32 11.11

4.84 21.54

8.73

0.93
62.23 836.89

99.02 529.73

41.71

0.98

8.87

6.93

100.00 3600.00 53.54

0.00

25.73

96.24 2983.32

83.74 1426.12 22.79

99.29

99.95 1385.91 11.48

0.00

0.00

42.45 1068.72 80.66

57.28 1041.23 12.62

83.74 3600.00

77.16

93.92

96.83 839.54

0

1,632

0

32,917

8,522

0

0

17828

25,753000

022,078149

12.26 0.60

1,547

34.86 0.63 48,9311

62.51 0.61

46.46 1.87

0.00 0.01

9.11 0.57

33.41 0.85

49.44 1.81

0.23 1.03 110,052

88.52 17.31

0.00 0.06

2.54 0.04

19.34 1.37

31.86 0.55

91.29 1.06

68.57

77.05 0.00

99.24 0.59

50.20 0.31

Ratio

0.46 2.28

0.02 179332.00

85.96 41.64

Time

CPLEX default

Gap (%) Closed Time Nodes(%)

72.66 2178.64 46.00

Gap (%) Closed Time(%)

CPLEX root

9.70 86.28

488.15 1.09

0.02 46.50

0.09 305.56

0.02 2206.50

0.08 44934.25

59.55 50.10

2.18 654.18

367.61 0.23

142.72 25.22

0.02 4355.00

0.03 4.00

149.47 7.15

165.15 6.30

52.17
0.38 9475.16

8.70

94.64

0.00

1.71 93.92

89.12

36.50 60.49

Closed Time(%)

DWR best gap

39.13

30.28

28.08

89.12

92.38

Gap
(%)

LP gap DWR auto gap
(%)

<10

40.43

9.89

14210.43

83592.00

LP<1

Average

1096560.00

130596.00

timtab2

tr12-30

Opt*

Instance

Table 3 continued

Automatic Dantzig–Wolfe reformulation

123

123
10.46

34.42
55.71

0.00

3.87
16.99

11.62
34.26

66.25

Closed
(%)
0.58

Time

0

Nodes

0.07

Time

196.43

Ratio

56.25

50.33

13.75

Gap (%)CPLEX default

31.25

14.61

36.66

0.00

Time

CPLEX root

<102

34.43

34.42

Closed(%)

Gap (%)

Time

Gap
(%)

Closed(%)

DWR best gap

DWR auto gap

6.25

7515.27

LP gap(%)<1011460.00

rococoC10001000
Average

LP

<1

Opt*

Instance

Table 3 continued

M. Bergner et al.

Automatic Dantzig–Wolfe reformulation

to “DWR auto,” “DWR best,” and CPLEX, respectively. For each we report the integrality gap, the improvement with respect to the LP relaxation bound (a value of zerois reported if the decomposition is not able to improve over the LP relaxation bound),the number of nodes used to attain that bound, and the time needed to compute thatbound. In bold we report the best configuration in terms of integrality gap.
The table is composed of two horizontal blocks, that refer to the 23 selected instancesof MIPLIB2003 and to the 16 selected instances of MIPLIB2010, respectively. The
last row of each block reports average values.
On 15 out of the 23 MIPLIB 2003 and 11 out of the 16 MIPLIB 2010 instances,the dual bound found by our DWR approach improves on CPLEX’s root node boundwith default settings, and in four more instances the bound is the same.
On the average, whereas “DWR best” clearly outperforms CPLEX, “DWR auto” isstill competitive, experimentally supporting the meaningfulness of our proxy measurefor the quality of a decomposition.

6.2 Overall performance comparison
Since DWR and CPLEX produce bounds of different quality with different computational efforts, in a second set of experiments we aim at comparing the trade-off
between computing time and quality of the dual bound given by the two methods.
This is an overall index designed to measure the potential of both methods in actuallysolving MIPs to optimality, and to help answering the ultimate question of whether
good dual bounds can be provided by DWR in reasonable time.
We remark that, apart from the computation of bounds, the remaining machinery inbranch-and-price and branch-and-cut algorithms is equivalent, as similar preprocessing, heuristics, constraint propagation, branching techniques, etc., can be implemented.
Also the re-optimization process in the nodes of the branch-and-bound tree is similar,
provided the same decomposition is kept along the search tree. In fact, state of the art
branch-and-price frameworks [16] keep pools of columns, and therefore the optimization of a node is likely to start from a restricted master problem containing at least thecolumns forming an optimal basis for the father node, as in standard branch-and-boundcodes, and very often many more high quality columns already generated during the
optimization of siblings.
At the same time, re-implementing the techniques included in a state-of-the-art
branch-and-cut solver was beyond of the scope of this paper. Therefore, in order toperform a comparison, we decided to consider the DWR giving best proxy measure in“DWR auto,” and to take as a performance index the ratio between the time requiredby our algorithm to obtain the dual bound at the root node, when using such a decomposition, and the time required by CPLEX with default settings to obtain, either at theroot node or through branching, the same bound. In the comparison we included also
“DWR best.”
The results are also reported in Table 3. The last vertical block for CPLEX withdefault settings shows the number of nodes and the CPU time needed to reach thesame bound as “DWR auto;” the ratio between “DWR auto” time and CPLEX time forobtaining such a bound. Below each vertical block we report the fraction of instances123

M. Bergner et al.

(a)

(b)

Fig. 15 Performance profiles comparing DWR auto, DWR best and CPLEX

for which “DWR auto” was faster than CPLEX (<1), within one order of magnitude
w.r.t. CPLEX (<10) and within two orders of magnitude (<102 ).
On average, CPLEX is much faster in solving the root node relaxation. At the sametime, on many instances in which good decompositions can be found, CPLEX needs
to explore many branching nodes, and to spend a high CPU time for matching thebound given by “DWR auto.” Finally, it is interesting to note that in more than halfof the instances (up to 52.17 % of MIPLIB2003 and up to 56.25 % of MIPLIB2010),“DWR auto” is still within a factor 102 from CPLEX, making us optimistic that, after
a suitable software engineering process, our approach would become computationallycompetitive on more instances.
6.3 Performance profiles
Finally, we compared “DWR auto” and “DWR best” to CPLEX using the methodologyof performance profiles [10]. We performed two different analyses, the first one is basedon time performance and the second one is based on bound performance, disregarding
the time needed to compute it.
The time analysis is displayed in Fig. 15a, where we show two different comparisons, the first one concerns “DWR auto” and CPLEX (black lines, dashed and solid),
the second one concerns instead “DWR best” and CPLEX (grey lines, dashed andsolid). For each pair of algorithms and for each instance, we consider as a performance index the ratio between the time needed by the “DWR” method or CPLEX tocompute the “DWR” bound (respectively “auto” or “best”) and the smaller between
the two values. Normalization is then performed with respect to the fastest algorithm.
For each value π on the horizontal axis, we report on the vertical axis the fraction of
the dataset for which the corresponding algorithm is at most π1 times slower than thefastest algorithm. The value on the vertical axis, corresponding to π = 1, indicates thefraction of the instances in which the considered method is the fastest among the two
considered. As far as “DWR auto” is concerned, it outperforms CPLEX in around 8 %of the instances, while “DWR best” in around 19 %. From the figure, we can also seethat in around 50 % of the instances “DWR auto” is at most two orders of magnitude
slower (see also discussion of Table 3 in Sect. 6.2), while “DWR best” only one orderof magnitude.

123

Automatic Dantzig–Wolfe reformulation

In Fig. 15b, we display the second performance profile based on the bounds of thethree methods. As a performance index we consider, for each instance and for eachmethod, the ratio between the bound obtained by such a method and the best boundobtained using either “DWR auto,” “DWR best,” or CPLEX. Normalization is thenperformed with respect to the index of the best algorithm among the three considered.
The picture reports the normalized performance index on the horizontal axis; for eachvalue π on the horizontal axis, we report on the vertical axis the fraction of the dataset
for which the corresponding algorithm closes at least a fraction π of the best achieved
gap, on the considered set of 39 instances. We always assume that the LP relaxation
value is available, thus the minimum closed gap for each instance and each algorithm is
set to 0. For instance, if we look for a method which is able to close at least about 20 %of the best achieved gap, then CPLEX is the best choice. However, if we increase thisrequirement to about 45 % both “DWR auto” and “DWR best” outperform CPLEX.
Finally, in almost 70 % of the cases “DWR best” provides best bound, while both
CPLEX and “DWR auto” give best bound in about 39 % of the instances.

7 Discussion
We have performed the first systematic investigation with an automatic Dantzig–Wolfe
type reformulation of arbitrary MIPs. Even though it is clear from theory that sucha reformulation can be used to improve the dual bound, it has not been considered agenerally useful computational tool in practice. Thus, the most unexpected outcomeof our study is that already a fairly basic implementation, combined with an automatic
choice of the decomposition, is actually capable of competing with a state-of-theart MIP solver. For some instances the dual bound computed by our method is sostrong that, given a heuristic solution of optimal value, optimality can be provenat the root node. Furthermore, on a relevant subset of MIPLIB instances, we could
automatically detect Dantzig–Wolfe type reformulations for which the decomposition
approach yields an overall better computing behavior than a state-of-the-art branchand-cut based general purpose solver.
The results even improve if we choose the decomposition by explicit computations,demonstrating that there is still potential. It turned out that different decompositions forthe same instance lead to, sometimes significantly, different dual bounds (see Figs. 13
and 14), and also to drastic differences in the computation times needed to solve theresulting relaxation. Thus, out of the many questions spawned by our work, the mostimportant one is, both from a theoretical and a practical point of view, to characterize agood decomposition. We believe that answers will be hard to find as they immediately
relate to the very core of computational integer programming: to better describe, in acomputational and efficient manner, the convex hull of integer feasible points. On theother hand, approaching this important topic from a decomposition point of view may
yield new insights previously overlooked.
Our experimental setup for detecting a matrix structure certainly can be improved;
it is just one out of probably many conceivable approaches to come up with a proof-ofconcept. In particular, the process of generating decompositions by simply changing
parameters of partitioning algorithms, and then selecting one with best proxy measure,123

M. Bergner et al.

(a)

(b)

(d)

(e)

(g)

(h)

(c)

(f)

(i)

Fig. 16 Detected matrix structures for selected MIPLIB2003 and MIPLIB2010 instances

can give place to better and more sophisticated approaches. Many experimental issues
on the design of good proxy measures are also left open, e.g., how the blocks’ balancing
and density impact on the overall performances.
Certainly, we will see alternatives in the future. We also alert the reader that theseeming omnipresence of arrowhead structure in MIPLIB instances (c.f. Fig. 16) may
either reproduce structures which were incorporated in the model by a human modeler,
accidentally or on purpose, or simply be an artifact of the model and algorithm we use
to detect/enforce this structure. In any case, only a part of variables and constraints
describe the logic and mechanism of the problem to be modeled. Another substantial
part is present only because of “technical purposes” and “modeling tricks.” Detecting
and exploiting this information in a decomposition may lead to new insights into how
a good MIP model should be formulated.
There are some possible immediate extensions concerning the implementation.
Only further experimentation can show whether the advantage in the root node can beretained throughout the search tree (it is also conceivable that an advantage becomesvisible only further down the tree). If one is only interested in a strong dual bound, theaddition of generic cutting planes is a natural next step (see [16]).
Of course, at the moment, our work is not intended to produce a competitive tool,
but to demonstrate that the direction is promising. Even in the future, we do not expectthat decomposition techniques will become the single best option approach to solve
MIPs. However, we hope that one can soon distinguish a priori, only based on the123Automatic Dantzig–Wolfe reformulation

instance, whether it can pay to apply a decomposition like DWR or not. Our results
indicate that promising instances are more than an exception.
Taking into account the fact that state-of-the-art solvers make successful use ofgeneric cutting planes for about 15 years now, it is clear that outer approximations of
the integer hull have a prominent headway in experience over inner approximations.
We hope to have inspired further research and experimentation with the second option;indeed, follow-up research [14,15,24,29] is already available.
A final word is in order on what to expect in terms of measuring the quality ofa decomposition. The selection of “a good set” of cutting planes from a large poolof available ones is a core topic in the computational cutting plane literature. Theory
suggests measures for the quality of single cutting planes, e.g., to prefer facet-defining
inequalities or deep(er) cuts etc. For selecting from a collection of cuts, however, atheoretical argumentation is much scarcer and less-founded, and appears only problemspecific in the literature.
In that light, a general theoretical a priori measure of what constitutes a good
Dantzig–Wolfe decomposition currently appears to be out of reach, and we may needto fall back to computational proxies like the one we propose.
Acknowledgments We sincerely thank an anonymous referee for thoughtful and motivating feedback,which led to a more meaningful experimental setup and a much improved organization of the material.

References
1. Achterberg, T., Koch, Th., Martin, A.: MIPLIB 2003. Oper. Res. Lett. 34(4), 361–372 (2006)
2. Aykanat, C., Pinar, A., Çatalyürek, Ü.V.: Permuting sparse rectangular matrices into block-diagonalform. SIAM J. Sci. Comput. 25, 1860–1879 (2004)
3. Bergner, M., Caprara, A., Furini, F., Lübbecke, M.E., Malaguti, E., Traversi, E.: Partial convexificationof general MIPs by Dantzig–Wolfe reformulation. In: Günlük, O., Woeginger, G.J., (eds.) Integer
Programming and Combinatorial Optimization, volume 6655 of Lect. Notes Comput. Sci., pp. 39–51,
Springer, Berlin (2011)
4. Borndörfer, R., Ferreira, C.E., Martin, A.: Decomposing matrices into blocks. SIAM J. Optim. 9(1),236–269 (1998)
5. Caprara, A., Furini, F., Malaguti, E.: Uncommon Dantzig–Wolfe reformulation for the temporal knapsack problem. INFORMS J. Comput. 25(3), 560–571 (2013)
6. Cleveland, W.S.: Robust locally weighted regression and smoothing scatterplots. J. Am. Stat. Assoc.
74(368), 829–836 (1979)
7. Colombo, M., Grothey, A., Hogg, J., Woodsend, K., Gondzio, J.: A structure-conveying modelling
language for mathematical and stochastic programming. Math. Program. Comput. 1, 223–247 (2009)
8. Desaulniers, G., Desrosiers, J., Solomon, M.M. (eds.): Column Generation. Springer, Berlin (2005)
9. Desrosiers, J., Lübbecke, M.E.: A primer in column generation. In: Desaulniers et al. [8], pp. 1–32.
10. Dolan, E.D., Moré, J.J.: Benchmarking optimization software with performance profiles. Math. Program. 91, 201–213 (2002)
11. du Merle, O., Villeneuve, D., Desrosiers, J., Hansen, P.: Stabilized column generation. Discret. Math.
194, 229–237 (1999)
12. Ferris, M.C., Horn, J.D.: Partitioning mathematical programs for parallel solution. Math. Program.
80(1), 35–61 (1998)
13. Fragnière, E., Gondzio, J., Sarkissian, R., Vial, J.-Ph.: A structure-exploiting tool in algebraic modelinglanguages. Manag. Sci. 46, 1145–1158 (2000)
14. Galati, M.V., Pratt, R.: The new decomposition algorithm in SAS/OR optimization. In: XXI International Symposium on Mathematical Programming (2012)123

M. Bergner et al.
15. Galati, M.V., Ralphs, T.K., Wang, J.: Computational experience with generic decomposition using the
DIP framework. In: Proceedings of RAMP 2012. COR@L Laboratory, Lehigh University (2012)
16. Gamrath, G., Lübbecke, M.E.: Experiments with a generic Dantzig–Wolfe decomposition for integer
programs. In: Festa, P. (ed.) Proceedings of the 9th Symposium on Experimental Algorithms (SEA),
volume 6049 of Lect. Notes Comput. Sci., pp. 239–252, Springer, Berlin (2010)
17. Geoffrion, A.M.: Lagrangean relaxation for integer programming. Math. Program. Stud. 2, 82–114(1974)
18. Karypis, G., Aggarwal, R., Kumar, V., Shekhar, S.: Multilevel hypergraph partitioning: application in
VLSI domain. IEEE Trans. VLSI Syst. 20(1), (1999)
19. Karypis, G., Kumar, V.: hmetis 1.5: a hypergraph partitioning package. Technical report, Departmentof Computer Science, University of Minnesota (1998)
20. Koch, Th., 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.:
MIPLIB 2010—mixed integer programming library version 5. Math. Program. Comput. 3(2), 103–163(2011)
21. Poggi de Aragão, M., Uchoa, E.: Integer program reformulation for robust branch-and-cut-and-pricealgorithms. In: Proceedings of the Conference Mathematical Program in Rio: A Conference in Honour
of Nelson Maculan, pp. 56–61 (2003)
22. Puchinger, J., Stuckey, P.J., Wallace, M.G., Brand, S.: Dantzig–Wolfe decomposition and branch-andprice solving in G12. Constraints 16(1), 77–99 (2011)
23. Ralphs, T.K., Galati, M.V.: DIP—decomposition for integer programming. https://projects.coin-or.
org/Dip (2009)
24. Ralphs, T.K., Galati, M.V., Wang, J.: DIP and DipPy: towards a decomposition-based MILP solver.
In: XXI International Symposium on Mathematical Programming (2012)
25. Tebboth, J., Daniel, R.: A tightly integrated modelling and optimisation library: a new framework for
rapid algorithm development. Ann. Oper. Res. 104(1–4), 313–333 (2001)
26. Tebboth, J.R.: A computational study of Dantzig–Wolfe decomposition. PhD thesis, University of
Buckingham (2001)
27. Vanderbeck, F.: BaPCod—a generic branch-and-price code. https://wiki.bordeaux.inria.fr/realopt/
pmwiki.php/Project/BaPCod (2005)
28. Vanderbeck, F., Wolsey, L.: Reformulation and decomposition of integer programs. In: Jünger, M.,
Liebling, Th.M., Naddef, D., Nemhauser, G.L., Pulleyblank, W.R., Reinelt, G., Rinaldi, G., Wolsey,
L.A. (eds.) 50 Years of Integer Programming 1958–2008. Springer, Berlin (2010)
29. Wang, J., Ralphs, T.K.: Computational experience with hypergraph-based methods for automatic
decomposition in integer programming. In: Gomes, C., Sellmann, M. (eds.) Integration of AI and
OR Techniques in Constraint Programming for Combinatorial Optimization Problems, volume 7874
of Lect. Notes Comput. Sci., pp. 394–402. Springer, Berlin (2013)
30. Weil, R.L., Kettler, P.C.: Rearranging matrices to block-angular form for decomposition (and other)algorithms. Manag. Sci. 18(1), 98–108 (1971)

123