Computers and Operations Research 84 (2017) 16–32

Contents lists available at ScienceDirect

Computers and Operations Research
journal homepage: www.elsevier.com/locate/cor

A strongly polynomial Contraction-Expansion algorithm for network
ﬂow problems
Jean Bertrand Gauthier a,∗, Jacques Desrosiers a, Marco E. Lübbecke bab

GERAD & HEC Montréal, 3000, chemin de la Côte-Sainte-Catherine, Montréal (Québec) Canada, H3T 2A7, Canada
RWTH Aachen University, Chair of Operations Research, Kackertstraße 7, D-52072 Aachen, Germanya r t i c l e

i n f oArticle history:
Received 23 March 2016
Revised 25 February 2017
Accepted 27 February 2017
Available online 1 March 2017
Keywords:
Network ﬂow problem
Residual network
Contracted network
Minimum mean cost cycle
Complexity analysis
Strongly polynomial algorithma b s t r a c t
This paper addresses the solution of the capacitated minimum cost ﬂow problem on a network containing n nodes and m arcs. Satisfying necessary and suﬃcient optimality conditions can be done on theresidual network although it can be quite time consuming as testiﬁed by the minimum mean cyclecanceling algorithm (MMCC). We introduce a contracted network which exploits these conditions on
a much smaller network. Since the construction of this contracted network is very ﬂexible, we study
its properties depending on the construction choice. A generic contraction algorithm is then produced
around the contracted network. Interestingly enough, it turns out it encapsulates both the MMCC andprimal network simplex algorithms as extreme cases. By guiding the solution using a particular expansion scheme, we are able to recuperate theoretical results from MMCC. As such, we obtain a strongly
polynomial Contraction-Expansion algorithm which runs in O(m3 n2 ) time. There is thus no improvementof the runtime complexity, yet the expansion scheme sticks to very practical observations of MMCC’s behavior, namely that of phases and jumps on the optimality parameter. The solution time is ultimately
signiﬁcantly reduced, even more so as the size of the instance increases.
© 2017 Elsevier Ltd. All rights reserved.

1. Introduction
The primal network simplex algorithm performs surprisingly
well considering that only a small fraction of iterations induce
non-degenerate pivots (see Ahuja et al., 1993, Fig. 18.7). In fact, thenetwork simplex and the cost scaling methods are typically preferred over competing alternatives for the solution of capacitatedminimum cost ﬂow problems (CMCF). Kovács (2015) suggests theuse of the former for smaller networks and the latter for larger
instances.
The minimum mean cycle-canceling algorithm (MMCC) is aprominent such alternative. Introduced by Goldberg and Tarjan(1989), this algorithm copes with degeneracy at the expenseof a more involved pricing problem able to identify improving
directions only, that is, with strictly positive step sizes. Despite
its strongly polynomial time complexity, the theoretical behaviorof MMCC is in practice no match for other methods. Radzik
and Goldberg (1994) even improve the complexity some ﬁve
years later and introduce the concept of phases. Two decades
further down the road brings another improvement due to∗

Corresponding author.
E-mail addresses: jean-bertrand.gauthier@hec.ca (J.B. Gauthier), jacques.
desrosiers@hec.ca (J. Desrosiers), marco.luebbecke@rwth-aachen.de (M.E. Lübbecke).
http://dx.doi.org/10.1016/j.cor.2017.02.019
0305-0548/© 2017 Elsevier Ltd. All rights reserved.

Gauthier et al. (2015) which combines phases with Cancel-andTighten (CT) presented alongside MMCC in the seminal paper as aself-standing algorithm.
This paper presents a Contraction-Expansion algorithm (CE)inspired by MMCC and its complexity proof. Since the visual aid
granted by the network ﬂow formulation gives a lot of perspective
to the theoretical analysis, we opt to present algorithmic choices ina constructive fashion. We adopt the deﬁnitions and nomenclatureof Ahuja et al. (1993).
The paper is organized as follows. Section 2 deﬁnes the network
problem and exposes the building block of this paper, namely theso-called contracted network. The contracted network is a ﬂexible
construction which gives rise to a generic contraction algorithm
whose properties are discussed in Section 3. A behavioral study ofthe so-called optimality parameter follows in Section 4. From theseobservations, an expansion scheme guiding the solution process
is drafted in Section 5, where the ensuing complexity analysis
expands upon theoretical results from the minimum mean cyclecanceling algorithm. Our ﬁnal thoughts can be found in Section 6.
We should also mention that the goal of this research is not tobeat leading algorithms on network problems. Our intent is a better understanding of algorithms with strictly positive step size at
every iteration. We believe that a more involved pricing problem
is worth investigating in methods where rather complex mathe-

J.B. Gauthier et al. / Computers and Operations Research 84 (2017) 16–32

17matical programs are already solved to guide the solution process,evoking in particular, the oracle of a column generation algorithm.
2. Network problem
Assume a capacitated directed network G := (N, A), where
N denotes the set of n nodes and A the set of m arcs. The arc
parametrization is captured by the cost vector c := [cij ](i, j) ∈ A
and non-negative bounds [ij , uij ], (i, j) ∈ A. A supply or demand,
respectively, deﬁned by a positive or negative value bi , i ∈ N, isassociated with each node such thati∈N bi = 0. Supported by
G and the vector of bounded ﬂow variables x := [xij ](i, j) ∈ A , aformulation of CMCF is given by:z∗ := min



ci j xi j(i, j )∈As.t.

xi j −

j:(i, j )∈A

i j ≤ xi j ≤ ui j ,

x ji = bi ,[ πi ] ,∀i ∈ N,j:( j,i )∈A∀(i, j ) ∈ A,

(1)

where π := [πi ]i∈N is the vector of dual variables, also known as
node potentials. We further assume that G does not contain any
multiarc, i.e., two or more arcs sharing both the same head andtail nodes.
When the right-hand side b := [bi ]i ∈ N is the vector of all zeros, we obtain a circulation problem. The latter is also the supportof the minimum mean cycle-canceling algorithm. Indeed, its corecomponent, namely the so-called residual network is a circulation
problem for which upper bounds are eventually neglected while
searching for an improved solution.
We recall the path and cycle deﬁnitions used by Ahuja et al.
(1993). A path from node i1 to node ir (abbreviated i1 ir ) ina directed graph G = (N, A ) is a sequence without repetition ofnodes and arcs i1 − a1 − i2 − a2 − . . . − ir−1 − ar−1 − ir , satisfying
the property that either ak = (ik , ik+1 ) ∈ A (forward arc) or ak =(ik+1 , ik ) ∈ A (backward arc) for all k, 1 ≤ k ≤ r − 1. The sequence
is typically given using nodes only. A directed path is an oriented
version of a path consisting of forward arcs only. A cycle is a path
together with either the forward arc (ir , i1 ) ∈ A or the backward
arc (i1 , ir ) ∈ A. A directed cycle is a directed path together withthe arc (ir , i1 ) ∈ A (which indeed imposes the orientation of thedirected cycle). The cost of a path or a cycle is the sum of the coston the forward arcs minus that on the backward arcs. Note that
there is no backward arc in a directed path nor a directed cycle.
With respect to any dual variable vector π := (πi )i∈N , thereduced cost of xij , (i, j) ∈ A, is deﬁned as c̄i j := ci j − πi + π j . The
reduced cost of a path or a cycle is then computed analogously
to its cost. Therefore, although Ahuja et al. (1993) state thefollowing fundamental property for a directed cycle, the same
telescopic summation argument proves to be enough regardless of
its directed nature.
Proposition 1. The cost and reduced cost of a cycle, directed or not,are equal.
In the following we deﬁne the traditional residual network
using which we can determine whether the current solution
x0 := [x0i j ](i, j )∈A can be improved or not. An optimality certiﬁcateis then provided which indeed amounts to the statement ofcycle-canceling algorithms. We ﬁnally move on to a contraction
manipulation which induces a so-called contracted network.
2.1. Residual network
The residual network allows a marginal construction of the
ﬂow that may traverse the network aside the current ﬂow x0(increasing ﬂow on certain arcs, possibly decreasing it on others).

Fig. 1. A change of variables.

The combination of x0 along with an optimal residual ﬂow computed on the residual network is optimal in the original network.
As depicted in Fig. 1, each arc (i, j) ∈ A can be replaced by two
residual arcs representing possible ﬂow increments or decrementsthat depend on the remaining capacities:• a forward arc (i, j) with cost dij := cij and residual ﬂow
0 ≤ yi j ≤ ri0j := ui j − x0i j ;• a backward arc (j, i) with cost d ji := −ci j and residual ﬂow
0 ≤ y ji ≤ r 0ji := x0i j − i j .
Denoted G(x0 ) := (N, A(x0 )), the residual network with respect
to x0 reﬂects the change of variables yi j − y ji := xi j − x0i j , (i, j) ∈ A.
It is based on the original nodes in N and the set of residual arcs
A(x0 ). Indeed, among the possible arc support A := {(i, j) ∪ (j, i)|(i,j) ∈ A} only those arcs with strictly positive residual capacities areof interest, i.e., A(x0 ) := {(i, j ) ∈ A | ri0j > 0}.
An equivalent formulation of (1) using the new y-variables
reads asz∗ = z0 + mins.t.



j:(i, j )∈A(x0 )0 ≤ yi j ≤ ri0j ,di j yi j(i, j )∈A(x0 )yi j −

y ji = 0,[ πi ] ,∀i ∈ N,j:( j,i )∈A(x0 )∀(i, j ) ∈ A(x0 ),(2)

