c 2003 Society for Industrial and Applied Mathematics


SIAM J. COMPUT.
Vol. 32, No. 4, pp. 833–840

A FASTER SCALING ALGORITHM FOR MINIMIZINGSUBMODULAR FUNCTIONS∗SATORU IWATA†
Abstract. Combinatorial strongly polynomial algorithms for minimizing submodular functions
have been developed by Iwata, Fleischer, and Fujishige (IFF) and by Schrijver. The IFF algorithm
employs a scaling scheme for submodular functions, whereas Schrijver’s algorithm achieves strongly
polynomial bound with the aid of distance labeling. Subsequently, Fleischer and Iwata have describeda push/relabel version of Schrijver’s algorithm to improve its time complexity. This paper combines
the scaling scheme with the push/relabel framework to yield a faster combinatorial algorithm for
submodular function minimization. The resulting algorithm improves over the previously best known
bound by essentially a linear factor in the size of the underlying ground set.
Key words. submodular function, discrete optimization, algorithm
AMS subject classiﬁcation. 90C27
DOI. 10.1137/S0097539701397813

1. Introduction. Let V be a ﬁnite nonempty set of cardinality n. A set function
f on V is submodular if it satisﬁes
f (X) + f (Y ) ≥ f (X ∩ Y ) + f (X ∪ Y )∀X, Y ⊆ V.

Submodular functions are discrete analogues of convex functions [14]. Examples ofsubmodular functions include cut capacity functions, matroid rank functions, andentropy functions.
The ﬁrst polynomial-time algorithm for submodular function minimization is due
to Grötschel, Lovász, and Schrijver [9]. A strongly polynomial algorithm has also
been described by Grötschel, Lovász, and Schrijver [10]. These algorithms rely on theellipsoid method, which is not eﬃcient in practice.
Recently, combinatorial strongly polynomial algorithms have been developed by
Iwata, Fleischer, and Fujishige (IFF) [13] and by Schrijver [16]. Both of these algorithms build on works of Cunningham [2, 3]. The IFF algorithm employs a scaling
scheme developed in capacity scaling algorithms for the submodular ﬂow problem[7, 11]. In contrast, Schrijver [16] directly achieves a strongly polynomial bound by
introducing a novel subroutine in the framework of lexicographic augmentation. Subsequently, Fleischer and Iwata [5, 6] have described a push/relabel algorithm using
Schrijver’s subroutine to improve the running time bound. In this paper, we combine
the scaling scheme with the push/relabel technique to yield a faster combinatorial
algorithm.
Let γ denote the time required for computing the function value of f and M
denote the maximum absolute value of f . The IFF scaling algorithm minimizes anintegral submodular function in O(n5 γ log M ) time. The strongly polynomial version∗ Received by the editors November 12, 2001; accepted for publication (in revised form) March 31,
2003; published electronically June 10, 2003. This research was supported in part by the Sumitomo
Foundation and a Grant-in-Aid for Scientiﬁc Research of the Ministry of Education, Science, Sports
and Culture of Japan.
http://www.siam.org/journals/sicomp/32-4/39781.html† Department of Mathematical Informatics, University of Tokyo, Tokyo 113-8656, Japan (iwata@mist.i.u-tokyo.ac.jp).

833

834

SATORU IWATA