where z0 := c x0 is the objective function value of the solution x0 .
Observe that using both forward and backward arcs between thesame pair of nodes can be simpliﬁed to sending the net ﬂow in asingle direction only. This is why we assume yi j y ji = 0, (i, j) ∈ A.
Fig. 2 exhibits the construction of the residual network G(x0 ).
Fig. 2a contains arc ﬂow variables x0i j , (i, j) ∈ A, at different values.
Each of these can be attributed a status depending on its actual
value: lower when x0i j = i j , upper when x0i j = ui j , or free wheni j < x0i j < ui j . When a variable is free, the ﬂow can be varied ineither direction, implying the presence of two anti-parallel residualarcs. However, when a variable is lower (resp. upper), this induces
only one arc oriented in the forward (resp. backward) direction.
An arc at its lower or upper bound is also said to be restricted.
Let the disjoint sets L(x0 ), U(x0 ), and F(x0 ) be a partition thatechoes the status of the original arcs at x0 . In order to simplify
the presentation, we denote this partition by {Lk , Uk , Fk }, where k≥ 0 refers to the solution xk in iteration k.
2.2. Optimality conditions
It is known that a ﬂow solution x∗ is optimal if and only if theresidual network G(x∗ ) contains no negative cost directed cycle,
shortened to negative cycle, (see Ahuja et al., 1993, Theorem 9.1).
By iteratively sending as much ﬂow as possible along negative
cycles, i.e., canceling negative cycles, and updating the residual
network accordingly, one obtains a generic cycle-canceling algorithm which terminates when there remains no negative cycle(Klein, 1967).
We recall here three equivalent necessary and suﬃcient
optimality conditions which can be interpreted in different perspectives with respect to network ﬂows (see Ahuja et al., 1993,

18

J.B. Gauthier et al. / Computers and Operations Research 84 (2017) 16–32

Fig. 2. Residual network construction where undirected edges represent two anti-parallel arcs.

Theorems 9.1, 9.3, and 9.4). The ﬁrst two refer to the residual
network whereas the third is based on the original network. In
regards to the dual condition perspective, with respect to any dual
variable vector π , the reduced cost of xij , (i, j) ∈ A, is deﬁned asc̄i j := ci j − πi + π j . Let the reduced cost d¯i j := di j − πi + π j of yij ,(i, j) ∈ A(x0 ), be computed in the same way.
Primal: G(x0 ) contains no negative cycle.
Dual: ∃π such that d¯i j ≥ 0, ∀(i, j ) ∈ A(x0 ).
Complementary slackness: ∃π such that, for every arc (i, j) ∈ A,
x0i j = i j if c̄i j > 0;x0i j = ui j if c̄i j < 0;c̄i j = 0 if i j < x0i j < ui j .

(3)
Although the notion of cycle cancellation can be seen as an
intuitive step size method, its foremost intention is the eliminationof said cycle from the current residual network. Indeed, passing astrictly positive ﬂow on any directed cycle W in G(x0 ) obviously
permits a transition between solutions, yet maintaining feasibility
is ensured by limiting the ﬂow to at most the smallest residual
capacity of the arcs forming the cycle, say ρ := min(i, j )∈W ri0j . A directed cycle is canceled when the step size is equal to ρ such thatat least one of the residual capacities ri0j , (i, j ) ∈ W, is saturated.
Note that any negative cycle qualiﬁes for an improvement withrespect to the objective function value since ρ > 0 by deﬁnition.
As such, various oracles capable of coming up with negative cycles
are conceivable, not all of these leading to strongly polynomialalgorithms (see Ahuja et al., 1993, Section 9.6).
Since negative cycles are canceled sequentially, the oracle
essentially deﬁnes a cycle-canceling algorithm. The minimum
mean cycle-canceling is one such algorithm famous for its strongly
polynomial runtime complexity, see Gauthier et al. (2015) for a
recent survey. MMCC uses an oracle which dismisses residual
capacities from the residual network and only identiﬁes negative
cycles of minimum mean cost. Technical aspects of this algorithm
are presented in Section 4.2 where we build upon these results
to show the runtime complexity of our Contraction-Expansion
algorithm. The contracted network permits the exploration ofalternative oracle constructions hoping to identify negative cycles

faster. Whether these cycles are directed or not on the residual
network remains to be seen.

2.3. Contracted network
A cycle free solution x0 does not contain any cycle of free arcs
(Ahuja et al., 1993). In primal simplex terminology, such and only
such solutions are basic. If one assumes a basic solution x0 , theset of free arcs F0 then deﬁnes a forest, that is, a collection of
node-disjoint trees.
We arbitrarily identify each tree with one of its nodes which
we call its root. For every node i ∈ N, let R(i ) be its associated root
(or tree identiﬁer) such that any two nodes i = j in N belonging
to the same tree must have the same root node R(i ) = R( j ).
Moreover, if the node i is a root, then R(i ) = i. For the example
in Fig. 3a, root nodes 1, 2, and 3 have been chosen such thatR(1 ) = R(4 ) = R(5 ) = R(6 ) = 1. Also, as a visual aid, we introduce the notion of a tree-layer to group the trees consisting offree arcs. Let G(F0 , x0 ) denote the residual network on x0 where
the tree-layer deﬁned with respect to the set of free arcs F0 is
superposed. Note that G(F0 , x0 ) ≡ G(x0 ). In fact, the arcs within theclouds in Fig. 3a are bidirectional such that one must imagine thetree-layer with respect to the original arcs in A. Based on the treelayer, we perform a contraction of every tree to its root node which
effectively hides all free arcs and lets the remaining restricted arcs
in L0 ∪ U0 be visible. The tail and head for each of these visible arcs
(i, j) is respectively redeﬁned to R(i ) and R( j ). This contraction
is likely to produce multiarcs. For instance, the arcs (2, 6), (7, 4),
and (8, 5) have the same head and tail in the contracted network.
Granted the actual cost computation has yet to be addressed, trivial cost dominance rules can be applied on such multiarcs. The end
result appears in Fig. 3b where the presence of the two dominatedarcs between root nodes 1 and 2 concerns only eﬃciency matters.
Recall that a tree T is a connected graph that contains no cycle.
In linear algebra terms, the coeﬃcient columns in (1) associatedwith the arcs of T are linearly independent. By extension, we calla subset of arcs linearly independent if no cycle can be formed.
Let us formalize the construction of the contracted network using
a general tree-layer deﬁnition. While a great deal of attention isJ.B. Gauthier et al. / Computers and Operations Research 84 (2017) 16–32

19

Fig. 3. Contracted network based on the set of free arcs, arcs within a cloud are contracted and become hidden while the others remain visible.

given to the set of free arcs F0 , the truth is that the tree-layer canbe deﬁned with respect to an arbitrary although linearly independent set of arcs, say P0 ⊂ A, where the superscript is omitted unless
iterate comparisons are required. The set P is then used to partition the arcs of A(x0 ) in two categories. As such, one can think of
the sets HP (x0 ) and VP (x0 ) as those arcs that are hidden and visible,
respectively, in the contracted network. These sets are formed by

HP (x ) :=0

(i, j )∈P

(i, j ), ( j, i ),(i, j ),( j, i ),if(i, j ) ∈ F 0
if(i, j ) ∈ L0
if(i, j ) ∈ U 0

VP (x0 ) := A(x0 ) \ HP (x0 ).

(4)

(5)

Deﬁnition 1. With respect to the set P, the tree-layer identiﬁes
root nodes in the set

NP (x0 ) := {i ∈ N | R(i ) = i},(6)

as well as the arc partition HP (x0 ) and VP (x0 ) of the residual
network G(x0 ).
The contracted network H(P, x0 ) is then obtained from the setsof root nodes and visible arcs:

H (P, x0 ) := (NP (x0 ), VP (x0 )).

(7)VP (x0 )where every visible arc (i, j) in the set
is remapped to(R(i ), R( j )) thus maintaining a bijection between arcs of thecontracted network and the residual arcs.
We present two alternative contraction examples which collectively cover all possibilities. The ﬁrst example uses a subset offree arcs P ⊂ F0 while the second combines free arcs with some restricted arcs such that P⊃F0 . In the most general case, one can mix
both possibilities by using some elements of F0 and some of AF0 .
Fig. 4 covers the ﬁrst example where the subset of free arcs
used is P = F 0 \ {(5, 6 )}. In Fig. 4a, the free arc (5, 6) has been
duplicated in both directions and consequently appears in thevisible set VP (x0 ) with both arcs now being at a lower bound ofzero in the residual network. We refer to this kind of manipulation

on free arcs as coerced degeneracy. This yields a larger contracted
network since there is now a forest with four trees to handle.
Fig. 4b portrays the consequent contracted network.
As additional free arcs are duplicated in both directions, one
eventually reaches a point where the set P⊆F0 is empty, such thatno free arcs are hidden at all, hence yielding a contracted network
whose arc set is the same as that of the residual network, i.e.,V∅ (x0 ) = A(x0 ). The reader is now invited to consider the other extreme case where the tree-layer consists of a single spanning tree.
Such is the content of Fig. 5 where the set P = F 0 ∪ {(7, 4 ), (3, 9 )}is still linearly independent and consists of the union of the freearcs along with two additional restricted arcs. Without loss of generality, assume these two arcs are basic degenerate in the primal
simplex sense such that P = B0 ≡ B(x0 ) corresponds to a set of basic arcs at x0 . The tree-layer G(B0 , x0 ) seen in Fig. 5a then splits
the arcs in two subsets: The nine basic arcs of the spanning tree
and the seven nonbasic arcs. The contracted network H(B0 , x0 ) appears in Fig. 5b, where the spanning tree is contracted to the soleroot node hiding all basic arcs in the process and leaving visible
the nonbasic ones. Each nonbasic arc becomes a loop, indeed, adirected cycle on H(B0 , x0 ).
The arc cost computations of the contracted network are coming in the midst of the following section where features and properties of the contracted network are derived according to the
choice of the set P, most notably the nature of the optimality conditions the oracle derived from the contracted network is able to
fulﬁll. Since the discussion revolves around the content of the set
P, it is worthwhile to underscore that only two cases are possible
for an arc contained in the latter: It is either free or it is not.
3. Contracted network properties
Although the selection of the set P is fairly arbitrary,
Section 3.1 addresses how easy it is to meet the linear independence requirement regardless of the current solution x0 being
basic or not. Section 3.2 then states that any contracted cycle, adirected cycle identiﬁed on the contracted network, is uniquely

20

J.B. Gauthier et al. / Computers and Operations Research 84 (2017) 16–32

Fig. 4. Contracted network with coerced degeneracy on the free arc (5, 6).

Fig. 5. Contracted network with a set of basic arcs (primal network simplex algorithm).

extended on the residual network G(x0 ). The analysis of thecontracted network’s arc cost is performed in Section 3.3 which
is followed with the so-called pricing problem and the algorithm
in Section 3.4. We then derive in Section 3.5 the optimality certiﬁcate whose nature depends on the selected set P. As the set P
inﬂuences the content of the contracted network, different known
algorithms are referenced in Section 3.6 by examining the possible
outcomes of the pricing problem. Finally, Section 3.7 shows thatonce the set P is selected, the remaining arbitrary decisions thatmust be made have no impact on the algorithm.
3.1. Nonbasic solution
The set F0 trivially fulﬁlls the linear independence assumption
when the current solution x0 is basic. We show that the basic status is not restrictive for the selection of the set P. Fig. 6a presents a
nonbasic solution where a cycle of free arcs is contained in a given
feasible solution. Figs. 6b and 6c handle the issue using two alternative mechanisms. With the ﬁrst mechanism, the cycle of freearcs is canceled, in either direction, yielding at least one restricted

arc within the cycle. Since the cancellation ultimately modiﬁes
the current solution, one may altogether prefer the improving
direction, say saturating the arc (5, 6) to its upper bound. With thesecond mechanism, the arc (4, 1) is coerced degenerate in G. Thisduplicating manipulation does not change the current solution
yet provides a fast way to eliminate the cycles of free arcs thus
allowing one to deﬁne the set P using only independent arcs of F0 .
Letting f := |F0 | ≤ m, Proposition 2 asserts that, regardless ofwhether x0 is basic or not, it is easy to maintain a tree-layer induced by a linear independent set P.
Proposition 2. A linear independent set P ⊂ A(x0 ) can be derived
from F0 either by removing any cycle of free arcs using at most f cycle
cancellations or by applying coerced degeneracy to at most f arcs.
Proof. It is trivial to verify that rendering at most f arcs coerced
degenerate means that there remains a suitable linearly independent subset of free variables capable of forming a tree-layer.
In fact, one may think of these coerced degenerate variables as
super-basic such that canceling a cycle of free arcs amounts to

J.B. Gauthier et al. / Computers and Operations Research 84 (2017) 16–32

21

Fig. 6. Working with a nonbasic feasible solution on network G.

Fig. 7. Rooted path R(i )  R( j ) associated with visible arc (i, j).

the recovery of an optimal [or not] basis which can be done in astraightforward manner according to Marsten et al. (1989). 
3.2. Uniqueness of the extended cycle
Observe that in the tree-layered residual network, any visiblearc (i, j) ∈ VP (x) connects two root nodes using ﬁrst a path fromR(i ) to i followed by arc (i, j) and second a path from j to R( j ).
Fig. 7 portrays such an alternated path-arc-path sequence, fromnow on called a rooted path. When performing the contraction, it
is convenient to let P (i ), i ∈ N, denote the path from node i toroot node R(i ). The rooted path associated with the visible arc(i, j) ∈ VP (x) can therefore be decomposed as {−P (i ), (i, j ), P ( j )},where −P (i ) is the reversed path P (i ) and the orientation of afree arc is the one given by (i, j).
Several remarks about the rooted path are noteworthy. The tworoot nodes can possibly be the same, see the arc compatibility
Deﬁnition 2. The two distinct paths are formed using only arcs
in HP (x0 ), that is, hidden arcs in tree structures. Either or both ofthese paths could be of length zero such that the deﬁning arc (i, j)∈ VP (x0 ) could support the rooted path by itself. Finally, by deﬁnition of a tree, any visible arc (i, j) ∈ VP (x0 ) induces a unique rooted
path.
Proposition 3. Any contracted cycle WH obtained on the network
H(P, x0 ) yields a unique extended cycle WH:G on the residual network
G(x0 ).
Proof. Since only arcs that are part of the bijection are used to
produce contracted cycles, the uniqueness of rooted paths guarantees that the extended cycle uniquely exists in the tree-layered
residual network. Once again, the orientation of the free arcs inthe extended cycle is given by that of the contracted cycle. 
Implementation. The remainder of this section is dedicated todescribing a contracted cycle obtained on the contracted network
H(P, x0 ) in more details. Let

WH := {(R(i1 ), R( j1 ) ), (R(i2 ), R( j2 )), . . . , (R(i|WH | ), R( j|WH | ))}be one such directed cycle on the contracted network. Observethat it produces a sequence of entry and exit nodes in the treelayer which ultimately cycle through the same root nodes, i.e.,

Rs := R( js ) = R(is+1 ), s ∈ {1, . . . , |WH |}, where i|WH |+1 abusively
equals to i1 . The extended cycle, seen on the residual network
G(x0 ), can then be expressed as a concatenation of rooted paths.
This may not immediately lead to an elementary cycle.
One may take a look at Fig. 8 should the following explanation
require visual support. Think of the path contained in a given
tree, say between the arcs (is , js ) and (is+1 , js+1 ), s ∈ {1, . . . , |WH |}.
Notice that both paths P ( js ) and −P (is+1 ) have at least one nodein common, namely the root node Rs . A cycle is formed during theconcatenation if and only if there exists some other common node
us = Rs . Consider for instance that P ( js ) follows the path js  us
and then us  Rs while −P (is+1 ) travels on Rs  us followed byus  is+1 . Since all nodes connecting us to Rs and vice versa arecommon to both paths, an elementary path is found by detecting
the node us ∈ P ( js ) closest to js thus eliminating any back andforth play across opposing arcs. At last, the extended cycle WH:G
extracted from the contracted cycle WH is given by

WH:G :=|WH |{(is , js ), P ( js ) \ P (us ), −P (is+1 ) \ −P (us )},

(8)

s=1

where some of the composing paths P ( js ) ∪ −P (is+1 ), ∀s ∈{1, . . . , |WH |}, may be truncated to obtain an elementary cycle.

3.3. Arc cost transfer policy
So far, it has been established that any contracted cycle WH
exists uniquely as WH: G on the residual network. In order to extend optimality conditions to the contracted network, we requirea somewhat opposite feature: If the residual network contains anegative cycle then so must the contracted network.
Since we assume that P contains only free arcs, optimality conditions imply that the reduced cost on all these arcs must be zeroby the complementary slackness conditions (3). This is in line withthe way dual variables are determined in the primal network simplex algorithm, that is, Ahuja et al. (1993, Chapter 11.4) arbitrarily
set to zero the dual variable of the root node allowing a uniqueestablishment for the remaining dual variables. It is a consequenceof this procedure that π i , i ∈ N, gives the cost of the path from i to
the root node (with the cost from the root to itself of zero). This inturn produces a reduced cost of zero for every arc in the spanning
tree. We reproduce the same scheme using the more generalarbitrary constant initialization of the root node in Procedure 1.
Procedure 1 (Compute dual variable values). The dual variableof a root node π i , i ∈ NP (x0 ), is ﬁxed to an arbitrary value. The
remaining dual variable values π i , i ∈ NNP (x0 ), are then the costof the path P (i ) plus πR(i ) .

22

J.B. Gauthier et al. / Computers and Operations Research 84 (2017) 16–32

Fig. 8. Elementary tree path detection between js and is+1 .

By construction, the reduced cost of every arc in the forest
induced by P is zero, i.e., d¯i j = 0, (i, j) ∈ HP (x0 ), indeed for all
those arcs that are hidden in the contracted network.
Proposition 4. Given π satisfying Procedure 1, the reduced costof the contracted cycle WH is equal to the cost of its extended
counterpart WH: G .
Proof. Hiding the zero-cost arcs of WH: G makes the reduced
cost of the contracted cycle WH equal to the reduced cost of its
extended version. By Proposition 1, as the reduced cost and thecost of a cycle (directed or not) are the same, theresult follows. Remark. Since the cost and reduced cost of a cycle are equal, theway dual variables are established in order to satisfy optimality
conditions effectively transfers cost information of hidden arcs in
HP (x0 ) to those that remain visible in the pricing problem VP (x0 ).
Indeed, as the reduced cost of the former set is zero, the reduced
cost of the latter set must essentially capture the original cost
information. Furthermore, in case of multiarcs in the contracted
network, it is possible to keep only a single arc with the smallest
reduced cost. Such an arc is non-dominated with respect to its
reduced cost.
3.4. Contraction algorithm(i, j )∈VP (x0 )s.t.

(i, j )∈VP (x0 )|R(i )=

yi j = 1,

yi j −

yi j = 0,[ π ] ,∀ ∈ NP (x0 ),(i, j )∈VP (x0 )|R( j )=

[ μ] ,

ρ k := mink rikj ≥ 0,(i, j )∈WH:G(10)kwhich is zero only if an arc of the extended cycle WH:
has a residG
ual capacity of zero. We then obtain the solution xk+1 where thek , that is,
ﬂow update, if any, is only performed on the arcs of WH:G⎧ kk⎨xi j + ρ k , ∀(i, j ) ∈ A | (i, j ) ∈ WH:
G
k+1
k
k
k
x
−
ρ
,
∀
(
i,
j
)
∈
A
|
(
j,
i
)
∈
W
xi j :=
H:G
ij
⎩xk ,
otherwise(11)