runs in O(n7 γ log n) time, whereas an improved variant of Schrijver’s algorithm runsin O(n7 γ + n8 ) time [6].
The time complexity of our new scaling algorithm is O((n4 γ + n5 ) log M ). Since
the function evaluation oracle has to identify an arbitrary subset of V as its argument,
it is natural to assume γ is at least linear in n. With this assumption, the newalgorithm is faster than the IFF algorithm by a factor of n. The strongly polynomial
version of the new scaling algorithm runs in O((n6 γ + n7 ) log n) time. This is an
improvement over the previous best bound by essentially a linear factor in n.
These combinatorial algorithms perform multiplications and divisions, although
the problem of submodular function minimization does not involve those operations.
Schrijver [16] asks if one can minimize submodular functions in strongly polynomial
time using only additions, subtractions, comparisons, and oracle calls for the function
values. Such an algorithm is called “fully combinatorial.” A very recent paper [12]settles this problem by developing a fully combinatorial variant of the IFF algorithm.
Similarly, we can implement the strongly polynomial version of our scaling algorithm
in a fully combinatorial manner. The resulting algorithm runs in O(n8 γ log2 n) time,
improving the previous O(n9 γ log2 n) bound by a factor of n.
This paper is organized as follows. Section 2 provides preliminaries on submodular
functions. In section 3, we describe the new scaling algorithm. Section 4 is devoted
to its complexity analysis. Finally, in section 5, we discuss its extensions as well as afully combinatorial implementation.
2. Preliminary. This section provides preliminaries on submodular functions.
See [8, 14] for more details and general background.

For a vector x ∈ RV and a subset Y ⊆ V , we denote x(Y ) = u∈Y x(u). We
also denote x− the vector in RV with x− (u) = min{x(u), 0}. For each u ∈ V , let χudenote the vector in RV with χu (u) = 1 and χu (v) = 0 for v ∈ V \{u}.
For a submodular function f : 2V → R with f (∅) = 0, we consider the basepolyhedron
B(f ) = {x | x ∈ RV , x(V ) = f (V ), ∀Y ⊆ V : x(Y ) ≤ f (Y )}.
A vector in B(f ) is called a base. In particular, an extreme point of B(f ) is called
an extreme base. An extreme base can be computed by the greedy algorithm ofEdmonds [4] and Shapley [17] as follows.
Let L = (v1 , . . . , vn ) be a linear ordering of V . For any vj ∈ V , we denote
L(vj ) = {v1 , . . . , vj }. The greedy algorithm with respect to L generates an extreme
base y ∈ B(f ) by
(2.1)

y(u) := f (L(u)) − f (L(u)\{u}).

Conversely, any extreme base can be obtained in this way with an appropriate linear
ordering.
Lemma 2.1. Let Q and R be disjoint subsets of V such that Q ∪ R forms aninterval in L. Let L be the linear ordering obtained from L by moving Q to the place
immediately after R without changing the orderings in Q and in R. Then the extreme
base y  generated by L satisﬁes y  (q) ≤ y(q) for q ∈ Q and y  (r) ≥ y(r) for r ∈ R.
Proof. For any q ∈ Q, we have L (q) ⊇ L(q). Therefore, the submodularity of f
implies y  (q) = f (L (q)) − f (L (q)\{q}) ≤ f (L(q)) − f (L(q)\{q}) = y(q). Similarly,L (r) ⊆ L(r) holds for r ∈ R. Then it follows from the submodularity of f thaty  (r) = f (L (r)) − f (L (r)\{r}) ≥ f (L(r)) − f (L(r)\{r}) = y(r).

SUBMODULAR FUNCTION MINIMIZATION835