zk+1 := zk + ρ k |WHk | μkH .
Observe that when ρ k = 0, the set P k+1 must be different from
Pk for otherwise the algorithm would not terminate. Degeneracy
and cycling phenomena from the primal simplex algorithm cometo mind. Indeed, unless rules for zero step size cycles are included,
the same minimum mean cycle is identiﬁed in the next iteration. A
pseudo-code is elaborated in Fig. 9 where this process is repeateduntil optimality is reached.
The nature of the optimality certiﬁcate in the generic contraction algorithm goes hand in hand with the non-degeneracy guarantee provided (or not) by the selected set P and the induced contracted network. Let us analyze this property.
3.5. Nature of the optimality conditions(i, j )∈VP (x0 )yi j ≥ 0,

By Proposition 4, if the residual network contains a negative
cycle, we identify one by solving the pricing problem (9). By
construction, the current solution xk , k ≥ 0, is optimal whenμkH ≥ 0, that is, when the contracted network does not containany negative cycle. Otherwise, given the identiﬁed contracted cycle
WHk , one computes the non-negative step sizeij

We are now ready to present a contraction algorithm. We aim
to detect a negative cycle using the contracted network H (P, x0 ) =(NP (x0 ), VP (x0 )). Like in the pricing oracle of the minimum mean
cycle-canceling algorithm, residual capacities are omitted and webound the ﬂow values by a normalizing constraint (Dantzig andThapa, 2003, Section 10.2). Proposition 5 states the kind of cycles
that are identiﬁed using the following pricing problem:d¯i j yi j
μH := minarcs in W form either a cycle or a collection of cycles. Solutions of
the latter category cannot be expressed as extreme points of (9).
In the former case, the normalizing constraint forces the ﬂow oneach arc of the single cycle to be equal to 1/|W|. Finally, cost coefﬁcients of the pricing problem (9) are d̄ = [di j − πi + π j ](i, j )∈VP (x0 ) .
Factoring out 1/|W| and using Proposition 1, we arrive at
μH = |W1 | (i, j )∈W d¯i j = |W1 | (i, j )∈W di j . ∀(i, j ) ∈ VP (x ),0(9)where dual variables appear on the right between brackets. Flowconservation constraints are deﬁned for each root node where thesummation indices simply specify all arcs from VP (x0 ) that are
related to the speciﬁed root nodes.
Proposition 5. An extreme point solution to the pricing problem(9) corresponds to a minimum mean cost directed cycle on H(P, x0 ),averaged over the number of arcs it contains.
Proof. Deﬁne the strictly positive arc support set W := {(i, j ) ∈VP (x0 ) | y∗i j > 0} of an optimal solution. By ﬂow conservation, theGiven that any contracted cycle is uniquely extended, theunderlying expectation is that it is possible to travel on the hiddenarcs, i.e., the extended cycle is directed on G(x0 ). Depending on
the actual status of the arcs in set P, this expectation could bechallenged. Verifying that any set P which consists of only free arcs
always provides a contracted cycle WH on which at least one ﬂow
unit can always be sent on the extended cycle WH: G is trivial: Allarcs unaccounted for can be used in either direction. As supporting
evidence, the ﬁrst three extended cycle extractions illustrated inFig. 10 are directed in G(x0 ). They come from contracted cycles
obtained on a contracted network using only free arcs (Fig. 3b orJ.B. Gauthier et al. / Computers and Operations Research 84 (2017) 16–32

23

Fig. 9. Generic contraction algorithm for network ﬂow problems.

Fig. 10. Extended cycle extractions.

4 b). The fourth extended cycle extraction comes from Fig. 5b andyields an undirected cycle inducing a zero step size.
Fig. 10a is obtained from a contracted cycle deﬁned by the arc
loop on root node 1 in the contracted network of Fig. 3b. Fig. 10bcomes from the same contracted network but the contracted cycle
uses more arcs: {(1, 3), (3, 9), (8, 5)}. Fig. 10c identiﬁes part of theprevious cycle although it is based on the contracted network of
Fig. 4b where more visible arcs are available. Finally, Fig. 10d isbased on Fig. 5b used with the primal network simplex algorithm.
The arc (7, 4) is basic degenerate and one of the nonbasic arcs
corresponds to the visible arc (2, 6), a loop on root node 1. The
corresponding extended cycle is deﬁned by {(2, 6), (6, 5), (5, 4),(4, 7), (7, 2)}, or in accordance with the threefold decompositionof WH: G given in (8) as{(2, 6 ), {{(6, 1 )} \ {(6, 1 )}},{{(1, 6 ), (6, 5 ), (5, 4 ), (4, 7 ), (7, 2 )} \ {(1, 6 )}}}.
Unfortunately, the arc (4, 7) is not a residual arc and reduces
the possible ﬂow on the identiﬁed extended cycle to zero. In the

primal simplex algorithm, pivoting on arc variable x26 induces adegenerate pivot.
When P⊆F0 and μ0H < 0, the contracted network H(P, x0 ) is
guaranteed to identify a directed extended cycle, that is, when thehidden arcs are free, any contracted negative cycle has a strictly
positive step size on G(x0 ). This refers to the primal optimality
conditions of Section 2.2. On the other hand, when PF0 , there isno such guarantee since the extended cycle could contain an arc
with a residual capacity of zero such that the associated step size
is non-negative. In that case, μ0H ≥ 0 is only a suﬃcient optimality
condition. This proves the following proposition.
Proposition 6. Given the set P ⊂ A of linearly independent arcs, if• P⊆F0 , the oracle (9) provides necessary and suﬃcient optimality
conditions for (1), i.e., x0 is optimal if and only if μ0H ≥ 0.
• PF0 , the oracle (9) provides suﬃcient optimality conditions for(1), i.e., x0 is optimal if μ0H ≥ 0.
Regardless of the choice of the set P, arcs (i, j) ∈ A can be categorized in two classes: Those arcs that link nodes of the same tree,

24

J.B. Gauthier et al. / Computers and Operations Research 84 (2017) 16–32

i.e., the head and tail of the residual arc have the same root nodeR(i ) = R( j ), and those that link different trees, i.e., R(i ) = R( j ).
Deﬁnition 2. An arc (i, j) ∈ A is compatible with P if and only if
the head and tail refer to the same root node in the contracted
network, i.e., R(i ) = R( j ). Otherwise, it is incompatible.
Compatible arcs are easy to identify as loops on the contracted
network thus forming directed cycles on their own. Examples
can be seen in Fig. 3a as arcs (4, 1) and (9, 10) or in Fig. 3b as
the two loops on root nodes 1 and 2. This also implies that theextended cycle associated with a compatible arc (i, j) ∈ VP (x0 ) is
not inﬂuenced by the content of other trees. The reduced cost of
a loop is therefore the reduced cost of the cycle. For incompatible
arcs, rooted paths have to be combined to form cycles in H(P, x0 ).
In this respect, we like to think of the reduced cost d¯i j , (i, j) ∈VP (x0 ), as a rooted cost.
3.6. Extremal point solution space
It is straightforward to verify that the extreme cases, P = ∅and P = B0 , respectively, correspond to the minimum mean cyclecanceling algorithm and the primal network simplex algorithm.
In the former case, the contracted network H(∅, x0 ) has exactly
the same structure as the residual network G(x0 ), thus yielding a
pricing problem equivalent to MMCC’s. In the latter case, consider
the contracted network of Fig. 5 where the arcs of P are basic arcs
(indeed forming the primal network simplex spanning tree) andall visible arcs are nonbasic. Also recall that the arc cost transfer
policy seen in Procedure 1 is a generalized version of the oneused in the primal network simplex algorithm where only one root
node is used. As each visible arc is compatible (appearing as a loopin H(B0 , x0 )), applying cost dominance trivially results in a single
contender for the pricing problem to identify which incidentallyreconciles us with the entering variable of Dantzig’s pivot rule.
Let CP denote the set of all directed cycles obtainable fromextreme points of the pricing problem (9). The cardinality of thisset can be measured in two parts, that is, the compatible and incompatible portions. The ﬁrst portion is insigniﬁcant and reduces
to a single element regardless of the set P while the second grows
exponentially as the set P reduces in cardinality.
Proposition 7. The cardinality of the set of extreme point solutions|CP | of the pricing problem (9) on the contracted network H(P, x0 ) is2
O(2(n−|P|+1 ) ).
Proof. For the sake of simplicity, the argument is carried witha dominance rule applied on the visible arcs. Furthermore, we
use a worst case type analysis on the density of these sets. As
an extreme point of (9) can only identify a single compatible arc,
dominance can be applied across all compatible arcs regardless of
their root node association.
The more interesting portion thus concerns incompatible
combinations which of course require at least two incompatible
arcs. Using combinatorial enumeration, let us count the cycles inrough numbers by randomly selecting arcs within the available
possibilities VP (x0 ) thus forming directed cycles of different sizes
ranging from 2 to n − |P | + 1, the number of nodes in the conn−|P|+1 |VP (x0 )|tracted network, i.e., |CP | = k=2. Then again, oncekdominance is applied, the set of visible arcs vastly overestimates
the actual number of arcs remaining in the contracted network
which for a complete graph amounts to (n − |P | + 1 )2 . Since basicnn22
calculus reduces
to the dominant term 2n , we obtaink
k=1|CP | ∈ O(2(n−|P|+1) ). 
2