For any base x ∈ B(f ) and any subset Y ⊆ V , we have x− (V ) ≤ x(Y ) ≤ f (Y ).
The following theorem shows that these inequalities are in fact tight for appropriately
chosen x and Y .
Theorem 2.2. For a submodular function f : 2V → R, we havemax{x− (V ) | x ∈ B(f )} = min{f (Y ) | Y ⊆ V }.
Moreover, if f is integer-valued, then the maximizer x can be chosen from among
integral bases.
This theorem is immediate from the vector reduction theorem on polymatroids
due to Edmonds [4]. It has motivated combinatorial algorithms for minimizing submodular functions.
3. A scaling algorithm. This section presents a new scaling algorithm forminimizing an integral submodular function f : 2V → Z.
The algorithm consists of scaling phases with a scale parameter δ ≥ 0. It keeps
a set of linear orderings {Li | i ∈ I} of the vertices in V . We denote v i u if
v precedes u in Li or v = u. Each linear ordering Li generates an extreme baseyi ∈ B(f ) by the greedy algorithm.
The algorithm also keeps a base x ∈ B(f ) as aconvex combination x = i∈I λi yi of the extreme bases. Initially, I = {0} with anarbitrary linear ordering L0 and λ0 = 1.
Furthermore, the algorithm works with a ﬂow in the complete directed graph on
the vertex set V . The ﬂow is represented as a skew-symmetric function ϕ : V ×V → R.
Each arc capacity is equal to δ. Namely, ϕ(u, v) + ϕ(v, u) = 0 and −δ ≤ ϕ(u, v) ≤ δhold for any pair of vertices u, v ∈ V . The boundary ∂ϕ is deﬁned by ∂ϕ(u) =v∈V ϕ(u, v) for u ∈ V . Initially, ϕ(u, v) = 0 for any u, v ∈ V .
Each scaling phase aims at increasing z − (V ) for z = x + ∂ϕ. Given a ﬂow
ϕ, the procedure constructs an auxiliary directed graph Gϕ = (V, Aϕ ) with arc set
Aϕ = {(u, v) | u = v, ϕ(u, v) ≤ 0}. Let S = {v | z(v) ≤ −δ} and T = {v | z(v) ≥ δ}.
A directed path in Gϕ from S to T is called an augmenting path.
Each scaling phase also keeps a valid labeling d. A labeling d : V → Z is valid ifd(u) = 0 for u ∈ S and v i u implies d(v) ≤ d(u) + 1. A valid labeling d(v) serves asa lower bound on the number of arcs from S to v in the directed graph GI = (V, AI )with the arc set AI = {(u, v) | ∃i ∈ I, v i u}.
If there is an augmenting path P , the algorithm augments the ﬂow ϕ along P
by ϕ(u, v) := ϕ(u, v) + δ and ϕ(v, u) := ϕ(v, u) − δ for each arc (u, v) in P . This
procedure is referred to as Augment(ϕ, P ). As a result of Augment(ϕ, P ), the initialvertex s of P may get rid of S and no new vertex joins S. Thus Augment(ϕ, P ) does
not violate the validity of d.
Let W be the set of vertices reachable from S in Gϕ , and let Z be the set ofvertices that attains the minimum labeling in V \W . A pair (u, v) of u ∈ W andv ∈ Z is called active for i ∈ I if v is the ﬁrst vertex of Z in Li and u is the last vertex
in Li with v i u and d(v) = d(u) + 1. A triple (i, u, v) is also called active if (u, v)
is active for i ∈ I. The procedure Multiple-Exchange(i, u, v) is applicable to an active
triple (i, u, v).
For an active triple (i, u, v), the set of vertices from v to u in Li is called an active
interval. The active interval is divided into Q and R by Q = {w | w ∈ W, v ≺i w i u}
and R = {w | w ∈ V \W, v i w ≺i u}.
The procedure Multiple-Exchange(i, u, v) moves the vertices in R to the place
immediately after u in Li , without changing the ordering in Q and in R. Then it
computes an extreme base yi generated by the new Li . By Lemma 2.1, this results

836

SATORU IWATA

in yi (q) ≥ yi◦ (q) for q ∈ Q and yi (r) ≤ yi◦ (r) for r ∈ R, where yi◦ denotes the previous
yi .
Consider a complete bipartite graphwith the vertex sets Q◦ and R. The algorithm
ﬁnds
a
ﬂow
ξ
:
Q×R
→
R
such
that
+
r∈R ξ(q, r) = yi (q)−yi (q) for each q ∈ Q and

◦
ξ(q,
r)
=
y
(r)
−
y
(r)
for
each
r
∈
R. Such a ﬂow can be obtained easily byiiq∈Q
the so-called northwest corner rule. Then the procedure computes α = min{λi , δ/β}with β = max{ξ(q, r) | q ∈ Q, r ∈ R} and moves x by x := x + α(yi − yi◦ ). In orderto keep z invariant, the procedure adjusts the ﬂow ϕ by ϕ(q, r) := ϕ(q, r) − αξ(q, r)and ϕ(r, q) := ϕ(r, q) + αξ(q, r) for every (q, r) ∈ Q × R. The resulting ϕ satisﬁes thecapacity constraints due to the choice of α, and the vertices in W remain reachable
from S in Gϕ .
If α = λi , Multiple-Exchange(i, u, v) is called saturating. Otherwise, it is called
nonsaturating. In a nonsaturating Multiple-Exchange(i, u, v), a new index k is added
to I. The associated linear ordering Lk is the previous Li . The coeﬃcient λk isdetermined by λk := λi − α, and then λi is replaced by λi := α. Thus the algorithm
continues to keep x as a convex combination x = i∈I λi yi .
Suppose the labeling d is valid before the algorithm applies Multiple-Exchange toan active triple (i, u, v). For any vertex w in the active interval, d(v) ≤ d(w) + 1 andd(w) ≤ d(u)+1 hold. These inequalities and d(v) = d(u)+1 imply d(v) ≤ d(w) ≤ d(u).
Note that d(v) ≤ d(r) holds for any r ∈ R ⊆ V \W . Hence we have d(r) = d(v) for
any r ∈ R. If Multiple-Exchange(i, u, v) adds a new arc (s, t) to AI , then s ∈ Q andt ∈ R. Therefore, we have d(t) = d(v) ≤ d(s) + 1. Thus Multiple-Exchange(i, u, v)
does not violate the validity of d.
Let h denote the number of vertices in the active interval. The number of function
evaluations required for computing the new extreme base yi by the greedy algorithm
is at most h. The northwest corner rule can be implemented to run in O(h) time,
and the number of arcs (q, r) with ξ(q, r) > 0 is at most h − 1. Thus the total timecomplexity of Multiple-Exchange(i, u, v) is O(hγ).
If there is no active triple, the algorithm applies Relabel to each v ∈ Z. The
procedure Relabel(v) increments d(v) by one. Then the labeling d remains valid.
The number of extreme bases in the expression of x increases by one as a resultof nonsaturating Multiple-Exchange. In order to reduce the complexity, the algorithmoccasionally applies a procedure Reduce(x, I) that computes an expression of x as aconvex combination of aﬃnely independent extreme bases chosen from the currently
used ones. This computation takes O(n2 |I|) time with the aid of Gaussian elimination.
We are now ready to describe the new scaling algorithm.
Step 0: Let L0 be an arbitrary linear ordering. Compute an extreme base y0 by thegreedy algorithm with respect to L0 . Put x := y0 , λ0 := 1, I := {0}, andδ := |x− (V )|/n2 .
Step 1: Put d(v) := 0 for v ∈ V , and ϕ(u, v) := 0 for u, v ∈ V .
Step 2: Put S := {v | z(v) ≤ −δ} and T := {v | z(v) ≥ δ}, where z = x + ∂ϕ. Let
W be the set of vertices reachable from S in Gϕ .
Step 3: If there is an augmenting path P , then do the following.
(3-1) Apply Augment(ϕ, P ).
(3-2) Apply Reduce(x, I).
(3-3) Go to Step 2.
Step 4: Compute - := min{d(v) | v ∈ V \W } and put Z := {v | v ∈ V \W, d(v) = -}.
If - < n, then do the following.
(4-1) If there is an active triple (i, u, v), then apply Multiple-Exchange(i, u, v).
(4-2) Otherwise, apply Relabel(v) for each v ∈ Z.

SUBMODULAR FUNCTION MINIMIZATION