Granted it is possible to order the cardinalities, |C∅ |  |CF | |CB |, the same cannot be done for the actual sets. Indeed, hidden
arcs and dominance rules means that directed cycles present in asmaller set are not necessarily present in a larger one. Consider for
instance Fig. 10b and c where the cycle uses the non-dominatedarc (8, 5) for the former and (7, 4) for the latter.
While this computation tremendously dramatizes the size ofthe extreme point solution set of the pricing problem (9), it gives
the general intuition that some balance can be achieved between
the workload offset transferred to the oracle and the simplicity ofobtaining undirected cycles that may induce a step size of zero.
Let us put this in perspective of the contracted networks seen in
Figs. 3b, 4b, and 5b. The formulation of the pricing problem (9) isalways the same, yet the contracted network reminiscing of the
primal network simplex algorithm (5b) has |CB | = 1 whereas many
more cycles are present in the contracted network using all freevariables (3b). While the density of the resulting contracted network converges towards the density of the original network suchthat the exponential bound becomes increasingly approximative
as |P| get smaller, it is nevertheless expected that the diﬃculty ofsolving the pricing problem (9) follows the burden of the extreme
point solution set.
Furthermore, for the exponential cardinality growth to be inline with Proposition 6, one must realize that additional cycles
available through a smaller set P must come from the combinationof cycles available in a larger set P. Consequently, since cycles
available using the set P = F are suﬃcient to ascertain optimality,
some of the cycle combinations available in the pricing problem
using set P ⊂ F must lead to nonbasic solutions. In this matter,the purpose of Section 3.1 which handles nonbasic solutions is
twofold. Sure enough, it manages the initial starting solution
but more importantly solutions encountered during the solution
process. That is, in particular there is no guarantee that MMCC
(with P = ∅) travels through extreme points and the same is true
for other cases where actual free arcs are rendered visible in VP ,that is, when P ⊂ F.

3.7. Root selection
Since the selection of root nodes is arbitrary, the reader might
wonder what impact, if any, a different set of root nodes would
have on the oracle and thus the algorithm’s course. As shownin the proof of Proposition 8, it turns out very little. Indeed, theoracle’s content is modiﬁed on a per component basis yet as awhole it is completely unaffected.
Proposition 8. The root selection has no inﬂuence on the pricing
problem (9). The compatible set, the extreme point solution set CP , themean cost evaluation μH and the nature of the optimality certiﬁcatein Proposition 6 are unaltered.
Proof. On the one hand, it is no surprise that, for all nodes i ∈ N,
the paths P (i ) and their costs π i are modiﬁed when an alternative
set of root nodes is used. This means that while each visible arc(i, j) of VP (x0 ) continues to exist on H(P, x0 ), it is now associatedwith a different reduced or rooted cost d¯i j = di j − πi + π j . The
ﬁrst observation is that the set of extreme point solutions CP
corresponding to directed cycles obtained by linear combinationof visible arcs is unaltered. Furthermore, by Proposition 4 theoriginal cost on each cycle is maintained regardless of the root
selection. As such, the mean cost of every extended cycle is also
maintained which obviously means as much for the nature of
the optimality certiﬁcate. Finally, observe that not only is thecompatibility status of any arc (i, j) persistent, the reduced orrooted cost of a compatible arc is also immune to change. 

J.B. Gauthier et al. / Computers and Operations Research 84 (2017) 16–32

This is a testament to the fact that a degenerate solution
grants degrees of freedom for dual variables, in fact one for eachdegenerate constraint in accordance with Procedure 1. Intuitively
speaking, the pricing problem (9) can be interpreted to optimize
over these degrees of freedom. Indeed, solving the former not
only provides an optimal cycle WH0 of minimum mean reduced
cost μ0H on the contracted network but also a new set of dual
variable values for the root nodes. The following is an adaptationof Goldberg and Tarjan (1989, Theorem 3.3) or, more directly interms of notation, of Gauthier et al. (2015, Proposition 3) which
also bases its argument on linear programming theorems.
Proposition 9. When solving the pricing problem (9) at iteration k≥ 0, there exists some dual variable values πk ,  ∈ NP (xk ), such thatkkkkd¯i j − πR(i ) + πR( j ) ≥ μH , (i, j) ∈ VP (x ). By complementary slackness,the latter in fact hold at equality for all strictly positive variables,i.e., (i, j ) ∈ WHk .
In other words, the initial dual variable values on root nodes
are irrelevant as much as the roots selection. Either way, thepricing problem (9) optimizes dual variables aiming to satisfy
optimality conditions.
4. Behavioral study
Supported by numerical results from a minimum cost ﬂow
problem instance containing 1 025 nodes and 92 550 arcs (referred to as Instance 1), in Section 4.1 we analyze the behavior ofa speciﬁc variant of the contraction algorithm opposite MMCC’s
established behavior. This study not only serves to grasp somemechanical aspects of these algorithms, but more importantly
draws attention on key points of the theoretical minimum mean
cycle-canceling complexity. Section 4.2 connects the dots with
Cancel-and-Tighten and digs into more advanced aspects.
4.1. A lower bound on the minimum mean cost
Our analysis opposes the mean cost of negative cycles obtainedon the contracted network to those that would be obtained in theframework of MMCC. At a given iteration k ≥ 0, let μkG be theminimum mean cost of the directed cycle WGk on the residual network as obtained by MMCC. We then denote by μkH the minimum
mean cost of the directed cycle WHk fetched on the contracted
network andμkH:G = μkH|WHk |k|WH:|Gmight be better understood by observing the evolution of thevarious μ during the solution process. For illustrative purposes, we
systematically use the contracted network H(Fk , xk ) at every iteration, where all and only the free arcs are hidden. The value of μkH
and its counterpart μkH:G are naturally obtained on top of which
we also poll for the value μkG as if to look for the minimum mean
cycle. Fig. 11 shows the evolution of these measures for Instance 1.
While all three plots show a general increasing tendency,
what is striking is how different the inequality μkH ≤ μkG is fromμkG ≤ μkH:G . The ﬁrst gap is fairly intuitive and goes back to basic
mathematics: The density of the contracted network produces
contracted cycles which use relatively few arcs compared to theirextended counterparts. As the cycle costs remain the same, thedenominator strongly inﬂuences the mean evaluation. The ordering (13) can then be interpreted as a deceptively small minimum
mean cost μkH whose value is corrected when accounting for thehidden arcs. Observe that the ordering (13) is equal throughout ifk |.
there are no hidden arcs in the contracted cycle, i.e., |WHk | = |WH:
G
The second gap is much more tight and deserves more attention. In this matter, Fig. 12 zooms on the evolution of theextended mean cycle cost μkH:G and minimum mean cycle cost μkG .
For the record, MMCC both searches for and applies a minimum
mean cycle at each iteration, yielding an algorithm which featuresa non-decreasing property on μkG , see Goldberg and Tarjan (1989,
Lemma 3.5) or Gauthier et al. (2015, Proposition 4). Take a close
look around the 10 0 0 0th iteration. This is not a graphical aberration showing that this property is indeed lost when negative
cycles are not canceled in the order suggested by MMCC. What
we think is surprising for this particular instance is how little thisphenomenon occurs, only 21 times within 14 258 iterations to be
exact.
Convergence. As the contraction algorithm identiﬁes a negative
cycle at every iteration, it is evident that convergence of the objective function to optimality is guaranteed. By default, there arethen at most O(mCU) such negative cycles, a weakly polynomial
result referring to the largest absolute cost (C := max (i, j) ∈ A |cij |)and interval bound values (U := max(i, j )∈A ui j − i j ).
Nevertheless, comparing the canceled cycle mean cost to that of
the minimum mean cycle cost is an enlightening exercise. Indeed,
the erratic behavior also appears in a strongly polynomial algorithm, namely Cancel-and-Tighten. Section 4.2 recalls some concepts from the minimum mean cycle-canceling algorithm (Gauthieret al., 2015) for which strongly polynomial runtime is establishedthrough the analysis of the so-called optimality parameter.

(12)

the mean cost of the extended cycle where hidden arcs are accounted for. The following proposition is veriﬁed by construction.
Proposition 10. Let xk , k ≥ 0, be a non-optimal solution. Given theset P⊆Fk , the following ordering of the minimum mean costs always
holds:

μkH ≤ μkG ≤ μkH:G .

25

(13)

k
Proof. The cycle WH:
is visible as is on the residual networkGk
G(x ) which means μkG ≤ μkH:G trivially holds by the nature of minimizing the mean cost in G(xk ). Furthermore, although contracted,

any cycle WGk that also appears in H(P, xk ) is evaluated with adifferent mean cycle cost which eventually uses less arcs such thatμkH ≤ μkG . 
Since we aim to establish complexity results in light of MMCC’s
analysis, it is quite natural to think of WGk as the reference minimum mean cost cycle. The inequalities literally state that the
optimal value of the minimum mean cost μkG may be underestimated by μkH and overestimated by μkH:G . The notion of estimation

4.2. Optimality parameter analysis
Strongly polynomial runtime complexity certiﬁes that an algorithm performs a polynomial number of operations measured
solely by the size of its input (see Cormen et al., 2009). For network ﬂow problems, this number should be a function of m and n
only. In introducing the minimum mean cycle-canceling algorithm,
Goldberg and Tarjan (1989) also provide such a complexity proof
which holds in two parts. First, an oracle capable of producing
a minimum mean cost directed cycle in O(mn) time. Second, an
optimal solution is reached by canceling at most O(m2 nlog n) suchnegative cycles. The ﬁrst part can be seen as the inner loop andthe second as the outer loop. Their product then yields a strongly
polynomial global complexity. Radzik and Goldberg (1994) reﬁnethe complexity analysis, reducing the number of cycle cancellations to O(m2 n). Although they also introduce the concept ofphases to analyze the behavior of the algorithm, Gauthier et al.
(2015) strongly insist on the latter to further improve the complexity result by combining the Cancel-and-Tighten strategy introduced
in Goldberg and Tarjan (1989) with the original algorithm.

26

J.B. Gauthier et al. / Computers and Operations Research 84 (2017) 16–32

Fig. 11. Comparison of μkH , μkG and μkH:G .

Fig. 12. Zoom on comparison of μkG and μkH:G .