837(4-3) Go to Step 2.
Step 5: Determine the set X of vertices reachable from S in GI . If δ ≥ 1/n2 , thenapply Reduce(x, I), δ := δ/2, and go to Step 1.
We now show that the last set X obtained by the scaling algorithm is a minimizerof f .
Lemma 3.1. At the end of each scaling phase, z − (V ) ≥ f (X) − n(n + 1)δ/2.
Proof. At the end of each scaling phase, d(v) = n for every v ∈ V \W . Since d(v)is a lower bound on the number of arcs from S to v in GI , this means there is no
directed path from S to V \W in GI . Thus we have X ⊆ W ⊆ V \T , which impliesz(v) ≤ δ for v ∈ X. It follows from S ⊆ X that z(v) ≥ −δ for v ∈ V \X. Since there
is no arcin GI emanating from X, we have yi (X) = f (X) for each i ∈ I, and hencex(X) = i∈I λi yi (X) = f (X). We also have ∂ϕ(X) ≥ −δ |X|·|V \X| ≥ −n(n−1)δ/2.
Therefore, we have z − (V ) = z − (X) + z − (V \X) ≥ z(X) − δ|X| − δ|V \X| = x(X) +∂ϕ(X) − nδ ≥ f (X) − n(n + 1)δ/2.
Lemma 3.2. At the end of each scaling phase, x− (V ) ≥ f (X) − n2 δ.
Proof. The set Y = {v | x(v) < 0} satisﬁes x− (V ) = x(Y ) = z(Y ) − ∂ϕ(Y ) ≥−z (V ) − ∂ϕ(Y ). Note that ∂ϕ(Y ) ≤ δ |Y | · |V \Y | ≤ n(n − 1)δ/2. Therefore, we havex− (V ) ≥ z − (V ) − n(n − 1)δ/2, which together with Lemma 3.1 implies x− (V ) ≥f (X) − n2 δ.
Theorem 3.3. At the end of the last scaling phase, X is a minimizer of f .
Proof. Since δ < 1/n2 in the last scaling phase, Lemma 3.2 implies x− (V ) >f (X) − 1. Then it follows from the integrality of f that f (X) ≤ f (Y ) holds for any
Y ⊆V.
4. Complexity. This section is devoted to complexity analysis of the new scaling
algorithm.
Lemma 4.1. Each scaling phase performs Augment O(n2 ) times.
Proof. At the beginning of each scaling phase, the set X obtained in the previous
scaling phase satisﬁes z − (V ) ≥ f (X) − 2n2 δ by Lemma 3.2. For the ﬁrst scaling
phase, we have the same inequality by taking X = ∅. Note that z − (V ) ≤ z(X) ≤f (X) + n(n − 1)δ/2 throughout the algorithm. Thus each scaling phase increasesz − (V ) by at most 3n2 δ. Since each augmentation increases z − (V ) by δ, each scaling
phase performs at most 3n2 augmentations.
Lemma 4.2. Each scaling phase performs Relabel O(n2 ) times.
Proof. Each application of Relabel(v) increases d(v) by one. Since Relabel(v) isapplied only if d(v) < n, Relabel(v) is applied at most n times for each v ∈ V in ascaling phase. Thus the total number of relabels in a scaling phase is at most n2 .
Lemma 4.3. The number of indices in I is at most 2n.
Proof. A new index is added as a result of nonsaturating Multiple-Exchange. In
a nonsaturating Multiple-Exchange(i, u, v), the arc (q, r) that determines β satisﬁes
ϕ(q, r) ≤ 0 after the update of ϕ, and the vertex r in R becomes reachable from S
in Gϕ . This means the set W is enlarged. Thus there are at most n applicationsof nonsaturating Multiple-Exchange between augmentations. Hence the number ofindices added between augmentations is at most n. After each augmentation, thenumber of indices is reduced to at most n. Therefore, |I| ≤ 2n holds throughout thealgorithm.
In order to analyze the number of function evaluations in each scaling phase,
we now introduce the notion of reordering phase. A reordering phase consists ofconsecutive applications of Multiple-Exchange between those of Relabel or Reduce. By
Lemmas 4.1 and 4.2, each scaling phase performs O(n2 ) reordering phases.

838

SATORU IWATA

Lemma 4.4. There are O(n2 ) function evaluations in each reordering phase.
Proof. The number of function evaluations in Multiple-Exchange(i, u, v) is at mostthe number of vertices in the active interval for (i, u, v). In order to bound thetotal number of function evaluations in a reordering phase, suppose the procedure
Multiple-Exchange(i, u, v) marks each pair (i, w) for w in the active interval. We now
intend to claim that any pair (i, w) of i ∈ I and w ∈ V is marked at most once in areordering phase.
In a reordering phase, the algorithm does not change the labeling d nor does itdelete a vertex from W . Hence the minimum value of d in V \W is nondecreasing.
After execution of Multiple-Exchange(i, u, v), there will not be an active pair for i
until the minimum value of d in V \W becomes larger. Let Multiple-Exchange(i, s, t)
be the next application of Multiple-Exchange to the same index i ∈ I. Then we haved(t) > d(v) = d(u) + 1, which implies v ≺i u ≺i t ≺i s in the linear ordering Li before
Multiple-Exchange(i, u, v). Thus a pair (i, w) marked in Multiple-Exchange(i, u, v) will
not be marked again in the reordering phase.
Since |I| ≤ 2n by Lemma 4.3, there are at most 2n2 possible marks without
duplications. Therefore, the total number of function evaluations in a reordering
phase is O(n2 ).
In order to ﬁnd an active triple eﬃciently in Step (4-1), we keep track of possible
candidates of active triples. For each i ∈ I and - = 1, . . . , n − 1, let ui denote thelast vertex u in Li such that u ∈ W and d(u) = - − 1. Similarly, vi denotes the ﬁrst
vertex v in Li such that v ∈ V \W and d(v) = -. Then (i, ui , vi ) is an active triple if- = min{d(v) | v ∈ V \W } and vi ≺ ui . At the beginning of each reordering phase,
we scan the linear orderings to ﬁnd all those candidates in O(n2 ) time.
In the rest of the reordering phase, we update the candidates whenever a newvertex is added to W . Let w be the vertex that is added to W . For each i ∈ I, ifui  w with d(w) = - − 1, then we replace ui by w. If w = vi , then we ﬁnd the
new vi by scanning Li . Thus it takes O(n2 ) time to update the candidates when a
new vertex is added to W . Since at most n vertices are added to W , each reordering
phase requires O(n3 ) fundamental operations.
Theorem 4.5. The algorithm performs O(n4 log M ) function evaluations and5O(n log M ) arithmetic computations.
Proof. Consider the set U = {u | x(u) > 0} for the initial base x ∈ B(f ). Then we
have x− (V ) = x(V ) − x(U ) ≥ f (V ) − f (U ) ≥ −2M . Therefore, the initial value of δ
satisﬁes δ ≤ 2M/n2 . Each scaling phase cuts the value of δ in half, and the algorithm
terminates when δ < 1/n2 . Thus the algorithm consists of O(log M ) scaling phases.
Since each scaling phase performs O(n2 ) reordering phases, Lemma 4.4 implies
that the number of function evaluations in a scaling phase is O(n4 ). In addition, eachreordering phase requires O(n3 ) steps to keep track of active triples. By Lemma 4.1,
each scaling phase performs O(n2 ) calls of Reduce, which requires O(n3 ) arithmetic
computations. Thus each scaling phase consists of O(n4 ) function evaluations andO(n5 ) arithmetic computations. Therefore, the total running time bound is O((n4 γ +n5 ) log M ).
5. Discussions. A family D ⊆ 2V is called a distributive lattice (or a ringfamily) if X ∩ Y ∈ D and X ∪ Y ∈ D for any pair of X, Y ∈ D. A compact
representation of D is given by a directed graph as follows. Let D = (V, F ) be adirected graph with the arc set F . A subset Y ⊆ V is called an ideal of D if no arc
enters Y in D. Then the set of ideals of D forms a distributive lattice. Conversely,
any distributive lattice D ⊆ 2V with ∅, V ∈ D can be represented in this way due toSUBMODULAR FUNCTION MINIMIZATION

839