In MMCC, the value μkG is coined as the optimality parameterbecause it converges without oscillations to 0 from below instrongly polynomial time. It might however be more appropriate to reserve this name for another value as illustrated by theupcoming analysis.
Type 2 (negative) cycles. Even though the non-decreasing
property of μkG across iterations in the minimum mean cyclecanceling algorithm is interesting, it has been argued by Gauthieret al. (2015) that the strictly increasing behavior observed acrossphases is more enlightening. The phase deﬁnition goes hand inhand with the proof of Gauthier et al. (2015, Proposition 5) which
distinguishes between Type 1 and Type 2 cycles. It shows that a
Type 2 cycle is attained within m cancellations or optimality isachieved, where a Type 2 cycle is obtained when there exists atleast one variable in the cycle with a non-negative reduced costcomputed with respect to a set of dual variables established at the
beginning of a phase. Given μhG at the beginning of phase h ≥ 0,such a Type 2 cycle W then serves to imply a strict increase on
μhG asμG ≥ 1 − 1/|W  | μhG ≥ (1 − 1/n )μhG

(14)

thus marking the end of the phase, i.e., the sequence of iterations
leading to this measurable jump factor.
Since each phase h ≥ 0 contains at most m cancellations, thenumber of cancellations can be interpreted as at most O(mnlog n)or O(mn) phases depending on the complexity point of view
(Gauthier et al., 2015). These points of view respectively refer to
Theorem 2 which uses the minimal increasing factor (1 − 1/n )and Theorem 3 rather exploiting the stronger factor (1 − 1/|W  | ).
While the concept of phases is useful for the complexity analysisof MMCC, it is not transparent at all in the implementation.
Indeed, in MMCC, the phase is purely a question of theoretical
existence; dual variables are never required to advance suchthat the solution process cares not about these cycle Types. The
enlightenment comes from the inclusion of Cancel-and-Tighten in
the analysis where phases are observed in actu. The latter ﬁxes
dual variables and depletes the residual network of Type 1 cycles(those formed with negative reduced cost arcs only). Once thisCancel-step is terminated, one must conclude that a Type 2 cycle
comes next thus implying a jump with respect to some lower
h on μh at the beginning of phase h.
bound μˆ CT
Gk at iteration k of the Type 1
Fig. 13 opposes the mean value μCT
cycles canceled in the Cancel-and-Tighten implementation withk ≤ μk proposed by Gauthier et al. (2015,
the lower bound value μˆ CTGk can clearly be
Proposition 14). While the erratic behavior of μCT
observed throughout the solution, a general pattern of increase canbe noted across the phases. The minimum (i.e., optimal) mean cycle value μkG is once again fetched as background information. The
4 900th iteration is worth a look: Again, a small decrease for μkG . In
total, four such occurrences within 7 294 cancellations containedin 357 phases. Furthermore, Cancel-and-Tighten maintains strongly
polynomial runtime despite the usage of Type 1 cycles goingagainst the non-decreasing property of μkG . The strictly increasing
h between phases obtained with the existence oflower bounds μˆ CTa Type 2 cycle marking the end of phase h is indeed where theproperties are established (Gauthier et al., 2015, Theorem 6).
Arc ﬁxing. Strongly polynomial time complexity is achieved bykeeping track of the number of phases through the concept ofarc ﬁxing as seen in the minimum mean cycle-canceling algorithm

J.B. Gauthier et al. / Computers and Operations Research 84 (2017) 16–32

27

kk
Fig. 13. Comparison of μCT, μkG and μˆ CT.

(Gauthier et al., 2015, Propositions 9 and 11). A relaxation of thearc ﬁxing condition is used in Proposition 11, where the two different arc ﬁxing conditions are in line with the complexity point ofviews, that is, O(mnlog n) or O(mn) phases. The proof is straightforward and recuperates the Cancel-and-Tighten proof adaptation
h are shown to mimic the bewhere the lower bound values μˆ CThavior of the true value μhG . Indeed, as mentioned in Gauthier et al.
(2015, p. 131), it is conceivable to rewrite the arc ﬁxing conditions
using a lower bound μˆ h on μhG instead.

5.1. End-phase markers

Proposition 11. Assuming μˆ h is a lower bound for μhG at the
beginning of phase h ≥ 0.

Deﬁnition 3. A phase h ≥ 0 is a sequence of cycle cancellations
terminated whenever a measurable jump is observed in strongly
polynomial metrics, that is, both the factor and the time required
to obtain it are expressed in strongly polynomial time. A phase solution xh is understood as the solution at the beginning of phase h.

• Arc ﬁxing for arc (i, j) ∈ A occurs or has already occurred if|c̄ihj | ≥ −2nμˆ h .
• Implicit arc ﬁxing for arc (i, j) ∈ A occurs or has already occurred
if |c̄i∗j | > −nμˆ h , where c̄i∗j is the reduced cost of arc (i, j) computedwith an optimal set of dual variables.
It is at this point important to distinguish between μkG at
iteration k and strongly polynomial runtime. Indeed, whether one
k or any other cycle-canceling method, the minlooks at μkH or μCT
imum mean cycle value μkG can always be fetched as secondaryk ≤ μk ≤ μk ). If
information (recall that μkH ≤ μkG ≤ μkH:G and μ
ˆ CT
G
CT
strongly polynomial properties are to be established on oscillating
values from one iteration to the next, it appears mandatory torather divert the analysis to a lower bound μˆ h over the phases
instead. In other words, while the portion above μkG (blue) may
be unpredictable, it is the portion below μkG (green) that should
be well-behaved. We are ready to propose an adaptation of thecontraction algorithm. Whether that adaptation is necessary toachieve a strongly polynomial time complexity is an open question
although, since the proof relies on a controllable behavior of μH ,we conjecture that it is.
5. Contraction-Expansion algorithm
Section 5.1 introduces a ﬂexible phase deﬁnition based on
so-called Type 3 cycles that serve as end-phase markers. This isfollowed in Section 5.2 by a discussion of an expansion scheme
which controls the visible and hidden arc sets, that is the content
of the set P, for the proposed algorithm. Section 5.3 argues that applying the contraction on the residual network and expanding thecontracted network using that speciﬁc expansion scheme as thealgorithm progresses produces a strongly polynomial algorithm.
In Section 5.4, we show that the expansion scheme is not uniquesuch that different strategies can be used to improve the algorithm. Finally, computational results are presented in Section 5.5.

Our contribution is to revisit the phase deﬁnition in order toextract the true pertinent information which allows convergencein strongly polynomial time. In this respect, the current phasedeﬁnition is built upon a weak jump condition which waits forthe identiﬁcation of a Type 2 cycle as the minimum mean cycle
to conﬁrm a jump on μhG at the end of phase h. Let us propose amore ﬂexible deﬁnition.

The factor proposed in the following cycle Type obviously
satisﬁes the strongly polynomial requirement whereas the timerequirement is shown in Proposition 12.
Deﬁnition 4. Let xh , h ≥ 0, be a non-optimal phase solution. At
an iteration t within the phase h, let us call a cycle Wt on G(xt )a Type 3 (negative) cycle if its underestimated mean cost μˆ t (withμˆ t ≤ μtG ≤ μt ) produces the measurable jumpμˆ t ≥ 1 −1|W t |


μhG .

(15)

Proposition 12. In MMCC, a Type 3 cycle is identiﬁed within at mostm cancellations.
Proof. Consider the phase h. The proof is trivial and connectswith the existence of a Type 2 cycle, say W , which must haveμG = μˆ  ≥ (1 − 1/|W  |)μhG . Therefore, the iteration t either happens simultaneously to  or earlier. 
Observe that whether the jump factor is obtained using an
actual Type 2 cycle or not is irrelevant: A phase is completed inaccordance with Deﬁnition 3. In other words, the only important
aspect of the Type 2 cycle is the measurable jump it procures on
μhG , a tactic which can incidentally also be veriﬁed against a lower
bound according to Proposition 11, that is,μˆ t ≥ 1 −1|W t |

μˆ h .

(16)h in Cancel-and-Tighten (Gauthier et al.,This is true for μˆh = μˆ CT
2015, Theorem 3) and is also used with μˆ h = μhH in the proposed
Contraction-Expansion algorithm.

28

J.B. Gauthier et al. / Computers and Operations Research 84 (2017) 16–32

Fig. 14. Contraction-Expansion algorithm.

5.2. Expansion scheme

5.3. Complexity analysis

The focus on phases should by now be realized by the reader.
Fig. 14 presents a pseudo-code for the Contraction-Expansion algorithm whereas the following paragraphs explain how the proposed
expansion scheme constrains the latter into producing such phases.
Let x0 be a non-optimal solution at iteration k = 0. Note thatx0 is also a non-optimal phase solution ensuring a lower boundμˆ 0 ≤ μ0G . Let the set P0 = F 0 . By Proposition 6, the extended cycle0WH:
is directed on the residual network G(x0 ) and the step size isGstrictly positive.
0
When WH:is canceled, the aftermath is hard to predict but
G0
one thing is for certain: Only the arcs part of WH:
are affected.
G0Some arcs that were free in x remain free in x1 in the next
iteration while all the other arcs have changed status either fromrestricted lower to restricted upper and vice versa, from free torestricted or from restricted to free for a total of six possibilities.
Four of these end up with a new status restricted whereas two
end up with a free status. Intuitively speaking, since the contraction happens around the free variables and the non-decreasing
μkG property is lost when the contraction is systematically applied,
let us concentrate on controlling these newly freed variables in
x1 , i.e., F1 F0 . By applying coerced degeneracy on these speciﬁcvariables, the hidden set then contains only those variables thatwere already free in x0 and are still free in x1 . This amounts toselecting the set P 1 = F 1 ∩ F 0 at iteration k = 1.
From there, the idea is simple, repeatedly apply cycle cancelk
lation with the extended cycle WH:
at iteration k and expandG
the contracted network using the intersection of variables free inboth the previous and new solution until the algorithm reaches
a Type 3 cycle, that is, a negative cycle producing a jump factork | ). The phase h = 0 is then terminated andof at least (1 − 1/|WH:G
a new phase begins. Since we may at this point re-apply the fullcontraction, it is worthwhile to eliminate any cycle of free arcs
from the phase solution to maximize the contraction beneﬁt.

The complexity analysis of the Contraction-Expansion algorithm
revolves around bringing the solution process of the latter to termswith MMCC’s behavior. As testiﬁed by Cancel-and-Tighten, theseterms refer to the strictly increasing phase markers. Fig. 15 might
help to get a feel for this endeavor. It displays the value of μkHusing the minimal expansion scheme P k+1 = F k+1 ∩ F k where eachtrail corresponds to a phase during the solution of Instance 1.
The upcoming analysis focuses on the outer loop, that is, we
examine the time required to reach the end of a phase and thetotal number of such phases.
Proposition 13. The Contraction-Expansion algorithm completes aphase in O(m2 n) time.
Proof. The proof is threefold with respect to the expansion
scheme portion of the algorithm. We prove that 1) the nondecreasing property of μkH during the phase is maintained, 2) atmost O(n) contracted cycle cancellations are required to reach acontracted network equivalent to the residual network, and 3)
ﬁnding the minimum mean reduced cost cycle on the contracted
network requires O(mn) time.
Gauthier et al. (2015, Proposition 4) state that a cycle cancellation on the residual network G(xk ) cannot introduce a minimum
mean cycle in G(xk+1 ) yielding μkG+1 < μkG in the following iteration. Here is an adaptation for the expansion scheme. Consider
the contracted network H(Pk , xk ) and an optimal contracted cycle
WHk of mean reduced cost μkH along with an optimal set of dual
variables πk ,  ∈ NP (xk ), on the root nodes (Proposition 9). From
there, ﬁx the root nodes’ dual variables to these new values andupdate the remaining dual variables accordingly (Procedure 1). Thekkkresidual network G(xk ) hence satisﬁes d¯i j − πR( i ) + π R ( j ) ≥ μH , ( i ,kk¯j) ∈ A(x ), and in particular di j = 0, (i, j) ∈ HP (x ). Upon cancelingkthe expanded cycle WH:
and obtaining xk+1 , every new residGk+1
ual arc in G(x) either has a reduced cost of 0 or −μkH , i.e.,J.B. Gauthier et al. / Computers and Operations Research 84 (2017) 16–32