Birkhoﬀ’s representation theorem [1, Theorem 2.5]. Moreover, contracting strongly
connected components of D to single vertices, we may assume that the directed graph
D is acyclic.
For minimizing a submodular function f on D, we apply the scaling algorithm witha minor modiﬁcation. The modiﬁed version uses the directed graph Gϕ = (V, Aϕ ∪ F )instead of Gϕ = (V, Aϕ ). The initial linear ordering L0 must be consistent with D;i.e., v i u if (u, v) ∈ F . Then all the linear orderings that appear in the algorithm
will be consistent with D. This ensures that the set X obtained at the end of eachscaling phase belongs to D. Thus the modiﬁcation of our scaling algorithm ﬁnds a
minimizer of f in D.
Iwata, Fleischer, and Fujishige [13] also describe a strongly polynomial algorithm
that repeatedly applies their scaling algorithm with O(log n) scaling phases. The
number of iterations is O(n2 ). Replacing the scaling algorithm by the new one, we
obtain an improved strongly polynomial algorithm that runs in O((n6 γ + n7 ) log n)time.
A very recent paper [12] has shown that the strongly polynomial IFF algorithm
can be implemented by using only additions, subtractions, comparisons, and oracle
calls for function values. Similarly, the new strongly polynomial scaling algorithm canbe made fully combinatorial as follows.
The ﬁrst step towards a fully combinatorial implementation is to neglect Reduce.
This causes growth of the number of extreme bases for convex combination. However,the number is still bounded by a polynomial in n. Since the number of indices added
between augmentations is at most n, each scaling phase yields O(n3 ) new extreme
bases. Hence the number of extreme bases through the O(log n) scaling phases isO(n3 log n).
The next step is to choose an appropriate step length in Multiple-Exchange so thatthe coeﬃcients should be rational numbers with a common denominator bounded bya polynomial in n. Let θ denote the value of δ in the ﬁrst scaling phase. Then
κ = θ/δ is an integer. For each i ∈ I, we keep λi = µi /κ with an integer µi . We
then modify the deﬁnition of saturating Multiple-Exchange. Multiple-Exchange(i, u, v)
is now called saturating if λi ξ(q, r) ≤ ϕ(q, r) for every (q, r) ∈ Q × R. Otherwise,
it is called nonsaturating. In nonsaturating Multiple-Exchange(i, u, v), let ν be theminimum integer such that νξ(q, r) > ϕ(q, r)κ for some (q, r) ∈ Q × R. Such an
integer ν can be computed by binary search. Then the new coeﬃcients λk and λi aredetermined by µk := µi − ν and µi := ν. Thus the coeﬃcients are rational numbers
whose common denominator is κ, which is bounded by a polynomial in n throughthe O(log n) scaling phases. Then it is easy to implement this algorithm using only
additions, subtractions, comparisons, and oracle calls for the function values.
Finally, we discuss time complexity of the resulting fully combinatorial algorithm.
The algorithm performs O(n2 ) iterations of O(log n) scaling phases. Since it keeps
O(n3 log n) extreme bases, each scaling phase requires O(n6 log n) oracle calls forfunction evaluations and O(n7 log n) fundamental operations. Therefore, the total
running time is O((n8 γ + n9 ) log2 n). This improves the previous O(n9 γ log2 n) boundin [12] by essentially a linear factor in n.
In order to reduce this time complexity, McCormick [15] suggests a more eﬃcient
implementation for ﬁnding active triples. For each i ∈ I and - = 1, . . . , n, let σidenote the last vertex s in Li with d(s) = - − 1. Similarly, τi denotes the ﬁrst vertex
t in Li with d(t) = -. Then there is an active triple (i, u, v) with d(v) = - only if
τi ≺i σi . At the beginning of each reordering phase we scan the linear orderings to840

SATORU IWATA

ﬁnd all σi and τi in O(n|I|) time. Note that within the reordering phase, σi andτi are invariant until the algorithm performs Multiple-Exchange(i, u, v) with d(v) = -.
Once such a Multiple-Exchange is applied, there will be no active triples for the samei and - in the rest of the reordering phase.
For a pair of i and - with τi ≺ σi , we may restrict the search for active triples
to the interval between τi and σi in Li . Since these intervals are disjoint, thetotal number of fundamental operations required for ﬁnding active triples is O(n|I|)in each reordering phase. This reduces the number of fundamental operations in ascaling phase to O(n6 log n). Thus the resulting fully combinatorial algorithm runs inO(n8 γ log2 n) time.
Acknowledgments. The author is grateful to Lisa Fleischer, Satoru Fujishige,
Yasuko Matsui, Tom McCormick, and Kazuo Murota for stimulating conversations
and helpful comments on the manuscript. The idea of Multiple-Exchange was originally suggested by Satoru Fujishige as heuristics to improve the practical performanceof the IFF algorithm. Lisa Fleischer kindly pointed out an error in an earlier versionof this paper, and Tom McCormick generously allowed me to include his idea for
implementing the algorithm that leads to the bound of the fully combinatorial version. Finally, the author wishes to thank two anonymous referees for many helpful
suggestions to improve the presentation of this paper.
REFERENCES[1] M. Aigner, Combinatorial Theory, Springer-Verlag, Berlin, New York, 1979.
[2] W. H. Cunningham, Testing membership in matroid polyhedra, J. Combin. Theory Ser. B, 36(1984), pp. 161–188.
[3] W. H. Cunningham, On submodular function minimization, Combinatorica, 5 (1985), pp. 185–
192.
[4] J. Edmonds, Submodular functions, matroids, and certain polyhedra, in Combinatorial Structures and Their Applications, R. Guy, H. Hanani, N. Sauer, and J. Schönheim, eds., Gordon
and Breach, New York, 1970, pp. 69–87.
[5] L. Fleischer and S. Iwata, Improved algorithms for submodular function minimization andsubmodular ﬂow, in Proceedings of the 32nd ACM Symposium on Theory of Computing,Portland, OR, 2000, ACM, New York, 2000, pp. 107–116.
[6] L. Fleischer and S. Iwata, A push-relabel framework for submodular function minimization
and applications to parametric optimization, Discrete Appl. Math., to appear.
[7] L. Fleischer, S. Iwata, and S. T. McCormick, A faster capacity scaling algorithm for submodular ﬂow, Math. Program., 92 (2002), pp. 119–139.
[8] S. Fujishige, Submodular Functions and Optimization, North-Holland, Amsterdam, 1991.
[9] M. Grötschel, L. Lovász, and A. Schrijver, The ellipsoid method and its consequences incombinatorial optimization, Combinatorica, 1 (1981), pp. 169–197.
[10] M. Grötschel, L. Lovász, and A. Schrijver, Geometric Algorithms and Combinatorial
Optimization, Springer-Verlag, Berlin, 1988.
[11] S. Iwata, A capacity scaling algorithm for convex cost submodular ﬂows, Math. Programming,
76 (1997), pp. 299–308.
[12] S. Iwata, A fully combinatorial algorithm for submodular function minimization, J. Combin.
Theory Ser. B, 84 (2002), pp. 203–212.
[13] S. Iwata, L. Fleischer, and S. Fujishige, A combinatorial strongly polynomial algorithm forminimizing submodular functions, J. ACM, 48 (2001), pp. 761–777.
[14] L. Lovász, Submodular functions and convexity, in Mathematical Programming — The Stateof the Art, A. Bachem, M. Grötschel, and B. Korte, eds., Springer-Verlag, Berlin, 1983,
pp. 235–257.
[15] S. T. McCormick, Submodular function minimization, in Handbook of Discrete Optimization,
K. Aardal, G. Nemhauser, and R. Weismantel, eds., Elsevier, New York, to appear.
[16] A. Schrijver, A combinatorial algorithm minimizing submodular functions in strongly polynomial time, J. Combin. Theory Ser. B, 80 (2000), pp. 346–355.
[17] L. S. Shapley, Cores of convex games, Internat. J. Game Theory, 1 (1971), pp. 11–26.