29

Fig. 15. Evolution of μkH with the Contraction-Expansion algorithm.

kkkk+1 ). Observe that every arc ind¯i j − πR(i ) + πR( j ) ≥ μH , (i, j ) ∈ A (xF k ∩ F k+1

already has a reduced cost of 0 such that the arc costtransfer policy already holds. The contraction is readily available
with every remaining visible arc evaluated at a reduced costgreater than or equal to μkH . The mean cost of the next contracted
cycle is then at least μkH .
Recall that the pricing problem in MMCC can be derived fromthe pricing problem (9) by making visible all residual arcs, that is,by setting P = ∅. Then, observe that the set P is updated by intersecting the sets of free variables of the previous solution with thecurrent one such that its size either stays the same or decreases
at each cancellation. Assuming the initial phase solution x0 is abasic solution, at most |F0 | ∈ O(n) cancellations yielding a decrease
are then obviously required to reach P = ∅. Furthermore, when acycle is canceled without modifying the set of free arc variables,
it means that one or several restricted variables change fromone bound to another. Unfortunately, this kind of phenomenon
gives meaning to the worst case complexity bound observed inpathological instances. An iteration counter limiting the numberof cancellations prior to reaching the residual network to at mostn is put in place to make the move directly should it be necessary.
Trivially, at most n cycle cancellations allow the expansion scheme
to reach a contracted network equivalent to the residual network.
From there, a Type 3 cycle is ensured within an additional m cycle
cancellations (Proposition 12). To sum up, at most n + m ∈ O(m )cycle cancellations are required to meet the required jump.
Finally, solving the pricing problem (9), which can still beaccomplished in O(mn) time using dynamic programming (Karp,
1978), dominates all the other operations performed at every iteration. Indeed, the data structure is maintained in O(m) time while
the extended cycle extraction, step size computation and solution
update is done in O(n) time. A phase is ultimately completed in at
most O(m) iterations each one requiring at most O(mn) for a total
phase runtime of O(m2 n). 
Theorem 1. The Contraction-Expansion algorithm runs in O(m3 n2 )strongly polynomial time.
Proof. We show that at most O(mn) phases can occur in accordance with Gauthier et al. (2015, Theorem 3). Since μkH ≤ μkG isindeed a lower bound for the true minimum mean cycle costvalue by Proposition 10, this is also true for any phase solution.
k | ), the phase h ends and the
As soon as μkH ≥ μˆ h (1 − 1/|WH:Glower bound is increased to μ
ˆ h+1 := μkH . By Proposition 11, arc
ﬁxing occurs on the lower bound value μˆ h as well such thatthe phase time complexities are valid. By Proposition 13, sinceeach phase runs in O(m2 n), the compound time is obtained. An

initial valid lower bound for μ0G can be trivially obtained withμˆ 0 := − max(i, j )∈A |ci j |. 

The proposed phase deﬁnition along with the Type 3 cycle not
only satisfy theoretical properties of the strongly polynomial timecomplexity of MMCC, they also express very practical observations.
The hope is that not only the Type 3 cycle occurs much earlier
than after m cancellations, it also happens while the phase still
exploits the contracted network.
5.4. Alternative end-phase markers and expansion schemes
So long as strongly polynomial phase time is maintained,
alternative expansion schemes may be used. Figs. 16–17 show theevolution of μkH for two simple variations. The ﬁrst variation (cycle
expansion) updates the set P with a more aggressive reduction,i.e., a faster expansion. The update writes as

P k+1 = F k+1 ∩ F kkk\ {(i, j ) ∈ A | (i, j ) ∈ WH:
G or ( j, i ) ∈ WH:G }(17)kcycle WH:
G

such that all arcs of the expandedthat are still free(thus common to Fk and F k+1 ) are also removed from P k+1 . The
second variation (cycle expansion & loops) uses the same update
and also uses a post cycle-cancellation heuristic which cancels
loops derived from all compatible variables with negative rooted
costs as well.
For anyone familiar with successful divide-and-conquer methods, the solution speed typically beneﬁts from the decomposition
at a higher rate than the cost of the latter. While postponing theend of a phase with a less aggressive expansion scheme appears toagree with this statement, expanding the contracted network fasteralso means that end-phase markers are reached faster at which
point the full contraction is re-applied. It seems that reaching said
end-phase markers as fast as possible is of particular interest.
Speaking of end-phase markers, a measurable jump could beestablished using alternative Type 3 cycle deﬁnitions. Consider for
instance the modiﬁcation of (16) as follows:μˆ t ≥ 1 −1|W t |

μˆ h ≥ 1 −1
maxk∈h |Wk|μˆ h ≥ 1 −1n
μˆ h .
(18)

where k ∈ h is understood as an iteration k within the phaseh. While the last criterion is reminiscing of Theorem 2 which
contents itself with the same jump every time thus obtaining theO(mnlog n) phases, the second criterion also tracks a cycle lengthand compromises on the desired jump. The latter is in fact thecriterion used in all plots of the Contraction-Expansion algorithm.

30

J.B. Gauthier et al. / Computers and Operations Research 84 (2017) 16–32

Fig. 16. Evolution of μkH with cycle expansion.

Fig. 17. Evolution of μkH with cycle expansion and cancellation of negative loops.

5.5. Computational experiments
This section is separated in two parts. The ﬁrst details computational results with respect to three speciﬁc instances and serves
to show CE indeed behaves as expected when identifying smaller
cycles on the contracted network. The second complements intuitive assertions made in the ﬁrst part with a computational proﬁleon benchmark instances from DIMACS. The latter contains several
network problems of different nature codenamed by family GOTO,GRIDGEN, GRIDGRAPH, NETGEN, ROAD, and VISION. Each family
is further divided into subcategories of different sizes each one
containing ﬁve versions A-E which vary only in the data set.
Algorithms are implemented in C++. All results are obtained using
an Intel i7-4770K@3.50GHz with 16GB of RAM, running
Windows 10.
Table 1 resumes the content of plots displayed in Figs. 15–
17 along with their computational times as well as averages ofcycle sizes and induced step sizes. The performances of MMCC
and CT are added for order comparison purposes whereas results
regarding CE also incorporate an average contraction level given
by α where α k = |NP (xk )|/n, k ≥ 0.
First of all, there is no denying that Cancel-and-Tighten isorders of magnitude faster than MMCC. The proposed ContractionExpansion algorithm sits somewhere in between although it isclear that the contraction boosts the speed of fetching negative cycles by a signiﬁcant amount even more so as the size ofthe instance gets larger. The same feel is palpable across otherbenchmark instances such as the signiﬁcantly larger problem gridgen_sr_13a from the GRIDGEN family which contains 8 191 nodes
and 745 381 arcs. It is important to understand that the beneﬁt of
the contraction comes from the degeneracy observed in encountered solutions. For instance, the problem road_paths_05_WI_a(n = 519 157 and m = 1 266 876) from the ROAD family is structured with uniform capacities at 1 such that there are never
any free variables, hence no contraction. We therefore just get
penalized with the contraction computational overhead.
The remainder of this section looks at the version A of the
GRIDGRAPH, NETGEN, and GRIDGEN families, in the process,bringing further testimonies to the previous claims. The computational results obtained on the ﬁrst family are presented in moredetails as follows. For each of the three categories of problems(long, square, wide), the number of arcs is roughly m ≈ 2n
whereas the node size n ranges from 28 to 217 , where the poweris given by the instance diﬃculty index i. Table 2 displays the average contraction level and CPU time CPUCE for CE. As expected thelatter increases with the instance size, yet the average contractionlevel of about 60% appears stable across this family. The middle
columns indicate the relative CPU time calculated with respect to
MMCC’s as β := CPUCE /CPUMMCC . For instance, CE requires 2 619 s
to solve grid_square_16 meaning that MMCC took 2 619/0.21 ≈12 471 s to terminate. The relative computing times β are plottedin Fig. 18 with respect to the instance diﬃculty index i. Onceagain, the relative performance of CE increases with the instancesize.
The NETGEN family is broken down into ﬁve categories (8,deg, lo_8, lo_sr, sr). The category deg is stable in node size at
n = 4 096 although it increases exponentially in terms of arcs
m = n2i . All other problems range from sizes 28 to 219 in termsof nodes n, and in terms of arcs either m = 8n for 8 and lo_8
or m = n2i/2 for lo_sr and sr problems. The mean of the averagecontraction levels α is about 93% although it is signiﬁcantly closer
to 100% for problems in the lo_8 and lo_sr categories. Said lack ofcontraction is apparent in Fig. 19 where β is mostly present above
1 represented by the horizontal dotted line.

J.B. Gauthier et al. / Computers and Operations Research 84 (2017) 16–32

31

Table 1
Computational results for variations of CE vs. MMCC and CT.
Solution

CPU (sec)k / cycles

h

|W |

ρ

α

Instance 1
MMCC
CE - minimal expansion
CE - cycle expansion
CE - cycle expansion/loops
CT

91.52
87.59
63.98
51.84
0.82

4 296
7 256
4 686
3 449/9 499
7 294

90
151
144
142
357

27.11
5.46
5.19
7.40/3.32
21.62

1.30
1.30
1.70
1.61/2.23
1.37

–
0.32
0.47
0.58
–

gridgen_sr_13a
MMCC
CE - minimal expansion
CE - cycle expansion
CE - cycle expansion/loopsCT

15 629.40
6 709.96
2 435.97
1 906.97
36.17

20 89325 2919 202
6 034/23 429
27 473

289
301
236
224
1 482

72.85
5.93
4.70
7.13/2.58
49.80

7.45
9.70
19.18
18.53/36.01
7.28

–
0.56
0.62
0.70
–

road_paths_05_WI_a
MMCC
CE - cycle expansion
CT

9 552.52
10 355.30
2 389.51

1 119
1 122
3 248

300
348
18 647

82.27
82.20
34.85

1.00
1.00
1.00

–
1.00
–

Table 2Computational results for CE on the GRIDGRAPH family, version A.
longsquarewide

iCPUCE (sec)βαCPUCE (sec)βαCPUCE (sec)

β

α

8
9
10
11
12
13
14
15
16
17

0.01
0.02
0.08
0.18
0.59
1.51
5.45
8.51
23.15
402.98

0.65
0.44
0.31
0.26
0.19
0.12
0.16
0.08
0.04
0.13

0.64
0.59
0.56
0.61
0.56
0.50
0.45
0.51
0.45
0.45

0.01
0.04
0.16
0.64
2.95
11.96
59.18
325.442 619.07
17 890.40

0.60
0.43
0.36
0.30
0.29
0.26
0.26
0.22
0.21
0.27

0.64
0.56
0.58
0.56
0.58
0.57
0.58
0.60
0.60
0.61

0.01
0.05
0.21
0.93
4.43
19.72
89.83
476.192 911.39
15 944.70

0.60
0.50
0.41
0.44
0.38
0.39
0.41
0.38
0.31
0.34

0.64
0.59
0.61
0.64
0.64
0.64
0.63
0.63
0.63
0.63

Fig. 18. Relative CPU time for CE when solving the GRIDGRAPH family, version A.

Fig. 19. Relative CPU time for CE when solving the NETGEN family, version A.

32

J.B. Gauthier et al. / Computers and Operations Research 84 (2017) 16–32

Fig. 20. Relative CPU time for CE when solving the GRIDGEN family, version A.

The GRIDGEN family is partitioned into three categories ofproblems (8, deg, sr). The category deg is once again stable atn = 4 097 and increases exponentially in terms of arcs m = n2i . The
two other categories range in node sizes n from 28 to 218 whereas
the number of arcs is given by m = 8n for the former category andm = n2i/2 for the latter. The mean of the average contraction levels
α is about 90% for these problems. Fig. 20 displays the experimentoutcome on this family and still features increasing relative performance as the instances grow in diﬃculty.
While we certainly did not cover the complete benchmark suite,
this paper does not pretend to propose a competitive algorithm
just yet. As a case in point, even larger instances up to 222 nodes
are available but MMCC’s solution proved to be far too demanding. We also believe the omitted instance versions and/or more results from other families would not signiﬁcantly contribute to thisframework. Finally, we recall that PS is a top performer for the solution of minimum cost ﬂow problems despite the fact that thevast majority of pivots performed, read over 70%, are degenerate(Ahuja et al., 1993, Fig. 18.7). With this in mind, the goal of thispaper is rather to understand better the mechanics of algorithms
with strictly positive step size improvements.
6. Conclusion
We start with a note addressed to users of the primal network
simplex algorithm. Observe that the spirit of the minimum mean
cycle-canceling algorithm is tangibly similar to that of PS. The
pricing step is home to the optimality certiﬁcate whereby thelatter is acquired unless an improving direction is identiﬁed. Itturns out this is not all that surprising since it has been shownthat PS and MMCC belong to (even constitute extreme cases of) amore generic contraction algorithm.
A variety of special cases inducing strictly positive step sizes
is also identiﬁed. By combining these with results from MMCC,
a strongly polynomial Contraction-Expansion algorithm which
behaves much better than the former, especially as the instance’s
size increases, comes to life. The reader is carried around thisbirth process by opposing the behavior of MMCC and Cancel-andTighten in a computational study. Both iteration and phase bases
are illustrated, although a strong emphasis on the latter concept issystematically done thus providing an alternative way of showing
strongly polynomial runtime.
Such a property can also be observed in this framework whenusing partial contraction. The latter is obtained by modifying the
choice of hidden arcs as the algorithm progresses. The selection ismade in such a way that it actually corresponds to an expansion of
the contracted network. Furthermore, the proposed ContractionExpansion algorithm sticks to practical observations otherwise
overlooked in MMCC. As such, phase markers are veriﬁed algorithmically rather than just existing for theoretical purposes. Iteven recuperates the idea that not all jumps are created equal
thus underlining the important aspect of Type 2 cycles, namely
the measurable jump. The Type 3 cycle is born. The strongly
polynomial argument uses both phases and Type 3 cycles on topof which the convergence of the original optimality parameter isneglected in favor of a lower bound.
Furthermore, while it is true that the time complexity is not
improved with respect to MMCC, we believe it is more interesting
to note that the same time complexity is achieved despite a more
complicated algorithm. Indeed, strongly polynomial algorithms
might beneﬁt from another look that aims to combine their timecomplexity with practical observations that make them behave
more eﬃciently.
There seems to be some arbitrage to be done between trying
to meet optimality conditions in a more aggressive manner andthe work required to do so. By contenting itself with a suﬃcientcondition, a signiﬁcant proportion of cancellations performed in PS
are degenerate whereas MMCC uses a rule whose computational
footprint is too high. That being said, the contracted network is
closer to the spirit of an oracle than is the residual network. By
this, we mean that it matters not to see all directed negative costcycles so long as at least one can be detected thus allowing to
improve and repeat.
Acknowledgements
Jacques Desrosiers acknowledges the Natural Sciences andEngineering Research Council of Canada for its ﬁnancial support.
References
Ahuja, R.K., Magnanti, T.L., Orlin, J.B., 1993. Network Flows: Theory, Algorithms, andApplications. Prentice Hall, Upper Saddle River, NJ, USA.
Cormen, T.H., Leiserson, C.E., Rivest, R.L., Stein, C., 2009. Introduction to Algorithms.
The MIT Press, Cambridge, MA, USA and London, England 3rd edition.
Dantzig, G.B., Thapa, M.N., 2003. Linear programming 2: Theory and extensions.
In: Springer Series in Operations Research and Financial Engineering (Book 2).
Springer, New York, NY, USA.
DIMACS,. Network ﬂows and matching: First DIMACS implementation challenge,
1990–1991. URL ftp://dimacs.rutgers.edu/pub/netﬂow.
Gauthier, J.B., Desrosiers, J., Lübbecke, M.E., 2015. About the minimum mean cyclecanceling algorithm. Discrete Appl. Math. 196, 115–134. Advances in Combinatorial Optimization doi: 10.1016/j.dam.2014.07.005
Goldberg, A.V., Tarjan, R.E., 1989. Finding minimum-cost circulations by cancelingnegative cycles. J. ACM 36 (4), 873–886. doi:10.1145/76359.76368.
Karp, R.M., 1978. A characterization of the minimum cycle mean in a digraph. Discrete Math. 23 (3), 309–311. doi:10.1016/0012-365X(78)90011-0.
Klein, M., 1967. A primal method for minimal cost ﬂows with applications to theassignment and transportation problems. Manage. Sci. 14 (3), 205–220.
Kovács, P., 2015. Minimum-cost ﬂow algorithms: an experimental evaluation. Optim.
Methods Softw. 30 (1), 94–127. doi:10.1080/10556788.2014.895828.
Marsten, R.E., Saltzman, M.J., Shanno, D.F., Pierce, G.S., Ballintijn, J.F., 1989. Implementation of a dual aﬃne interior point algorithm for linear programming.
ORSA J. Comput. doi:10.1287/ijoc.1.4.287.
Radzik, T., Goldberg, A.V., 1994. Tight bounds on the number of minimum-mean cycle cancellations and related results. Algorithmica 11 (3), 226–242. doi:10.1007/
BF01240734.