Network flow algorithms

00 CD qt INETWORK FLOW ALGORITHMS Andrew V. Goldberg 0 Eva Tardos Robert E. Tarjan CS-TR-216-89 iL ] aZ mtr 1989 . NOV 291989 A ;- ;v.- 7- i r i A z "! e e::e ' 80f 161' Network Flow Algorithms Andrf w V. Goldbcrgi Department of Computer Sciuc(' Stanford Tniversity Stanford, CA 94305 Tardos2 LT'a Department of Mathematics. NI.I.T. Cambridge, NI A 02139 , ,:,th .rF,, i and j ., Computer Science Department Eotvos University Budapest Ely Robrrt E. Jajan3' )epartment of (oniputer Sciice Princeton Universi v Princeton, NJ 0851.1 AT&T Bell Laboratories Murray 1lii, - . ... NJ 0797-1 April 1989 'Research Development 2 Research 3 Research partially supported by NSF Presidential Young Investigator Grant CCR-8858097, 11 M l'aciilty Award, ONII Contract, N00014-88 K -0166, and AT&T,Bell Laboratories. partially supported by Air Force contract AFOSlI-86 0078. partially supported by the National Science Foundation, Grant No. 1C)13-86059;1, and the Office of Naval Research, Contract No. N0001,1-87-K-0467. . Contents Introduc-tion 3 1 Preliminaries 9 2 1.1 Flows and Residual Graphs .. .. .. .. ... 1.2 The Maximum Flow Problem .. .. .. .. ... 1.3 The Minimum-Cost Circulation Problem. .. .. .. .. ... 1.4 The Transshipmnt Problem. .. .. .. .. ... ... ... ... ... ... .... ... 2 1.5 The Generalized _71ow Problem. .. .. .. ... ... ... ... ... ... .... .. 14 1.6 The Restricted Pr Nbem .. ... .... 1.7 Decomposition T heorems. .. .. .. ... ... .. ... ... .... ... ... ... ... ... ... .... ... ... ... ... ... ... ... .... 9 .. .... 10 10 ... ... ... ... .. 16 ... ... ... ... ... .. 17 .... ... The Maximum Flow Problem 19 2.1 Introduction. .. .. .. .. ... 2.2 A Generic Algorithm .. .. .. .. ... 2.3 Efficient Implementation....... 2.4 Excess Scaling. 2.5 Use of Dynamic Trees. ... .. .. .. .... ... ... ... ... ... .... . ... ... .... ... ... . ... . .. .... ... ... ... ... ... .. .. .. ... ... ... ... .... ... ... ... .... ... .. .... 22 .. .. .. .. .. .. ... ... ... ... .... ... ... Introduction. .. .. .. .. ... 3.2 Approximate Optimality. .. .. .. .. ... ... 3.3 The Generalized Cost-Scaling Framework. .. .. .. ... 3.4 A Generic Refinement Algorithm. .. .. .. ... ... ... 3.5 Efficient Implementation. .. .. .. .. ... ... .... 3.6 Refinement Using Blocking Flows. .. .. .. .. ... ... ... ... ... ... 2R .... 30 33 .... ... ... .... ... 5 ..... 3 The Minimum-Cost Circulation Problem: Cost-Scaling 3.1 19 ... ... ... ... .... ... ... ... ... ... ... ... ... .. .... 33 34 ... ... .. 35 ... ... .. 35 ... .... ... .... ... 39 .. 41 2 4 5 CONTENTS Strongly Polynomial Algorithms Based on Cost-Scaling 45 4.1 Introduction . . . . . . . . . . . . . .. 4.2 Fitting Price Functions and Tight Error Parameters ....................... 45 4.3 Fixed Arcs ......... 47 4.4 The Strongly Polynomial Framework ................................. 49 4.5 Cycle-Canceling Algorithms ....................................... 50 . . . . . . . . . .. . . . . . . . . . . . . . ... ......................................... Capacity- Scaling Algorithms 45 55 5.1 Introduction ......... 5.2 The Edmonds-Karp Algorithm ....... 5.3 Strongly Polynomial Algorithms ....... ......................................... .............................. ............................. 6 The Generalized Flow Problem 55 56 58 63 6.1 Introduction ......... 6.2 Simple Combinatorial Algorithms .................................... 64 6.3 Polynomial-Time Combinatorial Algorithms ............................ 66 Bibliography ......................................... 63 72 Introduction Network flow problems are central problems in operations research, computer science, and engineering and they arise in many real world applications. Starting with early work in linear programming and spurred by the classic book of Ford and Fulkersoni'26, the study of such problems has led to continuing improvements in the efficiency of network flow algorithms. In spite of the long histcry of this study, many substantial results have been obtained within the last several years. In this survey we examine some of these recent developments and the ideas behind them. We discuss the classical network flow problems, the maimum flow problem and the mininuwcost circulation problem, and a less standard problem, the generalized flow problem, sometimes called the problem of flows with losses and gains. The survey contains six chapters in addition to this introduction. Chapter 1 develops the terminology needed to discuss network flow problems. Chapter 2 discusses the maximum flow problem. Chapters 3, 4, and 5 discuss different aspects of the minimum-cost circulation problem, and Chapter 6 discusses the generalized flow problem. In the remainder of this introduction, we mention some of the history of network flow research, comment on some of the results to be presented in detail in later sections, and mention some results not covered in this survey. We are interested in algorithms whose running time is small as a function of the size of the network and the numbers involved (e.g. capacities, costs, or gains). As a measure of the network size, we use n to denote the number of vertices and m to denote the number of arcs. As measures of the number sizes, we use U to denote the maximum arc capacity, C to denote the maximum arc cost, and B (in the case of the generalized flow problem) to denote the maximum numerator or denominator of the arc capacities and gains. In bounds using U, C, or B, we make the assumption that the capacities and costs are integers, and that the gains in the generalized flow problem are rational numbers. We are most interested in polynomial-time algorithms. We make the following distinctions. An algorithm is pseudopolynomial if its running time has a bound that is polynomial in n,m. and the appropriate subset of U,C, and B. An algorithm is polynomial if its running time has a bound that is polynomial in n,rn, and the appropriate subset of log U, logC, and log B 1 . An 'All logarithms in this paper without an explicit base are base two. 3 INTRODUCTION 4 algorithm is strot,gly ;olynomial if its running tim' has a bound that is polynomial in n and m2 . When comparing polynomial algorithms to strongly polynomial ones we shall use the similarity assumption that log U, log C and log B are O(log n) (32]. We shall also be interested in strongly polynomial algorithms, however. The network flow problems discussed in this survey are special cases of linear programming, and thus they can be solved by general-purpose linear programming algorithms. However, the combinatorial structure of these problems makes it possible to obtain more efficient algorithms. We shall not discuss in detail algorithms that are based on general linear programming methods. We should mention, however, that the first algorithm designed for network flow problems was the network simplex method of Dantzig [19]. It is a variant of the linear programming simplex method designed to take advantage of the combinatorial structure of network flow problems. Variants of the simplex method that avoid cycling give an exponential bound on the complexity of " the network flow problems. (Cunningham 1181 gives an elegant anti-cycling strategy for the network simplex method. based on graph-theoretic properties of the minimum-cost circulation problem). Recently, Goldfarb and Hao [51] have designed a variant of the primal network simplex method for the maximum flow problem that runs in strongly polynomial time (see Section 2.1). Orlin [80] designed a variant of the dual nctwork simplex method for the minimum-cost circulation problem that runs in strongly polynomial time (see Chapter 5). For a long time, the network simplex method has been the method of choice in practice, in particular for the mirtimuiu-.ost circulation problem (see e.g. [53]); for large instances of hard problems, the new scaling algorithms are probably better, however. The first pseudopolynomial algorithm for the maximum flow problem is the augmenting path algorithm of Ford and Fulkerson [25, 26]. Dinic [20] and Edmonds and Karp [21] independently obtained polynomial versions of the augmenting path algorithm. Since then, several more-efficient algorithms have been developed. Chapter 2 presents the push/relabel method, recently proposed by Goldberg [39] and Goldberg and Tarjan [44], along with some of its more efficient variants. The first pseudopolynomial algorithm for the minimum-cost circulation problem is the outof-kilter method, which was developed independently by Yakovleva [104], Minty [76], and Fulkerson [31]. The first polynomial algorithm for the minimum-cost circulation problem is due to Edmonds and Karp [21). To develop this algorithm Edmonds and Karp introduced the technique of scaling, which has proved to be a useful tool in the design and analysis of algorithms for a variety of combinatorial optimization problems. Chapter 3 and Section 5.2 are devoted to scaling algorithms for the minimum-cost circulation problem. The maximum flow algorithms of Dinic [20] and Edmonds and Karp [21] are strongly polynomial, but the minimum-cost circulation algorithm of Edmonds and Karp [21] is not. The first strongly polynomial algorithm for the ninimum-cost circulation problem was designed by Tar2 For a more formal definition of pclynomial and strongly polynomial algorithms, see [54]. INTRODUCTION 5 dos [95]. Chapter 4 and Section 5.3 are devoted to recent strongly polynomial algorithms for the minimum-cost circulation problem. The first augmenting path algorithms for the generalized flow problem were developed independently by Jewell [61, 62] and Onaga ['(8]. Many pseudopolynomial minimum-cost circulation algorithms have been adapted for the generalized flow problem (see [101] for a survey). The first polynomial-time algorithm for the generalized flow problem was the ellipsoid method [69]. Kapoor and Vaidya [64] have shown how to speed up Karmarkar [65] - or Renegar [88] - type interiorpoint algorithms on network flow problems by taking advantage of the special structure of the matrices used in the linear programming formulations of these problems. Vaidya's algorithm [102] is the fastest currently known algorithm for the generalized flow problem. The first polynomial algorithms for the generalized flow problem that are not based on general-purpose linear programming methods are due to Goldberg, Plotkin, and Tardos [42]. These algorithms are discussed in Chapter 6. The existence of a strongly polynomial algorithm for the generalized flow problem is an interesting open question. Important special cases of network flow problems that will not be covered in this survey are the bipartite matching problem and its weighted version, the assignment problem. These prohlems can be stated as maximum flow and minimum-cost circulation problems, respectively, on networks with unit capacities and a special structure (see e.g. [23, 96]). Some of the efficient algorithms for the more general problems have evolved from efficient algorithms developed earlier for these simpler problems. Kbnig's t71) proof of a good characterization of the maximum size of a matching in a bipartite graph gives an 0(nm)-time algorithm for finding a maximum matching. The Ford-Fulkerson maximum flow algorithm can be viewed as an extension of this algorithm. Hopcroft and Karp [57] gave an O(v' m) algorithm for the bipartite matching problem. Even and Tarjan observed [24] that Dinic's maximum flow algorithm, when applied to the bipartite matching problem, behaves similarly to the Hopcroft-Karp algorithm arid runs in O(,/f'if) time as well. A variation of the Goldberg-Tarjan maximum flow algorithm (which can be viewed as a generalization of Dinic's algorithm) can be easily shown to lead to the same bound [5, 83]. In spite of recent progress on related problems, the O(v/'ffm) bound has not been improved. The first algorithm for the assignment problem is the Hungarianmethod of Kuhn [72]. The outof-kilter algorithm is an extension of this algorithm to the minimum-cost circulation problem. The Hungarian method solves the assignment problem in O(n) shortest path computations. Edmonds and Karp [21] and Tomizawa [100) have observed that the dual variables can be maintained so that these shortest path computations are on graphs with non-negative arc costs. Combined with the shortest path algorithm of [29], this observation gives an O(n(m + n log n)) bound for the problem. Gabow (32] used scaling to obtain an O(n3/4m log C) algorithm for the problem. Extending ideas of the Hopcroft-Karp bipartite matching algorithm and those of the Goldberg-Tarjan minimum-cost 6 INTRODUCTION Problem Date Bipartite Matching 1973 Assignment Discoverer Running Time Hopcroft and Karp O(./'nm) [57] 1955 1987 Kuhn Gabow and Tarjan O(n(m + n log n)) O(/rT11log(aC)) [72] [34] Maximum Flow 1986 1988 Goldberg and Tarjan Ahuja, Orlin, and Tarjan O(nm log(n 2 /m)) O(nm log(-- Vl Minimum-Cost Circulation 1972 1987 1987 Edmonds and Karp Goldberg and Tarjan Orlin O(m log U(m + n log n)) O(nmlog(n2 /m) log(nC)) O(m log n(m + n log n)) 1988 Ahuja, Goldberg, Orlin. and Tarjan O(nm log log U log(nC)) 1989 Vaidya O(n 2 m' Generalized Flow References I I 5 [44] [4] + 2)) [21] [46] [82] [1] log(nB)) I [102] I Table 1: Fastest currently known algorithms for network flow problems. circulation algorithm (discussed in Section 3), Gabow and Tarjan obtained an O(N/'_rmlog(71C)) algorithm for the assignment problem. A more recent pseudopolynomial algorithm for the assignment problem is the auction algorithm of Bertsekas [9] (first published in [10]). This algorithm contains some of the elements of the push/relabel algorithms discussed in Sections 2 and 3. Some versions of the network simplex method have been shown to solve the assignment problem in polynomial ime. In particular, Orlin [81] shows that a natural version of the primal simplex method runs in polynomial time, and Balinski [6] gives a signature method that is a dual simplex algorithm for the assignment problem and runs in strongly polynomial time. For discussion of parallel algorithms for the bipartite matching and assignment problems, see [10, 33, 43, 67, 77]. In this survey we would like to highlight the main ideas involved in designing highly efficient algorithms for network flow problems, rather than to discuss in detail the fastest algorithms. However, for easy reference, we summarize the running times of the fastest currently known algorithms in Table 1. For each problem we list all of the bounds that are best for some values of the parameters, but we only list the first algorithm achieving the same bound. Some of the bounds stated in -. • •I INTRODUCTION 7 the table depend on the O(m + nlog n) implementation of Dijkstra's shortest path algorithm [29]. Chapter 1 Preliminaries In this chapter we define the problems addressed in this survey and review fundamental facts about these problems. These problems are the maximum flow problem, the minimum-cost circulation problem, the transshipment problem, and the generalized flow problem 1.1 Flows and Residual Graphs A network is a directed graph G = (V, E) with a non-negative capacity function u : E - R, .1 We assume that G has no multiple arcs, i.e., E C V x V. If there is an arc from a vertex c to a vertex w, this arc is unique by the assumption, and we will denote it by (v, w). This assumption is for notational convenience only. We also assume, without loss of generality, that the input graph G is symmetric: (v,w) E E (w,v) E E. A flow network is a network with two distinguished vertices, the source s and the sink t. . A pseudoflow is a function f : E -- R that satisfies the following constraints: f(v, w) <_u(v, u,) V(v, w) E E (capacity constraint), (1.1) f(r,w) = - f(',v) V(v,w) E E (1.2) (flr,- qntisymmetry constraint). Remark: To gain intuition, it is often useful to think only about the non-negative components of a pseudoflow (or of a generalized pseudoflow, defined below). The antisymmetry co-nstraints r-fl-ct the fact that a flow of value x from r to w can be thought of as a flow of value (-x) from u, to r. The negative flow values are introduced only for notational convenience. Note, for example, that one does not have to distinguish between lower and upper capacity bounds: the capacity of the arc (v, w) represents a lower bound on the flow value on the opposite arc. = RU oo. 10 CHAPTER 1. PRELIMINARIES Given a pseuo," ,w f, we define the excess function ef : V - R by ef(v) = Euv f(u, v), the net flow into v. We will say that a vertex v has excess if ef(v) is positive, and has deficit if it is negative,. For a vertex v, we define the conservation constraint by ef(v) = 0 (flow conservation constraint). (1.3) Given a pseudoflow f, the residual capacity function uf : E -- R is defined by uf(v. w) u(v,w) - f(v,w). The residual graph with respect to a pseudoflow f is given by G = (',Ef. where Ef = {(v,w) E Eluj(vw) > 0. 1.2 The Maximum Flow Problem To introduce the maximum flow problem, we need the following definitions in addition to the definitions of the previous section. Consider a flow network (G, u, s, t). A prcflou is a pseudoflow f such that the excess function is non-negative for all vertices other than s and t. A flou f on G is a pseudoflow satisfying the conservation constraints for all vertices except s and t. The ,'aluc IfI of a flow f is the net flow into the sink ef(t). A maximum flow is a flow of maximum value (also called an optimal flow). The maximum flow problem is that of finding a maxinmum flow in a given flow network. Given a flow f, we define an augmenting path to be a source-to-sink path in the residual graph. The foowing theorem, due to Ford and Fulkerson, gives an optimality criterion for maximum flows. Theorem 1.2.1 1.3 [26] A flow is optimal if and only if its residual graph contains no augmenting path. The Minimum-Cost Circulation Problem A circulation is a pseudoflow with zero excess at every vertex. A cost fu~iction is a real-valued function on arcs c : E - R. Without loss of generality, we assume that costs are antisymmetric: c(r,w) = -c(u,v) V(r,wi) E E (cost antisymmetry constraint). (1.4) The cost of a circulation f is given by E~f f=-,'O "0 (t,u-)EE:f(t,u>O The minimum-cost cirrulation problem is that of finding a minimum-cost (optimal) circulation in an input network (G,u,c). 1.3. TilE MINIMUM-COST CIRCULATION PROBLEM 11 We have assumed that the capacities are non-negative. This assumption is no loss of generality. Given a circulation problem with some negative capacities one can find a feasible circulation f using one invocation of a maximum flow algorithm. (See [86], problem 11(e), p. 215.) The non-negative residual capacities u(v, w) - f(v, w) then define an equivalent problem. Next we state two criteria for the optimality of a circulation. Define the cost of a cycle r to bc the sum of the costs of the arcs along the cycle and denote it by c(F). The following theorem states the first optimality criterion. Theorem 1.3.1 [15] A :irculation is optimal if and only if its residual graph contains no negative-cost cycle. To state the second criterion, we need the notions of a price function and a reduced cost function. A price function is a vertex labeling p : V , R. The reduced cost function with respect to a price function p is defined by cp(v, u,) = c(v, w)+p(v)-p(w). These notions, which originate in the theory of linear programming, axe crucial for many iniiimum-cost flow algorithms. As linear programming duld variables, vertex prices have a natural economic interpretation as the current market prices of the commodity. We can interpret the reduced cost p(t,u,) as the cost of buying a unit of commodity at v, transporting it to w, and then selling it. Due to the conservation constraints, the reduced costs define an equivalent problem. Theorem 1.3.2 [26] A circulation f is optimal for the minimum-cost circulation problem (G.u.c) if and only if it is optimal for the problem (G, u, cp) for every price function p. The second optimality criterion is as follows. Theorem 1.3.3 [26] A circulation f is optimal if and only if there is a price function p such that, for each arc (v, w), cp(v,w) < 0 , f(,, a,) = u(v,w) (complementary slackness constraint). A minimum-cost circulation problem may have no optimal solution. characterizes when this is the case. (1.5) The following theorem Theorem 1.3.4 There exists a minimum-cost circulation if and only if the input network contains no negative-cost cycle consisting of arcs with infinite capacity. Note that the maximum flow problem is a special case o the minimum-cost circulation problem (see e.g. [26]). To see this, consider a flow network and add a pair of new arcs (s, t) and (ts) with u(.*,t) = 0 and u(t,s) = oo. Define the costs of all original arcs to be zero, and define c(s,t) = -c(t,q) = 1. A minimum-cost circulation in the resulting network restricted to the original network is a maximum flow in the original network. 12 1.4 CHAPTER 1. PRELIMINARIES The Transshipment Problem IT this section we define the (uncapacitated) transshipment problem. Although this problem is equivalent to the minimum-cost circulation problem, some algorithms are more natural in the context of the trahsshipment problem. For , more extensive discussion of this problem and the closely related transportation problem, see [73]. In the transshipment problem, all arc capacties are either zero or infinity. In addition, the input to the problem contains a demand function d : V -, R such that -VEV d(v) = 0. For the transshipment problem, the notion of the oxcess at a vertex is defined as follows: f(w,v) - dxv). ef(v) = (1.6) wEV A , seudoflow f is feasible if the conservation constraints (1.3) hold for all vertices. The transshipment problem is that cr finding a minimum-cost feasible pseudoflow in an input network. In t;,e special case of integer demands, we shall use D to denote the maximum absolute value of a demand. Theorems analogous to Theorem 1.3.1, 1.3.3 and 1.3.2. hold for the transshipment problem. and the analog of Theorem 1.3.4 holds for a transshipment problem that has a feasible pseudoflow. We make two simplifying assumptions about the transshipinent problem. First, we assume that, for the zero flow, all residual arcs have non-negative cost. This assumption can be validated by first checking whether the condition of Theorem 1.3.4 is satisfied (using a shortest path computation). In the case of a transshipment problem, the residual graph of the zero flow consists of the arcs with infinite capacity. If this graph has no negative-cost cycles, then define p(v) to be the cost of a minimum-cost path from a fixed vertex s to the vertex v. The costs cp(v,w) define an equivalent problem that satisfies the assumption. Second, we assume that the residual graph of the zero flow is strongly connected. (There is a path between each pair of vertices.) This condition can be imposed by adding two new vertices x and y and adding an appropriately high-cost arc from x to y with u(z,y) = o and u(y,x) = 0, and arcs between the new vertices and every other vertex in both directions such that u(v,x) = u(y,v) = oc and u(x,v) = u(v,y) = 0 for every vertex v. If the original problem has a feasible solution, then in every minimum-cost solution of thp transformed problem all of the dummy arcs introduced to make the graph strongly connected have zero flow. Next we show that the transshipment problem is equivalent to the minimum-cost circulation problem (see e.g. [26]). Given an instance of the transshipment problem, we construct an equivalent instance of the minimum-cost circulation problem as follows. Add two new vertices, x and y, and arcs (y,z) and (z,y) with u(z,y) = 0 and u(y,x) = ,V:d(v)>Od(v). Define the cost of (y, X) to be small enough so that any simple cycle containing (y,x) has negative cost; for example, define c(y,x) = -c(x,y) = -nC. For every vertex v E V with d(v) > 0, add arcs (x,v) and (v.x) and define u(x,v) = d(v), u(v,x) = 0, and c(x,v) = c(v,x) = 0. For every vertex v E V with d(') < 0, 1.4. THE TRANSSHIPMENT PROBLEM 13 add arcs (v,y) and (y,v) and define u(v,y) = -d(v), u(y,v) = 0, and c(v,y) = c(y,r) = 0. Consider an optimal solution f to the minimum-cost circulation problem in the resulting network. The transshipment problem is feasible if and only if f saturates all new arcs, and in this case the restriction of f to the original arcs gives an optimal solution to the transshipment problem. Next we reduce the minimum-cost circulation problem to the transshipment problem. The reduction uses the technique of edge-splitting. Consider a minimum-cost circulation problem (G, u. c). First we make sure that for every edge {v,w} of the network, either u(v,w) or u(w,v) is infinite. For every edge {v,w} that does not satisfy the above condition, we introduce a new vertex x(t, ,), and replace the arcs (v,w) and (w,v) by the arcs (v, x(t,,w)), (x(vw),w), (w,x(t,)), and (x(v,), v). Define the costs and capacities of the new arcs as follows: u(v' XV.W)) = u(v,w). u(x(,W,, ') = 00, u(w, X(,,.)) u(w, V), u(I(V,), v) = oc, c(v,x(VW)) c(v, W), c( ',) = 0, c(v,X(V,u)) = 0, c(X(V,W.), v) c(w, v). This defines an equivalent minimum-cost circulation problem in which the capacity of every edge is unbounded in at least one direction. To complete the transformation, we need to force the additional restriction that every finite capacity is zero. Initialize the demand function d to be zero on all vertices. Let {t w} be an edge such that u(w,v) = oo and u(v,w) is finite but not zero. Replacing u(v,u,) by zero, and d(v) and d(w) by d(v) + u(v, w) and d(w) - u(v, u,), respectively, defines an equivalent problem with u(v, U,) = O. Remark: This transformation increases the number of vertices: from a network with n vertices and m edges, we obtain a network with up to n + m vertices and 2m edges. However, the resulting network has a special structure. It is important for the analyses in Chapter 5 to notice that the newly introduced vertices can be eliminated for shortest path computations. Consequently, this blowup in the number of vertices will not affect the running time of shortest path computations in the residual graph. The capacitated transshipment problem generalizes both the minimum-cost circulation problem and the (uncapacitated) transshipment problem discussed above. It is the same as the uncapacitated transshipment problem without the assumption that all arc capacities are infinite or zero. The reductions above can be extended to show that the capacitated transshipment problem is equivalent to the minimum-cost circulation problem. Whereas the previous simpler formulations are better suited for designing algorithms, the more general form can be useful in applications. 14 1.5 CHAPTER 1. PRELIMINARIES The Generalized Flow Problem In this section we define generalized pseudoflows, generalized flows, and the generalized flow problem, and compare these concepts with their counterparts discussed in the previous sections. For alternative formulations of the problem, see e.g. [42, 73]. (This problem is also known as the problem of flows with losses and gains.) The generalized flow problem is given by a network with a distinguished vertex, the source s, and a gain function -y : E -+ R + 2 on the arcs. We assume (without loss of generality) that the gain function is antisymmetric: -(v, w) = I V(v,w) E E (gain antisymmetry constraints). (1.7) In the case of ordinary flows, if f(v,w) units of flow are shipped from v to w, f(v,w) units arrive at w. In the case of generalized flows, if g(v,w) units of flow are shipped from v to u,, 7(v, w)g(v, w) units arrive at u. A generalized pseudoflow is a function g : E -* R that satisfies the capacity constraints (1.1) and the generalized antisymmetry constraints: g(v,w) = -y(w, v)g(w.v) V(v,w) E E (generalized antisymmetry constraints). (1.8) The gain of a path (cycle) is the product of the gains of the arcs on the path (cycle). For a given generalized pseudoflow g, the residual capacity and the residual graph G. = (V, E_) are defined in the same way as for pseudoflows. The excess eg(v) of a generalized pseudoflow g at a vertex v is defined by eg(v)=- E g(v,w). (1.9) (v,w)E E We will say that a vertex v has excess if e2 (v) is positive and deficit if it is negative. A generalized flow is a generalized pseudoflow that satisfies the conservation constraints (1.3) at all vertices except at the source. The value of a generalized pseudoflow g is the excess e9 (s). The generalized flow problem is that of finding a generalized flow of maximum value in a given network. In our discussion of this problem, we shall assume that the capacities are integers, and that each gain is a rational number given as the ratio of two integers. Let B denote the largest integer used ;n the specification of the capacities and gains. A flow-generating cycle is a cycle whose gain is greater than 1, and a flow-absorbing cycle is a cycle whose gain is less than 1. Observe that if one unit of flow leaves a vertex v and travels along a flow-generating cycle, more than one unit of flow arrives at v. Thus we can augment the flow 'R+ denotes the set of positive reals. 1.5. THE GENERALIZED FLOW PROBLEM 15 along this cycle; that is, we can increase the excess at any vertex of the cycle while preserving the excesses at the other vertices, by increasing the flow along the arcs in the cycle (and correspondingly decreasing the flow on the opposite arcs to preserve the generalized antisymmetry constraints). The generalized flow problem has the following interpretation in financial analysis. The commodity being moved is money, nodes correspond to different currencies and securities, arcs correspond to possible transactions, and gain factors represent the prices or the exchange rates (see Section 6.1). From the investor's point of view, a residual flow-generating cycle is an opportunity to make a profit. It is possible to take advantage of this opportunity, however, only if there is a way to transfer the profit to the investor's bank account (the source vertex). This motivates the following definition. A generalized augmenting path (GAP) is a residual flow-generating cycle and a (possibly trivial) residual path from a vertex on the cycle to the source. Given a generalized flow and a GAP in the residual graph, we can augment the flow along the GAP, increasing the value of the current flow. The role of GAP's in the generalized flow problem is similar to the role of negative-cost cycles in the minimum-cost circulation problem - both can be used to augment the flow and improve the value of the current solution. Onaga [79 proved that the non-existence of GAP's in the residual graph characterizes the optimality of a generalized flow. Theorem 1.5.1 [79] A generalized flow is optimal if and only if its residual graph contains no GAP. Using the linear programming dual of the problem, it is possible to give an alternate optimality criterion, similar to the criterion in Theorem 1.3.3 for the minimum-cost circulation problem. A price function p is a labeling of the vertices by real numbers, such that p(s) = 1. As in the case of the minimum-cost circulation problem, vertex prices can be interpreted as market prices of the commodity at vertices. If a unit of flow is purchased at v and shipped to w, then Y(v,uw) units arrive at w. Buying the unit at v costs p(v), and selling y(v,w) units at w returns p(uw)-(v,u,). Thus the reduced cost of (v, w) is defined as c(v, w) = p(v) - p(w) (v, w). Linear programming duality theory provides the following optimality criterion. Theorem 1.5.2 [61] A generalized flow is optimal if and only if there exists a price function p such that the complementary slackness conditions (1.5) hold for each arc (v, w) E E. One generalization of the problem, namely the generalized flow problem with costs, is worth mentioning here. As the name suggests, in this problem each arc (v,w) has a cost c(v, u,) in addition to its gain. The goal is to find a generalized flow that maximizes C,(s) - c(g), where c(g) = F(v,,):g( w)>0 c(v, w)g(v, w). The introduction of costs enriches the combinatorial structure of the problem and allows the modeling of more complex problems, in particular economic processes. For 16 CHAPTER I. PRELIMINARIES example, a positive cost flow-generating cycle with a path leading to a negative cost flow-absorbing cycle may be used as a producer-consumer model. The generalized flow problem with costs has not been studied as much as the other problems discussed in this survey, and the only polynomial-time algorithms known for this problem are based on general linear programming methods [64, 102]. 1.6 The Restricted Problem Next we introduce a special variant of the generalized flow problem, and show that this variant is equivalent to the original problem. Consider for a moment the following variation of the generalized flow problem: given a flow network with a source s E V and a sink t E V, find a generalized pseudoflow with maxizr.:7n -,cess at t and zero excess at each vertex other than s and t. Onaga [78] suggested the stud) of the special case of this problem in which the residual graph of the zero flow has no flow-generating cycles. We shall consider the corresponding special case of the generalized flow problem in which the residual graph of the zero flow has the property that all flow-generating cycles pass through the source. (If there are no flow-generating cycles, the zero flow is the optimal.) We shall also assume that the residual graph of the zero flow is strongly connected. A generalized flow network in which the residual graph of the zero flow is strongly connected and which has no flow-generating cycles not passing through the source is called a restricted network. The restricted problem is the generalized flow problem on a restricted network. The restricted problem has a simpler combinatorial structure that leads to simpler algorithms. Moreover, it turns out that the restricted problem is equivalent to the generalized flow problem. All of the algorithms that we review solve the restricted problem. In the rest of this section we shall review basic properties of the restricted problem, and we outline the reduction. In Chapter t3 we will use the term "generalized flow problem" to mean the restricted problem. One of the nice facts about the restricted problem is that the optimality condition given by Theorem 1.5.1 simplifies in this case, and becomes very similar to Theorem 1.3.1. This characterization, which is also due to Onaga [78], can be deduced from Theorem 1.5.1 with the use of the following lemma. The lemma can be proved by induction on the number of arcs with positive flow. Lemma 1.6.1 Let g be a generalized pseudoflow in a restricted network. If the excess at every vertex is non-negative, then for every vertex v there exists a path from v to s in the residual graph G. Theorem 1.6.2 [79] A generalized flow g in a restricted network is optimal if aaid only if the residual graph of g contains no flow-generating cycles. Note that the condition in the theorem can be formulated equivalently as "if and only if the residual graph of g has no negative-cost cycles, where the cost function is given by c = - log I." 1.7. DECOMPOSITION THEOREMS 17 Theorem 1.6.3 [42] The generalized flow problem can he reduced to the restricted problem in O(nw) time. Proofsketch: Given an instance of the generalized flow problem, reduce it to an instance in which all the flow-generating cycles pass through s, as follows. Let h be the generalized pseudoflow that saturates all arcs of gain greater than 1 and is zero on all arcs of gain 1. Define a new generalized flow network containing each arc (v, w) of the original network, with new capacity u(v, w) - h(v, w). For each vertex v E V such that eh(v) > 0, add an arc (s, v) of capacity u(s, v) = eh(v)/a, and with gain -(s,v) = a, where a = B'. Also add a reverse arc (v,s) with u(v,s) = 0. For each vertex v E V with eh(v) < 0, add an arc (v,s) of capacity u(v,s) = eh(t') and with gain 'y(v,s) = a; also add a reverse arc (s,v) with u(v,s) = 0. Let 4 be an optimal generalized flow in the new network. Consider the generalized pseudoflow 4 + h restricted to the arcs of the original network. Using the Decomposition Theorem 1.7.3 below, one can show that the value of this generalized pseudoflow is equal to the value of the optimal generalized flow, and an optimal generalized flow can be constructed from this generalized pseudoflow in O(nm) time. Intuitively, the new arcs ensure that vertices having excesses with respect to h are supplied with an adequate amount of "almost free" flow, and vertices having deficits with respect to h can deliver the corresponding amount of flow very efficiently to the source. Having eliminated flow-generating cycles not containing the source, we can impose the strong connectivity condition by deleting all of the graph except the strongly connected component con1 taining the source. This does not affect the value of an optimum solution by Theorem 1.5.1. Note that this transformation increases the size of the numbers in the problem description. However, the polynomial-time algorithms described later depend only on T (< B"), the maximum product of gains along a simple path, and L (< Bm), the product of denominators of arc gains: these parameters do not increase significantly. 1.7 Decomposition Theorems A useful property of flows and generalized flows is the fact that they can be decomposed into a small number of "primitive elements". These elements depend on the problem under consideration. In this section we review decomposition theorems for circulations, pseudoflows, and generalized pseudoflows. In the case of circulations, the primitive elements of the decomposition are flows along simple cycles. Theorem 1.7.1 [26] For every circulation f, there exists a collection of k < m circulations gi... gk such that for every i, the graph induced by the set of arcs on which gi is positive is a simple cycle consisting of arcs (v, w) with f(v, w) > 0, and 18 CHAPTER 1. PRELIMINARIES f(vw) = g (1.10) Such a decomposition can be found in O(nm) time. For pseudoflows (or flows), the primitive elements of the decomposition are flows along simple cycles and flows along simple paths. Theorem 1.7.2 [26] For every pseudoflow f, there exists a collection of k < m pseudoflows g,.... , 9k, such that for every i, the graph induced by the set of arcs on which gi is positive is either a simple cycle or a simple path from a vertex with deficit to a vertex with excess; it consists of arcs (v, w) with f(v,w) > 0, and (1.10) is satisfied. Such a decomposition can be found in O(nr) time. The five primitive elements of the decomposition for generalized pseudoflows g are defined as follows. The set of arcs on which an element of the decomposition of g is positive is a subset of the arcs on which g is positive. Let T(g) denote the set of vertices with excesses and let S(g) denote the set of vertices with deficits. The five types of primitive elements are classified according to the graph induced by the set of arcs with positive flow. Type I: A path from a ,'ertex in S(g) to a vertex in T(g). Such a flow creates a deficit and an excess at the two ends of the path. Type II: A flow-generating cycle and a path leading from this cycle to a vertex in T(g). Such a flow creates excess at the end of the path. (If the path ends at the source, then this corresponds to a GAP.) Type III: A flow-absorbing cycle and a path from a vertex in S(g) to this cycle. Such a flow creates a deficit at the end of the path. Type IV: A cycle with unit gain. Such a flow does not create any' excesses or deficits. Type V A pair of cycles connected by a path, where one of the cycles generates flow and the other one absorbs it. Such a flow does not create any excesses or deficits. Theorem 1.7.3 [42, 52] For a generalized pseudoflow g there exist k < m primitive pseudoflows ga,... ,gk such that, for each i, gi is of one of the five types described above, and (1.10) is satisfied. Such a decomposition can be found in O(nm) time. Remark: In all three of these theorems the time to find the claimed decomposition can be reduced to O(m log n) by using the dynamic tree data structure that will be discussed in Section 2.3. Chapter 2 The Maximum Flow Problem 2.1 Introduction The maximum flow problem has been studied for over thirty years. The classical methods for solving this problem are the Ford-Fulkerson augmenting path method [25, 26], the closely related blocking flow mr.xthod of Dinic [20], and an appropriate variant of the ntwork sivnlex method (see e.g. (601). A push-relabel method for solving maximum flow problems has been recently proposed by Goldberg [39] and fully dcveloped by Goldberg and Tarjan [44, 48]. Recently, parallel and distributed algorithms for the maximum flow problem have been studied as well. Many maximum flow algorithm se scaling techniques and data structures to achieve better running time bounds. The scaling teciniques used by maximum flow algorithms are capacity scaling, introduced by Edmonds and Karp [21] in the context of the minimum-cost circulation problem, and the closely related excess scaling of Ahuja and Orlin [3]. The dynamic tree data structure of Sleator and Tarjan [93, 94] makes it possible to speed up many maximum flow algorithms. The technical part of this survey deals with the push-relabel method, which gives better sequential and parallel complexity bounds than previously known methods, and seems to outperform them in practice. In addition to describing the basic method, we show how to use excess scaling and dynamic trees to obtain faster implementations of the method. In this section we discuss the pr'Vious approaches and their relationship to the push-relabel method. The augmenting path algorithm of Ford and Fulkerson is based on the fact that a flow is maximum if and only if there is no augmenting path. The algorithm repeatedly finds an augmenting path and augments along it, until no augmenting path exists. Although this simple generic method need not terminate if the network capacities are reals, and it can run in exponential time if the capacities are integers represented in binary [26], it can be made -fficient by restricting the choice of augmenting paths. Edmonds and Karp [21] have shown that two strategies for selecting t!.e next augmenting path give polynomial time bounds. The first such strategy is to select, at each iteration, a shortest augmenting path, where the length of a path is the number of arcs on the path. 19 20 CHAPTER 2. THE MAXIMUM FLOW PROBLEM The second strategy is to select, at each iteration, a fattest path, where the fatness of a path is the minimum of the residual capacities of the arcs on the path. Independently, Dinic [201 suggested finding augmenting paths in phases, handling all augm,-nting paths of a given shortest length in one phase by finding a blocking flow in a layered network. The shortest augmenting path length increases from phase to phase, so that the blocking flow method terminates in n - 1 phases. Dinic's discovery motivated a number of algorithms for the blocking flow problem. Dinic proposed an O(nm) blocking flow algorithm. Soon thereafter, Karzanov [68] proposed an O(n 2 ) algorithm, which achieves the best known bound for dense graphs. Karzanov also introduced the concept of a preflow; this concept is used by many maximum flow algorithms. Simpler algorithms achieving an O(n 2 ) bound are described in [74, 97]. A sequence of algorithms achieving better and better running times on non-dense graphs has been proposed [17, 32, 35. 36, 46, 91, 94]. The algorithm of [32] uses capacity scaling; the algorithmns of [94] and [46] use dynamic trees. The fastest currently-known sequential algorithm for the blocking flow problem. due to Goldberg and Tarjan, runs in O(mlog(n2 /m)) time [46]. The push-relabel method [39, 40, 44, 46] replaces the layered network of the blocking flow method by a distance labeling of the vertices, which gives an estimate of the distance to the sink :_ .he iehiduai graph. The algorithm maintains a preflow and a distance labeling, and uses two simple operations, pushing and relabeling, to update the preflow and the labeling, repeating them until a maximum flow is found. The pushing operation is implicit in Karzanov's algorithm. The construction of a layered network at each iteration of Dinic's method can be viewed as a sequence of relabeling operations. Unlike the blocking flow method, the push-relabel methoe solves the maximum flow problem directly, without reducing the problem to a sequence of blocking flow subproblems. As a result, this method is more general and flexible than the blocking flow method. and leads to improved sequential and parallel algorithms. Ahuja and Orlin [3] suggested one way of recasting algorithms based on Dinic's approach into the push-relabel framework. Goldberg's FIFO' algorithm [39] runs in O(n 3 ) time. Goldberg and Tarjan [44, 48] showed how to use the dynamic tree data structure to improve the running time of this algorithm to O(nm log(n 2 /m)). They also gave a largest-labelvariant of the algorithm, for which they obtained the same time bounds - O(n 3 ) without dynamic trees and O(nmlog(n2 /m)) with such trees. Cheriyan and Maheshwari [16] showed that (without dynamic trees) the FIFO algorithm runs in f?(n 3 ) time, whereas the largest-label algorithm runs in 0(n 2v&7 ) time. Ahuja and Orlin [3] described an O(nm + n 2 log U)-time version of the push-relabel method based on excess scaling. Ahuja, Orlin, and Tarjan [4] gave a modification of the Ahuja-Orlin algorithm that runs in O(nm + n 2 V4AgU) time without the use of dynamic trees and in O(nm log(! v 4o + 2)) time with them. The primal network simplex method (see e.g. [60]) is another classical method for solving the maximum flow problem. Cunningham [18] gave a simple pivoting rule that avoids cycling due to 'FIFO is an abbreviation of first-in, first-out. 2.1. INTRODUCTION 21 degeneracy. Of course, if a simplex algorithm does not cycle, it must terminate in exponential time. Until recently, no polynomial time pivoting rule was known for the primal network simplex method. Goldfarb and Hao [51] have given such a rule. The resulting algorithm does 0(nm) pivots, which correspond to 0(nm) flow augmentations; it runs in 0(n 2 m) time. Interestingly, the blocking flow method also makes O(nm) flow augmentations. Unlike the blocking flow method, the GoldfarbHao algorithm does not necessarily augment along shortest augmenting paths. In their analysis, Goldfarb and Hao use a variant of distance labeling and a variant of the relabeling operation mentioned above. Dynamic trees can be used to obtain an 0(nm log n)-time implementation of their algorithm [41]. The best currently known sequential bounds for the maximal flow problem are O(n7 log(n 2 /M)) and 0(nm log(-2--T +- 2)). Note that, although the running times of the known algorithms come very close to an O(mn) bound, the existence of a maximum flow algorithm that meets this bound remains open. With the increa.sinf, inrpst in parallel rompiting, parallel and distributed agorithms for the maximum flow problem have received a great deal of attention. The first parallel algorithm for the problem, due to Shiloach and Vishkin [92], runs in 0(n 2 log n) time on an n-processor PRAM [27] and uses 0(n) memory per processor. In a synchronous distributed model of computation, this algorithm runs in 0(n 2 ) time using 0(n 3 ) messages and 0(n) memory per processor. The algorithm of Goldberg [39, 40, 48] uses less memory than that of Shiloach and Vishkin: 0(m) memory for the PRAM implementation and O(A) memory per processor for the distributed implementation (where A is the processor degree in the network). The time, processor, and message bounds of this algorithm are the same as those of the Shiloach-Vishkin algorithm. Ahuja and Orlin [3] developed a PRAM implementation of their excess-scaling algorithm. The resource bounds are [rm/n] processors, 0(m) memory, and 0(n 2 lognlog U) time. Cheriyan and Maheshwari [16] proposed a synchronous distributed implementation of the largest-label algorithm that runs in 0(2) time using O(A) memory per processor and O(n 2 -,/-m) messages. For a long time, the primal network simplex method was the method of choice in practice. A study of Goldfarb and Grigoriadis [50] suggested that the algorithm of Dinic [20] performs better than the network simplex method and better than the later algorithms based on the blocking flow method. Recent studies of Ahuja and Orlin (personal communication) and Grigoriadis (personal communication) show superiority of various versions of the push-relabel method to Dinic's algorithm. An experimental study of Goldberg [40] shows that a substantial speedup can be achieved by implementing the FIFO algorithm on a highly parallel computer. More efficient algorithms have been developed for the special case of planar networks. Ford and Fulkerson [25] have observed that the the maximum flow problem on a planar network is related to a shcrtczt path problem on the planar dual of the network. The algorithms in [8, 25, 98, 55, 56, 22 CHAPTER 2. THE MAXIMUM FLOW PROBLEM procedure generic (V, E, u); [initialization) V(v, w) E E do begin f(v, w) - 0; if v = s then f(s, w) -- u(s, w); if w = s then f(v, s) - -u(s, v); end; Vw E V do begin ef(W) - E.(,.,,)E f (V,W); if w = s then d(w) = n else d(u,) = 0; end; [loop] while 3 an active vertex do select an update operation and apply it; return(f); end. Figure 2.1: The generic maximum flow algorithm. 59, 63, 75, 87] make clever use of this observation. 2.2 A Generic Algorithm In this section we describe the generic push-relabel algorithm [40, 44, 48]. First, however, we need the following definition. For a given preflow f, a distance labeling is a function d from the vertices to the non-negative integers such that d(t) = 0, d(s) = n, and d(v) <_d(w) + 1 for all residual arcs (v, w). The intuition behind this definition is as follows. Define a distance graph G as follows. Add an arc (s,t) to Gf. Define the length of all residual arcs to be equal to one and the length of the arc (s,t) to be n. Then d is a "locally consistent" estimate on the distance to the sink in the distance graph. (In fact, it is easy to show that d is a lower bound on the distance to the sink.) We denote by dG; (v, w) the distance from vertex v to vertex w in the distance graph. The generic algorithm maintains a preflow f and a distance labeling d for f, and updates f and d using push and relabel operations. To describe these operations, we need the following definitions. We say that a vertex v is active if v V {s, t} and ef(v) > 0. Note that a preflow f is a flow if and only if there are no active vertices. An arc (v,w) is admissible if (v,w) E ES and d(v) = d(w) + 1. The algorithm begins with the preflow f that is equal to the arc capacity on each arc leaving the source and zero on all arcs not incident to the source, and with some initial labehng d. The algorithm then repetitively performs, in any order, the update operations,push and rclabl,described in Figure 2.2. When there are no active vertices, the algorithm terminates. A summary of the algorithm appears in Figure 2.1. 2.2. A GENERIC ALGORITHM 23 push(v, w). Applicability: v is active and (v, w) is admissible. Action: send 6 E (0, min(e! (v), u! (v, w))] units of flow from v to w. relabel(v). Applicability: either s or t is reachable from v in G! and Vw E V uj(r,w) = 0 or d(w) > d(r). Action: replace d(v) by rrn'n(vw)EE{d(u)} + 1. Figure 2.2: The update operations. The pushing operation updates the preflow, and the relab~ling operation updates the distance labeling. Except for the excess scaling algorithm, all algorithms discussed in this section push the maximum possible amount b when doing a push. The update operations modify the preflow f and the labeling d. A push from v to w increases f(v,w) and ef(w) by up to 6 = min{cf(v), uff(v, u,)}, and decreases f(w, v) and (ff(v) by the same amount. The push is saturatingif us(v,w) = 0 after the push and nonsaturating otherwise. A relabelingof v sets the label of v equal to the largest value allowed by the valid labeling constraints. There is one part of the algorithm we have not yet specified: the choice of an initial labeling d. The simplest choice is d(s) = n and d(v) = 0 for v C V - {s). A more accurate choice (indeed. the most accurate possible choice) is d(v) = dG (V,t) for v, G V, where f is the initial preflow. The latter labeling can be computed in 0(m) time using backwards breadth-first searches from the sillk and from the source in the residual graph. The resource bounds we shall derive for the algorithm arc correct for any valid initial labeling. To simplify the proofs, we assume that the a.]gorithm starts with the simple labeling. In practice, it is preferable to start with the most accurate values of the distance labels, and to update the distance labels periodically by using backward breadtih-first search. Remark: By giving priority to relabeling operations, it is possible to maintain the following invariant: Just before a push, d gives the exact distance to t in the distance graph. Furthermore. it is possible to implement the relabeling operations so that the total work done to maintain the distance labels is O(nm) (see e.g. [48]). Since the running time bounds derived in this section are Q(,,n). one can assume that the relabeling is done in this way. In practice, however, maintaining exact distances is expensive; a better solution is to maintain a valid distance labehng and periodically update it to the exact labehng. Next we turn our attention to the correctness and termination of the algorithm. Our proof of correctness is based on Theorem 1.2.1. The following lemma is important in the analysis of the algorithm. Lemma 2.2.1 If f is a preflow and v is a vertex with positive excess, then the source s is reachable from v in the residual graph Gf. Using this lemma and induction on the number of update operations, it can be shown that 24 CHAPTER 2. THE MAXIMUM FLOW PROBLEM one of the two update operations must be applicable to an active vertex, and that the operations maintain a valid distance labeling and preflow. Theorem 2.2.2 [48] Suppose that the algorithm terminates. Then the preflow f is a maximum flow. Proof: When the algorithm terminates, all vertices in V - {s, t} must have zero excess, because there are no active vertices. Therefore f must be a flow. We show that if f is a preflow and d is a valid labeling for f, then the sink t is not reachable from the source s in the residual graph Gf. Then Theorem 1.2.1 implies that the algorithm terminates with a maximum flow. Assume by way of contradiction that there is an augmenting path s = vo, vi, . .. . . = t.Then 1 < n and (vi, v+,1 ) E Ef for 0 < i < 1. Since d is a valid labeling, we have d(?,,) <_d(1',+) + I for 0 < i < 1. Therefore, we have d(:') < d(t) + I < n, since d(t)= 0, which contradicts d(s)= n. I The key to the running time analysis of the algorithm is the following lemma. which shows that distance labels cannot increase too much. Lemma 2.2.3 At any time during the execution of the algorithm, for any vertex v C ', d(r) < 2n - 1. Proof: The lemma is trivial for v = s and v = t. Suppose v G V - {.t}. Since the algorithm changes vertex labels by on!y by means of the rclabeling operation. it is enough to prove the lemma for a vertex v such that s or t is reachable from v in Gf. Thus there is a simple path from 1'to s or t in GC. Let r = t0 , ..... 7,1be such a path. The length I of the path is at most i - 1.Since d is a valid labe!ing and (r,, v,+) E F n or 0.we have d(c) = d(v) < , we have d(t, ) < d(z',+ + 1.Therefore. since d(v'j) iseit her d(ri) + I< n + (n - 1) = 2n - 1. I Lemma 2.2.3 limits the number of relabeling operations, and allows us to amortize the work done by the algorithm over increases in vertex labels. The next two lemmas bound the number of relabelings and the number of saturating pushes. Lemma 2.2.4 The number of relabeling operations is at most 2n - 1 per vertex and at most (2n 1)(n - 2) < 2n 2 - overall. Lemma 2.2.5 The number of saturating pushes is at most nm. Proof: For an arc (, u,) E E, consider the saturating pushes from v to u,. After one such push. uf(r. it,) = 0, and another such push cannot occur until d(u) increases by at least 2,a push from it to v occurs, and d(v) increases by at least 2. If we charge each saturating push from r to u'except the first to the preceding label increase of v, we obtain an upper bound of n on the number of such pushes. U 2.3. EFFICIENTIMPLEMENTATION 25 The most interesting part of the analysis is obtaining a bound on the number of nonsaturating pushes. For this we use amortized analysis and in particular the potential function technique (see e.g. [98]). Lemma 2.2.6 The number of nonsaturating pushing operations is at most 2n 2 m. Proof: We define the potential 4t of the current preflow f and labeling d by the formula 4 -- {vlv is active) d(v). We have 0 < 4) < 2n 2 = by Lemma 2.2.3. Each nonsaturating push, say from a vertex v to a vertex w, decreases 4t by at least one, since d(w) = d(v) - 1 and the push makes 1' inactive. It follows that the total number of nonsaturating pushes over the entire algorithm is at most the sum of the increases in 4D during the course of the algorithm, since 4 = 0 both at the beginning and at the end of the computation. Increasing the label of a vertex v by an amount k increases 4) by k. The total of such increases over the algorithm is at most 2n 2 . A saturating push can increase 4) by at most 2n - 2. The total of such increases over the entire algorithi is at most (2n - 2)nm. Summing gives a bound of at most 2n nonsaturating pushes. 2 + (2n - 2)nm < 2n 2 m on the number of I Theorem 2.2.7 [48] The generic algorithm terminates after 0(n 2 m) update operations. Proof: Immediate from Lemmas 2.2.4, 2.2.5, and 2.2.6. 1 The running time of the generic algorithm depends upon the order in which update operatiMs are applied and on implementation details. In the next sections we explore these issues. First we give a simple implementation of the generic algorithm in which the time required for the nonsaturatinv pushes dominates the overall running time. Sections 2.3 and 2.A specify orders of the update, operations that decrease the number of nonsaturating pushes and permit 0(n 3 ) and 0(0in ;.2 log U)-time implementations. Section 2.5 explores an orthogonal approach. It shows how to u-e sophisticated data structures to reduce the time per nonsaturating push rather than the number of such pushes. 2.3 Efficient Implementation Our first step toward an efficient implementation is a way of combining the update operations locally. We need some data structures to represent the network and the preflow. We call an E E an edge of G. We associate the three values u( ',,r). with each edge {v,w}. Each vertex ' has a list of the iticidlet unordered pair {r,u'} such that (v,wu) u(u,t'), and f(ru)(= -f(u,,v)) {rw}, in fixed but arbitrary order. Thus each edge {r,u'} appears in exactly two list. the one for v and the one for w. Each vertex v hal a current (dgr {0%w }, which is the current candidate edges 26 CHAPTER 2. THE MAXIMUM FLOW PROBLEM discharge(v). Applicability: v is active. Action: let {v, w} be the current edge of v; time-to-relabel - false; repeat if (v, w) is admissible then push(v, w) else if {v, t,) is not the last edge on the edge list of v then replace {v, w} as the current edge of v by the next edge on the list else begin make the first edge on the edge li," of v the current edge; time-to-relabel - true; end; until ef(v) = 0 or time-to-relabel; if time-to-relabel then relabel(v); Figure 2.3: The discharge operation. for a pushing operation from v. Initially, the current edge of v is the first edge on the edge list of v. The main loop of the implementation consists of repeating the discharge operation described in Figure 2.3 until there are no active vertices. (We shall discuss the maintenance of active vertices later.) The discharge operation is applicable to an active vertex v. This operation iteratively attempts to push the excess at v through the current edge {v,w} of v if a pushing operation is applicable to this edge. If not, the operation replaces {v, w} as the current edge of v by the next edge on the edge list of v; or, if {v, w} is the last edge on this list, it makes the first edge on the list the current one and relabels v. The operation stops when the excess at v is reduced to zero or v is relabeled. The following lemma shows that discharge does relabeling correctly; the proof of the lemma is straightforward. Lemma 2.3.1 The discharge operation does a relabeing only when the relabeling operation is applicable. Lemma 2.3.2 The version of the generic push/relabel algorithm based on discharging runs in O(nm) time plus the total time needed to do the nonsaturating pushes and to maintain the set of active vertices. Any representation of the set of active vertices that allows insertion, deletion, and access to some active vertex in 0(l) time results in an O(n 2 m) running time for the discharge-based algorithm. by Lemmas 2.2.6 ana 2.3.2. (Pushes can be implemented in 0(1) time per push.) By processing active vertices in a more restricted order, we obtain improved performance. Two natural orders were suggested in [44, 481. One, the FIFO algorithm, is to maintain the set of active vertices as a queue, always selecting for discharging the front vertex on the queue, I adding newly 2.3. EFFICIENT IMPLEMENTATION 27 procedure process-vertez remove a vertex v from Bb; old-label - d(v); discharge(v); add each vertex w made active by the discharge to Bd(w); if d(v) $ old-label then begin b - d(v); add v to Bb; end else if Bb = end. then b -- b - 1; Figure 2.4: The process-vertex procedure. active vertices to the rear of the queue. The other, the largest-labcl algorithm, is to always select for discharging a vertex with the largest label. The FIFO algorithm runs in 0(n and the largest-label algorithm runs in O(n 2lv/) 3 time [16]. We shall derive an 0(n ) time [4-1. 4] 3 ) time bound for both algorithms, after first describing in a little more detail how to implement largest-label selection. The implementation maintains an array of sets Bi, 0 < i < 2n - 1, and an index b into the array. Set Bi consists of all active vertices with label i, represented as a doubly-linked list, so that insertion and deletion take 0(1) time. The index b is the largest label of an active vertex. During the initialization, when the arcs going out of the source are saturated, the resulting active vertices are placed in B 0 , and b is set to 0. At each iteration, the algorithm removes a vertex from Bb, processes it using the discharge operation, and updates b. The algorithm terminates when b becomes negative, i.e., when there are no active vertices. This processing of vertices, which implements the while loop of the generic algorithm, is described in Figure 2.4. To understand why the process-vertex procedure correctly maintains b, note that dischargc(v) either relabels v or gets rid of all excess at v, but not both. In the former case, v is the active vertex with the largest distance label, so b must be increased to d(v). In the latter case, the excess at v has been moved to vertices with distance labels of b - 1, so if Bb is empty, then b must be decreased by one. The total time spent updating b during the course of the algorithm is 0(n 2 ). The bottleneck in both the FIFO method and the largest-label method is the number of nonsaturating pushes. We shall obtain an O(n a ) bound on the number of such pushes by dividing the computation into phases, defined somewhat differently for each method. For the FIFO method. phase 1 consists of the discharge operations applied to the vertices added to the queue by the initialization of f; phase i+ 1, for i < 1, consists of the discharge operations applied to the vertices added to the queue during phase i. For the largest-label method, a phase consists of a maximad interval of time during which b remains constant. Lemma 2.3.3 The number of phases during the running of either the FIFO or the largest-label algo- CHAPTER 2. THE MAXIMUM FLOW PROBLEM 28 rithm is at most 4n 2 . Proof: Define the potential 4P of the current f and d by 4) = max{vlv is &tive}d(v), with the maximum taken to be zero if there are no active vertices. (In the case of the largest-label algorithm, 4) = b except on termination.) There can be only 2n 2 phases that do one or more relabelings. A phase that does no relabeling decreases 4t by at least one. The initial and final values of 4D are zero. Thus the number of phases that do no relabeling is at most the sum of the increases in 4) during the computation. The only increases in 4D are due to label increases; an increase of a label by k can cause 4) to increase by up to k. Thus the sum of the increases in 4) over the computation is at most 2n 2 , and so is the number of phases that do no relabeling. I Theorem 2.3.4 [48] Both the FIFO and the largest-label algorithm run in 0(n3) time. Proof: For either algorithm, there is at most one nonsaturating push per vertex per phase. Thus by Lemma 2.3.3 the total number of nonsaturating pushes is 0(n 3 ), as is the running time by Lemma 2.3.2. 1 Cheriyan and Maheshwari [16], by means of an elegant balancing argument, were able to improve the bound on the number of nonsaturating pushes in the largest-label algorithm to 0(n 2VM-). giving the following result: Theorem 2.3.5 [16] The largest-label algorithm runs in O(n 2 v'iM) time. 2.4 Excess Scaling A different approach to active vertex selection leads to running time U of the largest capacity as well as on the graph size. This approach, by Ahuja and Orlin [3] and developed further by Ahuja, Orlin, and in detail a slight revision of the original excess-scaling algorithm, bounds dependent on the size excess scaling, was introduced Tarjan [4]. We shall describe which has a running time of 0(nm + n 2 logU). For the termination of the excess-scaling method, all arc capacities must be integral; hence we assume throughout this section that this is the case. The method preserves integrality of the flow throughout the computation. It depends on a parameter A that is an upper bound on the maximum excess of an active vertex. Initially A - 2 PogU1. The algorithm proceeds in phases; after each phase, A is halved. When A < 1, all active vertex excesses must be zero, and the algorithm terminates. Thus the algorithm terminates after at most log 2 U + 1 phases. To maintain the invariant that no active vertex excess exceeds A, the algorithm does not always push the maximum possible amount when doing a pushing operation. Specifically, when pushing from a vertex r to a 2.4. EXCESS SCALING 29 vertex w, the algorithm moves an amount of flow 6 given by 6 = min{ef(v), uj(v,w), A -ej(u,)) if w 5 t,b = min{ef(v), uj(v,w)} if w = t. That is, 6 is the maximum amount that can be pushed while maintaining the invariant. The algorithm consists of initialization of the preflow, the distance labels, and A, followed by a sequence of process-vertex operations of the kind described in Section 2.3. Vertex selection for process-vertex operations is done by the large excess, smallest label rule: process an active vertex v with ef(v) > A/2; among all candidate vertices, choose one of smallest label. When every active vertex v has e1 (v) < A/2, A is halved and a new phase begins; when there are no active vertices, the algorithm stops. Since the excess-scaling algorithm is a version of the generic process-vertex-based algorithm described in Section 2.3, Lemma 2.3.2 applies. The following lemma bounds the number of nonsaturating pushes: Lemma 2.4.1 The number of nonsaturating pushes during the excess-scaling algorithm is O(n 2 log U). Proof: We define a potential 4) by 4 = vjl, . acti,,} e1 (v)d(v)/A. Since ef(v)/A < I for any active vertex, 0 < 4D < 2n 2 . Every pushing operation decreases 4D. Consider a nonsatiirating push from a vertex v to a vertex w. The large excess, smallest label rule guarantees that before the push et(v) > A/2 and either ej(u)) < A/2 or w = t. Thus the push moves at least A/2 units of flow, and hence decreases 4Dby at least 1/2. The initial and final values of 4) are zero, so the total number of nonsaturating pushes is at most twice the sum of the increases in (Dover the course of the algorithm. Increasing the label of a vertex by k can increase 4) by at most k. Thus relabelings account for a total increase in 4D of at most 2n 2 . A change in phase also increases 4), by a factor of two, or at most n 2 . Hence the sum of increases in 4) is at most 2n 2 + n 2 (log 2 U + 1). and therefore the number of nonsaturating pushes is at most 4n 2 + 2n 2 (log 2 U + 1). 1 The large excess, smallest label rule can be implemented by storing the active vertices with excess exceeding A/2 in an array of sets, as was previously described for the largest-label rule. In the former case, the index b indicates the nonempty set of smallest index. Since b decreases only 2 when a push occurs, and then by at most 1, the total time spent updating b is 0 (nm + n log U). From Lemmas 2.3.2 and 2.4.1 we obtain the following result: Theorem 2.4.2 [31 The excess-scaling algorithm runs in 0 (nm + n 2 log U) time. A more complicated version of excess scaling, devised by Ahuja, Orlin, and Tarjan [4]. has a running time of 0 (nm + n 2 Vg U). This algorithm uses a hybrid vertex selection rule that combines a stack-based mechanism with the "wave" approach of Tarjan [97]. CHAPTER 2. THE MAXIMUM FLOW PROBLEM 30 make-tree(v): Make vertex v into a one-vertex dynamic tree. Vertex v must be in no other tree. find-root(v): Find and return the root of the tree containing vertex v. find-value(v): Find and return the value of the tree arc connecting v to its parent. If v is a tree root, return infinity. find-min(v): Find and return the ancestor w of v such that the tree arc connecting w to its parent has minimum value along the path from v to find.root(v). In case of a tie, choose the vertex w closest to the tree root. If v is a tree root, return v. change-value(v,z): Add real number x to the value of every arc along the path from v to find-rool(v). link(v, w, z): Combine the trees containing v and w by making w the parent of v and giving the new tree arc joining v and w the value x. This operation does nothing if v and u; are in the same tree or if v is not a tree root. cut(v): Break the tree containing v into two trees by deleting the arc from v to its parent. This operation does nothing if v is a tree root. Figure 2.5: The dynamic tree operations. 2.5 Use of Dynamic Trees The two previous sections discussed ways of reducing the number of nonsaturating pushes by restricting the order of the update operations. An orthogonal approach is to reduce the time per nonsaturating push rather than the number of such pushes. The idea is to perform a succession of pushes along a single path in one operation, using a sophisticated data structure to make this possible. Observe that immediately after a nonsaturating push along an arc (v, w), (v, w) is still admissible, and we know its residual capacity. The dynamic tree data structure of Sleator and Tarjan [93, 941 provides an efficient way to maintain information about such arcs. We shall describe a dynamic tree version of the generic algorithm that has a running time of 0 (nm log n). The dynamic tree data structure allows the maintenance of a collection of vertex-disjoint rooted trees, each arc of which has an associated real value. We regard a tree arc as directed toward the root, i.e., from child to parent. We denote the parent of a vertex v by parrnt(v), and adopt the convention that every vertex is both an ancestor and a descendant of itself. The data structure supports the seven operations described in Figure 2.5. A sequence of I tree operations on trees of maximum size (number of vertices) k takes 0(1 log k) time. In the dynamic tree algorithm, the arcs of the dynamic trees are a subset of the admissible arcs and every activw vertex is a tree root. The value of an arc (v,w) is its residual capacity. Initially each vertex is made into a one-vertex dynamic tree using a make-tree operation. The heart of the dynamic tree implementation is the tree-push operation described in Figure 2.6. This operation is applicable to an admissible arc (v, w) such that v is active. The operation adds (v, w) to the forest of dynamic trees, pushes as much flow as possible along the tree from v to the root of the tree cc-itaining w, and deletes from the forest each arc that is saturated by this flow 2.5. USE OF DYNAMIC TREES 31 Tree-push (v, w) Applicability: v is active and (v, w) is admissible. Action: link (v,w, u1 (v, w)); parent (v) --w; 6 ,- min {e (v), find-value(find-min (v)); change-value (v,-6); while v $ find-root(v) and find-value(find-min(v)) = 0 do begin z *- find-min(v); cut (z); f(z, parent(z)) - u(z, parent(z)); f(parent(z), z) ---u(z, parent(z)); end. Figure 2.6: The tree-push operation. change. The dynamic tree algorithm is just the generic algorithm with the push operation replaced by the tree-push operation, with the initialization modified to make each vertex into a dynamic tree. and with a postprocessing step that extracts the correct flow value for each arc remaiing in a dynamic tree. (This postprocessing takes one find-value operation per vertex.) It is straightforward to show that the dynamic tree implementation is correct, by observing that it maintains the invariant that between tree-pushing operations every active vertex is a dynamic tree root and every dynami, tree arc is admissible. The following lemma bounds the number of tree-push operations: Lemma 2.5.1 The number of the tree-push operations done by the dynamic tree implementation is O(nm). Proof: Each tree-push operation either saturates an arc (thus doing a saturating push) or decreases the number of active vertices by one. The number of active vertices can increase, by at most one, as a result of a tree-push operation that saturates an arc. By Lemma 2.2.5 there are at most nm saturating pushes. An upper bound of 2nm + n on the number of tree-push operations follows. since initially there are at most n active vertices. I If we implement the dynamic tree algorithm using the discharge operation of Section 2.3 with push replaced by tree-push, we obtain the followirg result (assuming active vertex selection takes 0(1) time): Theorem 2.5.2 The discharge-basedimplementation of the dynamic tree algorithm runs in O(nm log n) time. Proof: Each tree-pushingoperation does 0(1) dynamic tree operations plus 0(1) per arc saturated. 32 CHAPTER 2. THE MAXIMUM FLOW PROBLEM The theorem follows from Lemma 2.2.5, Lemma 2.3.2, and Lemma 2.5.1, since the maximum size of any dynamic tree is O(n). I The dynamic tree data structure can be used in combination with the FIFO, largest-label, or excess-scaling method. The resulting time bounds are O(nmlog(n 2/M)) for the FIFO and largest-label methods [48] and O(nmlog(aVx/ U ) + 2)) for the fastest version of the excess- scaling method [4]. In each case, the dynamic tree method must be modified so that the trees are not too large, and the analysis needed to obtain the claimed bound is rather complicated. Chapter 3 The Minimum-Cost Circulation Problem: Cost-Scaling 3.1 Introduction Most polynomial-time algorithms for the minimum-cost circulation problem use the idea of scaling. This idea was introduced by Edmonds and Karp [21], who used it to develop the first polynomialtime algorithm for the problem. Scaling algorithms work by obtaining a sequence of feasible or almost-feasible solutions that are closer and closer to the optimum. The Edmonds-Karp algorithm scales capacities. R6ck [89] was the first to propose an algorithm that scales costs. Later, Bland and Jensen [13] proposed a somewhat different cost-scaling algorithm, which is closer to the generalized cost scaling method discussed in this chapter. The cost-scaling approach works by solving a sequence of problems, P0 , P1 ..... Pk, on the network with original capacities but approximate costs. The cost function c i for P, is obtained by taking the i most significant bits of the original cost function c. The first problem P0 has zero costs, and therefore the zero circulation is optimal. An optimal solution to problem P- 1 can be used to obtain an optimal solution to problem Pi in at most n maximum flow computations [13. 89]. Note that for k = [log 2 C[, ck = c. Thus the algorithm terminates in O(nlogC) maximum flow computations. Goldberg and Tarjan [40, 45, 47] proposed a generalized cost-scaling approach. The idea of this method (which is described in detail below) is to view the maximum amount of violation of the complementary slackness conditions as an error parameter, and to improve this parameter by a factor of two at each iteration. Initially the error is at most C, and if the error is less than 1/n, then the current solution is optimal. Thus the generalized cost-scaling method terminates in O(log(nC)) iterations. The computation done at each iteration is similar to a maximum flow computation. The traditional cost-scaling method of Rdck also improves the error from iteration to iteration, but it does so indirectly, by increasing the precision of the current costs and solving the 33 34 CHAPTER 3. THE MINIMUM-COST CIRCULATION PROBLEM: COST-SCALING resulting problem exactly. Keeping the original costs, as does the generalized cost-scaling approach, makes it possible to reduce the number of iterations required and to obtain strongly-polynomial running time bounds. Chapter 4 discusses a strongly-polynomial version of the generalized costscaling method. For further discussion of generalized versus traditional cost-scaling, see [46]. Time bounds for the cost-scaling algorithms mentioned above are as follows. The algorithms of Rhck [89] and Bland and Jensen [13] run in O(n log(C)M(n, m, U)) tine, where M(n, m, U) is the time required to compute a maximum flow on a network with n vertices, m arcs, and maximum arc capacity U. As we have seen in Chapter 2, M = O(nm log min{n 2 /m,.av4u + 21). The fastest known implementation of the generalized cost-scaling method runs in O(nmlog(n2 /M)log(nC)) time [45]. It is possible to combine cost scaling with capacity scaling. The first algorithm that combines the two scaling techniques is due to Gabow and Tarjan [34]. A different algorithm was proposed by Ahuja et al. [1]. The latter algorithm runs in O(nmloglogUiog(nC)) time, which makes it the fastest known algorithm for the problem under the similarity assumption. 3.2 Approximate Optimality A key notion is that of approximate optimality, obtained by relaxing the complementary slackness constraints in Theorem 1.3.3. For a constant c > 0, a pseudoflow f is said to be c-optimal with respect to a price function p if, for every arc (v, w), we have f(v,w) < u(v,w) . cp(v,w) >_-E (c-optimality constraint). (3.1) A pseudoflow f is c-optimal if f is c-optimal with respect to some price function p. An important property of c-optimality is that if the arc costs are integers and Cis small enough. any c-optimal circulation is minimum-cost. The following theorem, of Bertsekas [11] captures this fact.. Theorem 3.2.1 [11] If all costs are integers and c < 1/n, then an c-optimal circulation f is minimum- cost. The r-optimality constraints were first published by Tardos [95] in a paper describing the first strongly polynomial algorithm for the minimum-cost circulation problem. Bertsekas [11] proposed a pseudopolynomial algorithm based upon Theorem 3.2.1; his algorithm makes use of a fixed C< 1/n. Goldberg and Tarjan [40, 46, 47] devised a successive approximation scheme that produces a sequence of circulations that are c-optimal for smaller and smaller values of c; when C is small enough, the scheme terminates with an optimal circulation. We discuss this scheme below. 3.3. THE GENERALIZED COST-SCALING FRAMEWORK 35 procedure min-cost(V, E, u, c), [initialization] f -C; Vv, p(v) -- 0; V(v, w) EE, f(v, w) .--0; [loop] while c > 1/n do (f, f,p) -- refine(c, f,p); return(f); end. Figure 3.1: The generalized cost-scaling method. 3.3 The Generalized Cost-Scaling Framework Throughout the rest of this chapter, we assume that all arc costs are integral. We give here a highlevel description of the generalized cost-scaling method (see Figure 3.1). The algorithm maintains an error parameter c, a circulation f and a price function p, such that f is c-optimal with respect to p. The algorithm starts with c = C (or alternatively c = 2 P'9 2 C1), with p(v) = 0 for all v E V, and with the zero circulation. Any circulation is C-optimal. The main loop of the algorithm repeatedly reduces the error parameter c. When c < 1/n, the current circulation is minimum-cost, and the algorithm terminates. The task of the subroutine refine is to reduce the error in the optimality of the current circulation. The input to refine is an error parameter c, a circulation f, and a price function p such that f is E-optimal with respect to p. The output from refine is a reduced error parameter c, a new circulation f, and a new price function p such that f is c-optimal with respect to p. The implementations of refine described in this survey reduce the error parameter c by a factor of two. The correctness of the algorithm is immediate from Theorem 3.2.1, assuming that refine is correct. The number of iterations of refine is O(log(nC)). This gives us the following theorem: Theorem 3.3.1 [47] A minimum-cost circulation can be computed in the time required for O(log(nC)) iterations of refine, if refine reduces c by a factor of at least two. 3.4 A Generic Refinement Algorithm In this section we describe an implementation of refine that is a common generalization of the generic maximum flow algorithm of Chapter 2.2 and the auction algorithm for the assignment problem [9] (first published in (101). We call this the generic implemcntation. This implementation. proposed by Goldberg and Tarjat [471, is essentially the same as the main loop of the minimumcost circulation algorithm of Bertsekas [11], which is also a common generalization of the maximum 36 CHAPTER 3. THE MINIMUM-COST CIRCULATION PROBLEM: COST-SCALING procedure refine(, [initialization. f, p); C/2; V(v,w) E E do if cp(v,w) < 0 then begin f(v,w) ,- u(v,w);f(w,v) -- -u(v,w); end; [loop] f4- while 3 an update operation that applies do select such an operation and apply it; return(e, f,p); end. Figure 3.2: The generic refine subroutine. push(v, wv). Applicability: v is active and (t, u)) is admissible. Action: send 6 = min(ej(v),uj(v,tw)) units of flow from v to w. relabel(v). Applicability: v is active and Vw E V (u,(v, u,) = 0 or cp(v, u,) Action: replace p(v) by max(Vu,)E {p(w) - c(v, w) - }. 0). Figure 3.3: The update operations for the generic refinement algorithm. Compare with Figure 2.2. flow and assignment algorithms. The ideas behind the auction algorithm can be used to give an alternative interpretation to the results of (40, 471 in terms of relaxation methods; see (12]. As we have mentioned in Section 3.3, the effect of refine is to reduce c by a factor of two while maintaining the f-optimality of the current flow f with respect to the current price function p. The generic refine subroutine is described on Figure 3.2. It begins by halving C and saturating every arc with negative reduced cost. This conve-ts the circulation f into an c-optimal pseudoflow (indeed, into a 0-optimal pseudoflow). Then the subroutine converts the (-optimal pseudufliw into an c-optimal circulation by applying a sequence of the update operations push and relabel, each of which preserves e-optimality. The inner loop of the generic algorithm consists of repeatedly applying the two update operations, described in Figure 3.3, in any order, until no such operation applies. To define these operations, we need to redefine admissible arcs in the context of the minimum-cost circulation problem. Given a pseudoflow f and a price function p, we say that an arc (v, u,) is admissiblc if (v, w) is a residual arc with negative reduced cost. A push operation applies to an admissible arc (v,w) such that vertex v is active. It consists of pushing 6 = min{ej(v),uj(v,w)} units of flow from v to w,, thereby decreasing e (v,) and f(w. r) by 6 and increasing ef(w) and f(v, w) by 6. The push is saturating if uf(v, w) = 0 after the push and nonsaturatingotherwise. 3.4. A GENERIC REFINEMENT ALGORITHM 37 A relabd operation applies to an active vertex v that has no exiting admissible arcs. It consists of decreasing p(v) to the smallest value allowed by the c-optimality constraints, namely maX(vw)EEI {-c(v, w) + p(w) - C}. If an c-optimal pseudoflow f is not a circulation, then either a pushing or a relabeling operation is applicable. It is easy to show that any pushing operation preserves (-optimality. The next lemma gives two important properties of the relabeling operation. Lemma 3.4.1 Suppose f is an c-optimal pseudoflow with respect to a price function p and a vertex t, is relabeled. Then the price of v decreases by at least c and the pseudoflow f is -optimal with respect to the new price function p'. Proof: Before the relabeling, cp(v, w) _>0 for all (v,w) E Ef, i.e., p(r) > p(w) - cp(t,,) for all (v,w) E Ef. Thus p'(v) = max(v,w)EEf {p(w) - c(v,w) - } _< p(v) - f. To verify c-optimality, observe that the only residual arcs whose reduced costs are affected by the relabeling are those of the form (v, u,) or (u, v). Any arc of the form (U, v) has its reduced cost increased by the relabeling, preserving its c-optimality constraint. Consider a residual arc (v. ?r). By the definition of p', p'(v) >_p(w) - c(v,w) - c. Thus c, (v,w) = c(v,U!) + p'(v) - p(w) which means that (v, w) satisfies its c-optimality istraint. -, C 3 Since the update operations preserve c-optimality, and since some update operation applies if f is not a circulation, it follows that if refine terminates and returns (c,f,p), then f is a circulation which is c-optimal with respect to p. Thus refine is correct. Next we analyze the number of update operations that can take place during an execution of refine. We begin with a definition. The admissible graph is the graph GA = (I',EA) such that EA is the set of admissible arcs. As refine executes, the admissible graph changes. An important invariant is that the admissible graph remains acyclic. Lemma 3.4.2 Immediately after a relabeling is applied to a vertex v, no admissible arcs enter v. Proof: Let (u,v) be a residual arc. Before the relabeling, c,(u,v) > -( by (-optimality. By Lemma 3.4.1, the relabeling decreases p(v), and hence increases cp(u,v), by at least c. Thus cp(u,v) > 0 after the relabeling. I Corollary 3.4.3 Throughout the running of refine, the admissible graph is acyclic. Proof: Initially the admissible graph contains no arcs and is thus acyclic. Pushes obviously preserve acyclicity. Lemma 3.4.2 implies that reiabelings also preserve acyclicity. Next we derive a crucial lemma, which generalizes Lemma 2.2.1. I 38 CHAPTER 3. THE MINIMUM-COST CIRCULATION PROBLEM: COST-SCALING Lemma 3.4.4 Let f be a pseudoflow and f' a circulation. For any vertex v with e1 (v) > 0, there is a vertex w with ef(w) < 0 and a sequence of distinct vertices v = v 0 ,v 1 ,.. . ,vj..i,v1 = w such that (vi, vi+1 ) E Ef and (vi+j,v,) E E 1 , for 0 < i < 1. Proof: Let v be a vertex with ej(v) > 0. Define G+ = (VE+), where E+ = {(z,y)jf'(x,y) > f(z,y)}, and define G_ = (V,E_), where E_ = {(x,y)lf(x,y) > f'(x,y)). Then E+ C Ef, since (x,y) E E+ implies f(x,y) < f'(x,y) < u(x,y). Similarly E_ C Ef'. Furthermore (x,y) E E+ if and only if (y, x) E E_ by antisymmetry. Thus to prove the lemma it suffices to show the existence in G+ of a simple path v = vo, v, ... , vi with e 1 (t,) < 0. Let S be the set of all vertices reachable from v in G+ and let S = V - S. (Set S may be empty.) For every vertex pair (x,y) E S x 3, f(x,y) >_f'(x,y), for otherwise y E S. We have 0 f'(x, y) since f' is a circulation < Z(x,Y)E(Sx-)nE f(x, y) holds term-by-term = Z(Zi.)E(Sx3)nE by antisymmetrv by definition of S by antisymmetry. = Z(,y)E(S×S)nE f(X, y) + Z( ,Y)E(SxS)nE f(X, y) = Yz, = - )E(SxV)nE f(x, y) E Se(X) But v E S. Since ef(v) > 0, some vertex u e S must have ef(w) < 0. 3 Using Lemma 3.4.4 we can bound the amount by which a vertex price can decrtase during an invocation of rrfinc. Lemma 3.4.5 The price of any vertex v decreases by at most 3n( during an execution of refine. Proof: Let f2, aind P2, be the circulation and price functions on entry to rfinc. Suppose a relabeling cauoes the price of a vertex v to decrease. Let f be the pseudoflow and p the price function just after the relabeling. Then f.(v) > 0. Let v = v 0 ,r 1 ,.. .,t' = w with ef(w) < 0 be the vertex sequencp satisfying Lemma 3.4.4 for f and f' = f2,. The (-optimabty of f implies 1-1 - 1h < E cp(V,, V,+) i=0 1-1 = p(v) - p(w) + E c(',, t'i+i). i=0 (3.2) The 2c-optimality of f2, implies i-I - 21c < i-I -p 2 ,(r,+,,Vi) = P2,(u') - p 2 ((') + c(Vt+ 1 , v,). i=0 =0 (3.3) 3.5. EFFICIENTIMPLEMENTATION 39 But Ec(vi,vi+i) = - E,-c(vi+,vi) by cost antisymmetry. Furthermore. p(w) = p2,(u') since during refine, the initialization step is the only one that makes the excess of some vertices negative, and a vertex with negative excess has the same price as long as its excess remains negative. Adding inequalities (3.2) and (3.3) and rearranging terms thus gives p(v) >_p2,(v) - 31c > P2,(V) - 3nc. I Now we count update operations. The following lemmas are analogous to Lemmas 2.2.4. 2.2.5 and 2.2.6. Lemma 3.4.6 The number of relabelings during an execution of refine is at most 3n per vertex and 3n(n - 1) in total. Lemma 3.4.7 The number of saturating pushes during an execution of refine is at most 3nm. Proof: For an arc (v,wu), consider the saturating pushes along this arc. Before the first such push can occur, vertex v must be relabeled. After such a push occurs, v must be relabelcd before another such push can occur. But v is relabeled at most 3n times. Summing over all arcs gives the desired bound. I Lemma 3.4.8 The number of nonsaturating pushes during one execution of refine is at most 3n2( 71 + 71 Proof: For each vertex r, let iD(v) be the number of vertices reachable from v in the current admissible graph GA. Let it = 0 if there are no active vertices, and let 4) = {4(v)vc is active} otherwise. Throughout the running of rrfine, 4) > 0. Initially -t < n, since GA has no arcs. Consider the effect on 4) of update operations. A nonsaturating push decreases 4) by at least one. since GA is always acyclic by Corollary 3.4.3. A saturating push can increase 4Dby at most n. since at most one inactive vertex becomes active. If a vertex v is relabeled, 4Dalso can increase by at most n, since 4 (w) for w $ v can only decrease by Lemma 3.4.2. The total number of nonsaturating pushes is thus bounded by the initial value of 4D plus the total increase in 4) throughout the algorithm. i.e., by n + 3n 2 (n - 1) + 3n 2 ni < 3n 2 (m + n). 3.5 | Efficient Implementation As in the case of the generic maximum flow algorithm, we can obtain an especially efficient version of refine by choosing the order of the update operations careftly. Local ordering of the basic operations is achieved using the data structures and the dis.rhar( operation of Section 2.3. The discharge operation is the same as the one in Figure 2.3, but uses the 40 CHAPTER 3. TIlE MINIMUM-COST CIRCULATION PROBLEM: COST-SCALING procedure first-active; let L be a list of all vertices; let v be the first vertex in L; while 3 an active vertex do begin if v is active then begin discharge(v); if the discharge has relabeled v then move v to the front of L; end; else replace v by the vertex after v on L, end; end. Figure 3.4: The first-active method. minimum-cost circulation versions of the update operations (see Figure 3.3). As in the maximum flow case, it is easy to show that discharge applies the relabeling operation correctly. The overall running time of refine is O(nm) plus 0(1) per nonsaturating push plus the time needed for active vertex selection. In the maximum flow case, one of the vertex selection methods we considered was the largestlabel method. In the minimum-cost circulation case, a good method is first-active [46], which is a generalization of the largest-label method. The idea of the method is to process vertices in topological order with respect to the admissible graph. The first-active method maintains a list L of all the vertices of G, in topological order with respect to the current admissible graph GA, i.e., if (v,w) is an arc of GA, v appears before u' in L. Initially L contains the vertices of G in any order. The method consists of repeating the following step until there are no active vertices: Find the first active vertex on L, say v, apply a dischar~q operation to v, and move v to the front of L if the discharge operation has relabeled v. In order to implement this method, we maintain a current vertex v of L. which is the next candidate for discharging. Initially ., is the first vertex on L. The implementation, described in Figure 3.4, repeats the following step until there are no active vertices: If v is active, apply a discharge operation to it, and if this operation relabels v, move v to the front of L: otherwise (i.e., v is inactive), replace v by the vertex currently after it on L. Because the reordering of L maintains a topological ordering with respect to GA, no active vertex precedes v on L. This implies that the implementation is correct. Define a phase as a period of time that begins with v equal to the first vertex on L and ends when the next relabeling is performed (or when the algorithm terminates). Lemma 3.5.1 The first-active method terminates after 0(n 2 ) phases. Proof: Each phase except the last ends with a relabeling operation. I 3.6. REFINEMENT USING BLOCKING FLOWS 41 Theorem 3.5.2 [47] The first-active implementation of refine runs in 0(n O(n 3 3 ) time, giving an log(nC)) bound for finding a minimum-cost circulation. Proof: Lemma 3.5.1 implies that there are O(n 3 ) nonsaturating pushes (one per vert,-x per phase) during an execution of refine. The time spent manipulating L is 0(n) per phase, for a total of 0(n 3 ). M1 other operations require a total of O(nm) time. I A closely related strategy for selecting the next vertex to process is the wave method [40. 46, 47], which gives the same 0(n 3 ) running time bound for refine. (A similar pseudopolynomial algorithm, without the use of scaling and missing some of the implementation details, was developed independently in [11].) The only difference between the first-active method and the wave method is that the latter, after moving a vertex v to the front of L, replaces v by the vertex after the old position of v; if v is the last vertex on L, v is replaced by the first vertex on L. As in the maximum flow case, the dynamic tree data structure can be used to obtain faster implementations of refine. A dynamic tree implementation of the generic version of refine analogous to the maximum flow algorithm discussed in Section 2.5 runs in 0(nm log n) time [47]. A dynamic tree implementation of either the first-active method or the wave method runs in O(nm log(n 2 /m)) time [46]. In the latter implementation, a second data structure is needed to maintain the list L. The details are somewhat involved. 3.6 Refinement Using Blocking Flows An alternative way to implement the refine subroutine is to generalize Dinic's approach to the maximum flow problem. Goldberg and Tarjan [46, 47] showed that refinement can be carried out by solving a sequence of 0(n) blocking flow problems on acyclic networks (i.c., on networks for which the residual graph of the zero flow is acyclic); this extends Dinic's result, which reduces a maximum flow problem to n - 1 blocking flow problems on layered networks. In this section we describe the Goldberg-Tarjan algorithm. At the end of this section, we make a few comrnaents about blocking flow algorithms. To describe the blocking flow version of refine we need some standard definitions. Consider a flow network (G,u,s,t). A flow f is blocking if any path from s to t in the residual graph of zero flow contains a saturated arc, i.e., an arc (v,w) such that u 1 (v,w) = 0. A maximum flow is blocking, but not conversely. A directed graph is layered if its vertices can be assigned integer layers in such a way that layer(v) = layer~w) + 1 for every arc (v,w). A layered graph is acyclic but not conversely. An observation that is crucial to this section is as follows. Suppose we have a pseudoflow f and a price function p such that the vertices can be partitioned into two sets, S and S, such that no admissible arc leads from a vertex in S to a vertex in S; in other words, for every residual arc 42 CHAPTER 3. THE MINIMUM-COST CIRCULATION PROBLEM: COST-SCALING procedure refine(c, f,p); [initialization] c -/2; V(v, w) E E do if cp(v, w) < 0 then f(v, w) ,- u(v, w); [loop] while f is not a circulation do begin S .- {v E VI3u E V such that ef(u) > 0 and v is reachable from u in GA}: Vv E S, p(v) .-- p(v) - c; let N be the network formed from GA by adding a source s, a sink t, an arc (s,v) of capacity e1 (v) for each v E V with ef(v) > 0, and an arc (v,t) of capacity -ef(v) for each v E V with e1 (v) < 0; find a blocking flow b on N; V(v, w) E EA, f(v, w) -- f(v, w) + b(v, w); end; return((, f, p); end. Figure 3.5: The blocking refine subroutine. (v,w) E Ef such that v E S and w E 3, we have cp(v,w) > 0. Define p' to be equal top on S and to p - c on S. It is easy to see that replacing p by p' does not create any new residual arc with reduced cost less than -e. The blocking flow method augments by a blocking flow to create a partitioning of vertices as described above, and modifies the price function by replacing p by p'. Figure 3.5 describes an implementation of refine that reduces e by a factor of two by computing 0(n) blocking flows. This implementation reduces E by a factor of two, saturates all admissible arcs, and then modifies the resulting pseudoflow (while maintaining e-optimality with respect to the current price function) until it is a circulation. To modify the pseudoflow, the method first partitions the vertices of G into two sets S and S, such that S contains all vertices reachable in the current admissible graph GA from vertices of positive excess. Vertices in S have their prices decreased by E. Next, an auxiliary network N is constructed by adding to GA a source s, a sink t, an arc (s,v) of capacity es(v) for each vertex v with ej(v) > 0, and an arc (v,t) of capacity -cS(r) for each vertex with e1 (v) < 0. An arc (v,w) E EA has capacity u 1 (v,w) in N. A blocking flow b on N is found. Finally, the pseudoflow f is replaced by the pseudoflow f'(v, w) = f(v, w) + b(v, w) for (v, w) E E. The correctness of the blocking flow method follows from the next lemma, which can be proved by induction on the number of iterations of the method. Lemma 3.6.1 The set S computed in the inner loop contains only vertices v with e1 (v) > 0. At the beginning of an iteration of the loop, f is an E-optimal pseudoflow with respect to the price function p. Decreasing the prices of vertices in S by E preserves the c-optimality of f. The admissible graph remains acyclic throughout the algorithm. 3.6. REFINEMENT USING BLOCKING FLOWS 43 The bound on the number of iterations of the method follows from Lemma 3.4.5 and the fact that the prices of the vertices with deficit remain unchanged, while the prices of the vertices with excess decrease by c during every iteration. Lemma 3.6.2 The number of iterations of the inner loop in the blocking flow implementation of refine is at most 3n. Since the running time of an iteration of the blocking flow method is dominated by the time needed for a blocking flow computation, we have the following theorem. Theorem 3.6.3 [47] The blocking flow implementation of refine runs in O(nB(n,m)) time, giving an O(nB(n, m)log(nC)) bound for finding a minimum-cost circulation, where B(n, m) is the time needed to find a blocking flow in an acyclic network with n vertices and m arcs. The fastest known sequential algorithm for finding a blocking flow on an acychc network is due to Goldberg and Tarjan [46] and runs in O(mlog(n2 /m)) time. Thus, by Theorem 3.6.3, we obtain an O(nmlog(n2 /m)log(nC)) time bound for the minimum-cost circulation problem. This is the same as the fastest known implementation of the generic refinement method. There is a crucial difference between Dinic's maximum flow algorithm and the blocking flow version of refine. Whereas the former finds blocking flows in layered networks, the latter must find blocking flows in acyclic networks, an apparently harder task. Although for sequential computation the acyclic case seems t) be no harder than the layered case (the best known time bound is O(mlog(n2 /m) for both), this is not true for sequential computation. The Shiloach-Vishkin PRAM blocking flow algorithm [92] for layered networks runs in O(n 2 log n) time using 0(n2 ) memory and n processors. The fastest known PRAM algorithm for acyclic networks, due to Goldberg and Tarjan [49], runs in the same 0(n 2 log n) time bound but uses 0(nm) memory and m processors. 44 CHAPTER 3. THE MINIMUM-COST CIRCULATION PROBLEM: COST-SCALING Chapter 4 Strongly Polynomial Algorithms Based on Cost-Scaling 4.1 Introduction The question of whether the minimum-cost circulation problem has a strongly polynomial algorithm was posed in 1972 by Edmonds and Karp [21] and resolved only in 1984 by Tardos [95]. Her result led to the discovery of a number of strongly polynomial algorithms for the problem [30, 37, 82]. In this chapter we discuss several strongly polynomial algorithms based on cost scaling, in the next, we explore capacity-scaling algorithms, including strongly polynomial ones. Of the known strongly polynomial algorithms, the asymptotically fastest is the capacity-scaling algorithm of Orlin [82]. We begin by describing a modification of the generalized cost-scaling method that makes it strongly polynomial [46]. Then we describe the minimum-mean cycle-canceling algorithm of Goldberg and Tarjan [45]. This simple algorithm is a specialization of Klein's cycle-canceling algorithm [70]; it does not use scaling, but its analysis relies on ideas related to cost scaling. 4.2 Fitting Price Functions and Tight Error Parameters In order to obtain strongly polynomial bounds on the generalized cost-scaling method, we need to take a closer look at the notion of E-optimality defined in Section 3.2. The definition of -optimality motivates the following two problems: 1. Given a pseudoflow f and a constant E > 0, find a price function p such that f is E-optimal with respect to p, or show that there is no such price function (i.e., that f is not C-optimal). 2. Given a pseudoflow f, find the the smallest c > 0 such that f is c-optimal. For this c. we say that f is c-tight. 45 46 CHAPTER 4. STRONGLY POLYNOMIAL ALGORITHMS BASED ON COST-SCALING The problem of finding an optimal price function given an optimal circulation is the special case of Problem 1 with c = 0. We shall see that the first problem can be reduced to a shortest path problem, and that the second problem requires the computation of a cycle of minimum average arc cost. To %.ddressthese p'robleivr, we need some resuilts abut shortest paths and shortest pat' .rc. (see e.g. [96]). Let G be a directed graph with a distinguished source vertex s from which every vertex is reachable and a cost c(v,w) on every arc (v,w). For a spanning tree T rooted at s, the tree cost function d : V --+ R is defined recursively as follows: d(s) = 0, d(v) = d(parent(v)) + c(parent(v), v) for v E V - {s}, where parent(v) is the parent of v in T. A spanning tree T rooted at s is a shortest path tree if and only if, for every vertex v, the path from s to v in T is a minimum-cost path from s to v in G, i.e., d(v) is the cost of a minimum-cost path from s to v. Lemma 4.2.1 (see e.g. [96]) Graph G contains a shortest path tree if and only if G does not contain a negative-cost cycle. A spanning tree T rooted at s is a shortest path tree if and only if c(v, w) + d(v) _ d(w) for every arc (v, w) in G. Consider Problem 1: given a pseudoflow f and a nonnegative c, find a price function p with respect to which f is c-optimal, or show that f is not c-optimal. Define a new cost function c(' ) : E --, R by c(')(v, w) = c(v, w) + c. Extend the residual graph Gf by adding a single vertex s and arcs from it to all other vertices to form an auxiliary graph Gaux = (,aur, Eau) = (0 U {S), Ef U ({s) x V)). Extend c(') to Gaux by defining c( )(s, v) = 0 for every arc (s, v), where v E V. Note that every vertex is reachable from s in Gaux. Theorem 4.2.2 Pseudoflow f is c-optimal if and only if Gaux (or equivalently Gf) contains no cycle of negative c(O)-cost. If T is any shortest path tree of Gaux (rooted at s) with respect to the arc cost function c(' ) , and d is the associated tree cost function, then f is c-optimal with respect to the price function p defined by p(v) = d(v) for all v E V. Proof: Suppose f is c-optimal. Any cycle in Gaux is a cycle in Gf, since vertex s has no incoming arcs. Let F be a cycle of length I in Gar. Then c(F) > -1c, Therefore G., which implies c(,)(F) = c(r) + hc > 0. contains no cycle of negative c()-cost. Suppose Gau,. contains no cycle of negative c(')-cost. Then, by Lemma 4.2.1, Gauj has some shortest path tree rooted at s. Let T be any such tree and let d be the tree cost function. By Lemma 4.2.1, c(')(v,w)+d(v) >_d(w) for all (v,w) E Ef, which is equivalent to c(v,u,)+d(z')-d(w) > -c for all (v,w) E Ef . But these are the c-optimality constraints for the price function p = d. Thus f is c-optimal with respect to p. I 4.3. FIXED ARCS 47 Using Theorem 4.2.2, we can solve Problem 1 by constructing the auxiliary graph G,,,, and finding either a shortest path tree or a negative-cost cycle. Constructing G,,, takes 0(m) time. Finding a shortest path tree or a negative-cost cycle takes O(nm) time using the Bellman-Ford shortest path algorithm (see e.g. [96]). Let us turn to Problem 2: given a pseudoflow f, find the c s-ch that f is c-tight. We need a definition. For a directed graph G with arc cost function c, the minimum cycle mean of G, denoted by p(G,c), is the minimum, over all cycles F in G, of the mean cost of F, defined to be the total arc cost c(F) of r divided by the number of arcs it contains. The connection between minimum cycle means and tight error parameters is given by the following theorem, which was discovered by Engel and Schneider [22] and later by Goldberg and Tarjan [47]: Theorem 4.2.3 [22) Suppose a pseudoflow f is not optimal. Then f is c-tight for c = -u(Gf, c). Proof: Assume f is not optimal. Consider any cycle F in Gf. Let the length of F be 1. For any C, let c(' ) be the cost function defined above: c(')(v, w) = c(v, w) + C for (v, w) E Ef. Let " be such that f is t-tight, and let I = li(Gf,c). By Theorem 4.2.2, 0 < c(W)(F) = c(F) +le, i.e., c(F)/l > -c. Since this is true for any cycle F, p _>-c, i.e., e > -p. Conversely, for any cycle F, c(F)/ which implies c(-,)(r) > 0. By Theorem 4.2.2, this implies -A _>c. _ p, I Karp [66] observed that the minimum mean cycle can be computed in 0(nm) time by extracting information from a single run of the Bellman-Ford shortest path algorithm. This gives the fastest known strongly polynomial algorithm for computing the minimum cycle mean. The fastest scaling algorithm is the 0(yFnmlog(nC))-time algorithm of Orlin and Ahuja [84]. Since we are interested here in strongly polynomial algorithms, we shall use Karp's bound of 0(nm) as an estimate of the time to compute a minimum cycle mean. The following observation is helpful in the analysis to follow. Suppose f is an c-tight pseudoflow and , > 0. Let p be a price function such that f is c-optimal with respect to p, and let F be a cycle in G f with mean cost -c. Since -e is a lower bound on the reduced cost of an arc in G, every arc of F must have reduced cost exactly -e. 4.3 Fixed Arcs The main observation that leads to strongly polynomial cost-scaling algorithms for the minimumcost circulation problem is the following result of Tardos [95]: if the absolute value of the reduced cost of an arc is significantly greater then the current error parameter c, then the value of any optimal circulation on this arc is the same as the value of the current circulation. The following theorem is a slight generalization of this result (to get the original result, take c' = 0). This theorem can be used in two ways. The first is to drop the capacity constraint for an arc of large reduced 48 CHAPTER 4. STRONGLY POLYNOMIAL ALGORITHMS BASED ON COST-SCALING cost. This approach is used in [95]. The second, discussed below, is tc consider the arcs that have the same flow value in every c-optimal circulation for the current value of the ell-or parameter C and to notice that the flow through these arcs will not change. This approach is used in [46, 45]. Theorem 4.3.1 [95] Let c > 0, c' > 0 be constants. Suppose that a circulation f is c-optimal with respect to a price function p, and that there is an arc (v,w) E E such that jcp(v,w) >_n(c +e). Then, for any c-optimal circulation f', we have f(v,w) = f'(v,w). Proof: By antisymmetry, it is enough to prove the theorem for the case cp(v,w) > n(c + '). Let f' be a circulation such that f'(v,w) 0 f(v,w). Since cp(v,w) > E, the flow f through the arc (v,w) must be as small as the capacity constraints allow, namely -u(w,v), and threfore f'(v,w) 5 f(v,w) implies f'(v,w) > f(v,w). We show that f' is not C'-optimal, and the theorem follows. Consider G> = {(x, y) E Elf'(x, y) > f(x, y)}. Note that G> is a subgraph of Gf, and (v, u,) is al. arc of G>. Since f and f' are circulations, G> must contain a simple cycle F that passes through (v, w). Let I be the length of F. Since all arcs of F are residual arcs, the cost of F is at least cp(v, w) - (l 1)c > n(c + c') - (n - 1)c > nc'. Now consider the cycle F obtained by reversing the arcs on F. Note that F is a cycle in G< = {(x, y) E Elf'(x, y) < f(x, y)} and is therefore a cycle in Gf,. By antisymmetry, the cost of I is less than -nt' and thus the mean cost of T is less than -c'. But Theorem 4.2.3 implies that f' is not c-optimal. I To state an important corollary of Theorem 4.3.1, we need the following definition. We say that an arc (v, w) E E is c-fixed if the flow through this arc is the same for all c-optimal circulations. Corollary 4.3.2 [46] Let c > 0, suppose f is c-optimal with respect to a price function p, and suppose that (v,w) is an arc such that Icp(v,w)l > 2nc. Then (v,w) is c-fixed. Define F, to be the set of #-fixed arcs. Since the generalized cost-scaling method decreases C, an arc that becomes c-fixed stays c-fixed. We show that when E decreases by a factor of 2n, a new arc becomes c-fixed. Lemma 4.3.3 Assume c' < -L. Suppose that there exists an c-tight circulation f. Then F,, properly contains F,. Proof: Since every c'-optimal circulation is c-optimal, we have F, C F,. To show that the containment is proper, we have to show that there is an c'-fixed arc that is not (-fixed. 4.4. THE STRONGLY POLYNOMIAL FRAMEWORK 49 Since the circulation f is (-tight, there exists a price function p such that f is c-optimal with respect to p, and there exists a simple cycle r in Gf every arc of which has reduced cost -c. (See Section 4.2.) Since increasing f along r preserves c-optimality, the arcs of F are not c-fixed. We show that at least one arc of F is e-fixed. Let f' be a circulation that is '-optimal with respect to some price function p'. Since the mean cost of F is -c, there is an arc (v, w) of F with cp,(v, w) < -c < -2n'. By Corollary 4.3.2, the arc (v, w) is e-fixed. I In the nexi section we show how to use this lemma to get a strongly polynomial bound for a variation of the generalized cost-scaling method. 4.4 The Strongly Polynomial Frariework The minimum-cost circulation framework of Section 3.3 has the disadvantage thM ' the number of iterations of refine depends on the magnitudes of the costs. If the costs are huge intege.:, the method need not run in time polynomial in n and m; if the costs are irrational, the method need not even tcrminat.. In this section we show that a natural modification of the generalized costscaling approach produces strongly polyn&...A,- algorithms. The running time bounds we derived for algori..hms based on the approach of Section 3.3 remain valid for the niodified approach presented 'a 1ii section. The main idea of this modification is to improve e periodically by finding a price function that fits the current circulation better than the current price function. This idea can also be used to improve the practical performance of the method. The changes needed to make the generalized cost-scaling approach strongly polynomial, suggested by Lemma 4.3.3, are to add an extra computation to the main loop of the algorithm and to change the termination condition. Before calling refine to reduce the error parameter E, the new method computes the value A and a price function pA such that the current circulation f is A-tight with respect to p,\. The strongly polynomial method is described on Figure 4.1. The value of A and the price function pA in line (*) are computed as described in Section 4.2. The algorithn terminates when the circulation f is optimal, i.e., \ = 0. The time to perform line (*) is O(nm). (See Section 4.2.) Since all the implementations of refine that we have considered have a time bound greater than O(nm), the time per iteration in the new version of the algorithm exceeds the time per iteration in the original version by less than a constant factor. Since each iteration at least halves c, the bound of O(log(nC)) on the number of iterations derived in Chapter 3 remains valid, assuming that the costs are integral. For arbitrary real-valued costs, we shall derive a bound of O(mlog n) on the number of iterations. Theorem 4.4.1 O(m log n). 145] The total number of iterations of the while loop in procedure "Zin-cost is Proof: Consider a time during the execution of the algorithm. During the next O(log n) iterations, 50 CHAPTER 4. STRONGLY POLYNOMIAL ALGORITHMS BASED ON COST-SCALING procedure min-cost(V, E, u, c); [initialization] c - C; Vv, p(v) - 0; V(v,w) E E, f(v,w) .- 0; [loop] whilf- c > 0 do begin (*) find A and px such that f is A-tight with respect to pX; if A > 0 then (,f, p) - refine(A,f, px) else return(f); end. Figure 4.1: The strongly polynomial algorithm. either the algorithm terminates, or the error parameter is reduced by a factor of 2n. In the latter case, Lemma 4.3.3 implies that an arc becomes fixed. If all arcs become fixed, the algorithm terminates in one iteration of the loop. Therefore the total number of iterations is O(m log n). I The best strongly polynomial implementation of the generalized cost-scaling method [46], based on the dynamic tree implementation of refine, runs in O(nm 2 log(n 2 /m)log n) time. 4.5 Cycle-Canceling Algorithms The ideas discussed in Sections 4.2 - 4.4 are quite powerful. In this section we use ti.ese ideas to show that a simple cycle-canceling algorithm of Klein [70] becomes strongly polynomia ' a carpful choice is made among possible cycles to cancel. Klein's algorithm consists of repeatedy finding a residual cycle of negative cost and sending as much flow as possible around the cycle. This algorithm can run for an exponential number of iterations if the capacities and costs are integers, and it need not terminate if the capacities are irrational [26]. Goldberg and Tarjan [45] showed that if a cycle with the minimum mean cost is canceled at each iteration, the algorithm becomes strongly polynomial. We call the resulting algorithm the minimum-mean cycle-canceling algorithm. The minimum-mean cycle-canceling algorithm is closely related to the shortest augmenting path maximum flow algorithm of Edmonds and Karp [21]. The relationship is as follows. If a maximum flow problem is formulated as a minimum-cost circulation problem in a standard way, then Klein's cycle-canceling algorithm corresponds exactly to the Ford-Fulherson maximum flow algorithm, and the minimum-mean cycle-canceling algorithm corresponds exactly to the Edmonds-Karp algorithm. The minimum-mean cycle-canceling algorithm can also be interpreted as a steepest descent method using the L 1 metric. For a circulation f we define ((f) to be zero if f is optimal and to be the unique number (' > 0 such that f is c-tight otherwise. We use c(f) as a measure of the quality of f. Let f be an arbitrary 4.5. CYCLE-CANCELING ALGORITHMS 51 circulation, let c = c(f), and let p be a price function with respect to which f is c-optimal. Holding c and p fixed, we study the effect on E(f) of a minimum-mean cycle cancellation that modifies f. Since all arcs on a minimum-mean cycle have negative reduced cost with respect to p. cancellation of such a cycle does not introduce a new residual arc with negative reduced cost, and hence ((.f) does not increase. Lemma 4.5.1 A sequence of m minimum-mean cycle cancellations reduces (f) to at most (1 - I/n)c, i.e., to at most 1 - 1/n times its original value. Proof: Let p a price function such that f is c-tight with respect to p. Holding c and p fixed, we study the effect on the admissible graph GA (with respect to the circulation f and price function p) of a sequence of rn minimum-mean cycle cancellations that modify f. Initially every arc (v, u) e EA satisfies cp(v, w) > -. Canceling a cycle all of whose arcs are in EA adds only arcs of positive reduced cost to EJ and deletes at least one arc from EA. We consider two cases. Case 1: None of the cycles canceled contains an arc of nonnegative reduced cost. Then each cancellation reduces the size of EA, and after m cancellations EA is empty, which implies that f is optimal, i.e., c(f) = 0. Thus the lemma is true in this case. Case 2: Some cycle canceled contains an arc of nonnegative reduced cost. Let F be the first such cycle canceled. Every arc of F has a reduced cost of at least -C, one arc of F has a nonnegative reduced cost, and the number of arcs in F is at most n. Therefore the mean cost of F is at least -(1 - 1/n)c. Thus, just before the cancellation of F, E(f) _<(1 - 1/n)c by Theorem 4.2.3. Since c(f) never increases, the lemma is true in this case as well. I Lemma 4.5.1 is enough to derive a polynomial bound on the number of iterations, assuming that all arc costs are integers. Theorem 4.5.2 [45] If all arc costs are integers, then the minimum-mean cycle-canceling algorithm terminates after O(nrn log(nC)) iterations. Proof: The lemma follows frm Lemmas 3.2.1 and 4.5.1 and the observation that the initial circulation is C-optimal. I To obtain a strongly polynomial bound, we use the ideas of Section 4.3. The proof of the next theorem uses the following inequality: (1 - I),,(Inn+l) < - for n > 2. Theorem 4.5.3 [45] For arbitrary real-valued arc costs, the minimum-mean cycle-canceling algorithm terminates after O(nm2 log n) cycle cancellations. 52 CHAPTER 4. STRONGLY POLYNOMIAL ALGORITHMS BASED ON COST-SCALfING Proof: Let k = m(n Pn n + 11). Divide the iterations into groups of k consecutive iterations. We claim that each group of iterations fixes the flow on a distinct arc (v, w), i.e., iterations after those in the group do not change f(v, w). The theorem is immediate from the claim. To prove the claim, consider any group of iterations. Let f be the flow before the first iteration of tba group, f' the flow after the last iteration of the group, c = c(f), c = c(f'), and let p' be a price function for which f' satisfies the c'-optimality constraints. Let F be the cycle canceled in th,, first iteration of the group. By Lemma 4.5.1, the choice of k iniplier that c _ (1 Since the mean cost of F is -c, some arc on F, say (v,w), must have cp,(v,w) - n),l1,+11 K -c < -2nc. By Corollary 4.3.2, the flow on (v,w) will not be changed by iterations after those in the group. But f(v, w) is changed by the first iteration in the group, which cancels F. Thus each group fixes the flow on a distinct arc. I Theorem 4.5.4 [45] The minimum-mean cycle-canceling algorithm runs in O(n 2 m 31og,) time on networks with arbitrary real-valued arc costs, and in O(n 2 M 2 min{log(nC),m log n}) time on networks with integer arc costs. Proof: Immediate from Theorems 4.5.2 and 4.5.3. I Alth ugh the minimum-mean cycle-canceling algorithm seems to be of mostly theoretical inLerest, it has a variant that is quite efficient. This variant maintains a price function and. instead of canceling the minimum-mean-cost residual cycle, cancels residual cycles composed entirely of negative reduced cost arcs; if no such cycle exists, the algorithm updates the price function to improve the error parameter c. An implementation of the algorithm using the dynamic tree data structure runs in O(nm log nlog(nC)) time. See [45] for details. The techniques used to analyze the minimum-mean cycle-canceling algorithm also provide an approach to making the primal network simplex algorithm for the minimum-cost circulation problem more efficient. Tarjan [99] has discovered a pivot rule that gives a bound of O(n'Iog1+(Z )) on the number of pivots and the total running time. This is the first known subexponential time bound. For the extended primal network simplex method in which cost-increasing as well as costdecreasing pivots are allowed, he has obtained a strongly polynomial time bound. In contrast, the dual network simplex method is known to have strongly polynomial variants, as mentioned in Chapter 5. Barahona and Tardos [7] have exhibited another cycle-canceling algorithm that runs in polynomial time. Their algorithm is based on an algorithm of Weintraub [103], which works as follows. Consider the improvement in the cost of a circulation obtained by canceling a negative cycle. Since a symmetric difference of two circulations can be decomposed into at most m cycles, canceling the cycle that gives the best improvement reduces the difference between the current and the optimal values of the cost by a factor of (1 - I/m). If the input data is integral, only a polynomial number 4.5. CYCLE-CANCELING ALGORITHMS 53 of such improvements can be made until an optimal solution is obtained. Finding the cycle that gives the best improvement is NP-hard, however. Weintraub shows how to find a collection of cycles whose cancellaion reduces the cost by at least as much as the best improvement achievable by canceling a single cycle. His method requires a superpolynomial number of applications of an algorithm for the assignment problem. Each such application yields a minimum-cost collection of vertex-disjoint cycles. Barahona and Tardos [71 have shown that the algorithm can be modified so that the required collection of cycles is found in at most m assignmen, computations. The resulting minimum-rost circulation algorithm runs in polynomial time. 54 CHAPTER 4. STRONGLY POLYNOMIAL ALGORITHMS BASED ON COST-SCALING Chapter 5 Capacity-Scaling Algorithms 5.1 Introduction In this chapter, we survey minimum-cost circulation algorithms based on capacity scaling. We concentrate on two algorithms: the first polynomial-time algorithm, that of Edmonds and Karp [21]. who introduced the idea of scaling, and the algorithm of Orlin [82]. The latter is an extension of the Edmonds-Karp algorithm and is the fastest known strongly polynomial algorithm. in this chapter we shall consider the (uncapacitatedi transshipment problem. which is equivalent to the minimum-cost circulation problem. This simplifies the presentation considerably. We begin by describing a generic augmenting path algorithm, which we call the minimum-cost augmentation algorithm, that is the basis of most capacity-scaling algorithms. The algorithm is due to Jewel [61], Busacker and Gowen [14], and Iri [58]. The use of dual variables as described below was proposed independently by Edmonds and Karp [211 and Tomizawa [100]. The idea of the algorithm is to maintain a pseudoflow f that satisfies the complementary slackness constraints while repeatedly augmenting flow to gradually get rid of all excesses. Two observations justify this method: (1) augmenting flow along a minimum-cost path preserves the invariant that the current pseudoflow has minimum cost among those pseudoflows with the same excess function; (2) a shortest path computation suffices both to find a path along which to move flow and to find price changes to preserve complementary slackness. The algorithm maintains a pseudoflow f and a price function p such that cp(v,w) >_ 0 for every residual arc (v, w). The algorithm consists of repeating the auqmentation step, described in Figure 5.1, until every excess is zero. Then the current pseudoflow is optimal. The correctness of this algorithm follows from the observation that the price transformation in the augment step preserves complementary slackness and makes the new reduced cost of any minimum-cost path from s to t equal to zero. One iteration of the augment step takes time proportional to that required by a single-source shortest path computation on a graph with arcs of non-negative cost. The fastest known strongly polynomial algorithm for this computation is 55 CHAPTER 5. CAPACITY-SCALING ALGORITHMS 56 augment(s, t). Applicability: e(s) > 0 and e(t) < 0. Action: For every vertex v, compute 7r(v), the minimum reduced cost of a residual path from s to v. For every vertex v, replace p(v) by p(v) + 7r(v). Move a positive amount of flow from s to t along a path of minimum reduced cost. Figure 5.1: The augmentation step. Fredman and Tarjan's implementation [29] of Dijkstra's algorithm, which runs in O(nlogn + m) time. The problem can also be solved in O(m log log C) time [85] or O(nV'/Wg_ + m) time [2]. where C is the maximum arc cost. Since we are mainly interested here in strongly polynomial algorithms, we shall use O(n log n + m) as our estimate of the time for each augment step. Remark: The only reason to maintain prices in this algorithm is to simplify the minimum-cost path computations by guaranteeing that each such computation is done on a graph all of whose arc costs are non-negative. If prices are not maintained in this way, each shortest path computation takes O(nm) time using the Bellman-Ford algorithm (see e.g. [96]). 5.2 The Edmonds-Karp Algorithm The overall running time of the augmentation algorithm depends on the number of augmentation steps, which cannot be polynonially bounded without specifying their order more precisely. To impose an efficient order, Edmonds and Karp introduced the idea of capacity scaling, which Orlin [80] in his description of Edmonds-Karp algorithm reinterpreted as excess scaling. We shall present a modification of Orlin's version of the Edmonds-Karp algorithm. The algorithm maintains a scaling parameter A such that the flow on every arc is an integer multiple of A. For a given pseudoflow f and value of the scaling parameter A, we denote the sets of vertices with large excesses and large deficits as follows: S 1 (A) = {v E V: e 1 (v) > A}; Tf(A) = {v E V ej(v) < -A}. (5.1) The algorithm consists of a number of scaling phaseF. During a A-scaling phase the residual capacity of every arc is an integer multiple of A. The algorithm chooses vertices s E Sf(A/2), t E Tj(A/2), performs augnient(s, t), which sends A units of flow along a minimum-cost path from s to t,and repeats such augmentations until either Sf(A/2) or T(A/2) is empty. Note that this can create a deficit in place of an excess, and vice versa. Vertices with the new excesses and deficits, however, are not in Sj(A/2) UTs(A/2). Then the algorithm halves A and begins a new scaling phase. The initial value of A is 2 [IogD ] . The algorithm terminates after the 1-scaling phase. assuming all supplies are integral. Figure 5.2 provides a detailed description of a phase of the 5.2. THE EDMONDS-KARP ALGORITHM while S1 (A/2) 96 0 and Tf (A/2) :0 0 do begin 57 Choose s E. S1(A/2 ) and f E Tf (A/2). Perform augment (s,tf), sending A units of flow along the augmenting path selected. end; Figure 5.2: A phase of the Edmonds-Karp Algorithm. algorithm. The following lemma follows from the description of the algorithm. Lemma 5.2.1 During a A-scaling phase, the residual capacity of every arc is an integer multiple of A. Lemma 5.2.2 Let f' be the pseudoflow at the end of a 2A-scaling phase. If Sf4(A) = 0, then there are at most IS,(A/2)I flow augmentations during the A-scaling phase. An analogous statement is true if Tf,(A) = 0 at the end of a 2A-scaling phase. Proof: By assumption, Sf,(A) = 0. Therefore, for every vertex v E V, ej(v) < A during the A-scaling phase. This implies that pushing A units of flow along a path from s e S1 (A/2) to t E T(A/2) removes s from S (A/2). Note that since ej(v) < -A/2 for v E T(A/2), vertices with 1 new excesses are not in S 1 (A/2). There are [log D1 phases, each of which consists of up to n single-source shortest path computations on networks with non-negative arc costs. Thus we have the following results for networks with integral supplies: Theorem 5.2.3 [21] The Edmonds-Karp algorithm solves the transshipment problem in O((n log D) (n log n + m)) time, and the minimum-cost circulation problem in O((m log U)(n log n + n)) time. The excess-scaling idea can be combined with the idea of augmenting along approximately minimum-cost paths rather than exactly minimum-cost paths. Approximately minimum-cost paths can be defined using the c-optimality notion discussed in Chapter 3. An appropriate combination of these ideas yields the double scaling algorithm of Ahuja, Goldberg, Orlin, and Tarjan [1]. An implementation of this algorithm based on dynamic trees has a time bound of O(nm(loglog U) log(nC)) for the minimum-cost circulation problem. 58 CHAPTER 5. CAPACITY-SCALING ALGORITHMS 5.3 Strongly Polynomial Algorithms After Tardos [95] discovered the first strongly polynomial algorithm, Fujishige [30] and independently Orlin [80] developed more efficient strongly polynomial algorithms. Orlin's algorithm [821, that is the main focus of this section, is an extension of these algorithms. These algorithms are based on the method of Section 5.2 and the following dual version of Theorem 4.3.1. Theorem 5.3.1 [30, 80] Let f be a pseudoflow, and let p be a price function such that the complementary slackness conditions are satisfied and f(v, w) > ev:es(v)> Then every optimal price 0 ef(v)j. function p" satisfies cp.(v, w) = 0. Proof: Let p" be an optimal price function. and let f" bc a corresponding optimal pseudoflow. Assume by way of contradiction that cp.(w, v) 00. Define the pseudoflow f' by f'(x, y) = f(x, y)f*(x,y) (i.e., we can obtain an optimal pseudoflow by augmenting f by f'). The Decomposition Theorem (Theorem 1.7.2) imlies that f' r:n h,- drompospd into flows along a collection of cycl-; and a collection of paths from vertices with excess to vertices with deficits (both excesses and deficits are with respect to f). Furthermore, the Decomposition Theorem also implies that for any arc (x, y) that is not on one of the cycles, If'(x, y) I : v:e(v)>O ef(v). The arcs on the cycles of the decomposition are in Ef, and the arcs opposite to the ones on the cycles are in Ef.. This implies that the cycles have zero cost, and that every arc on the cycles must have zero reduced cost with respect to both p and p*. We have assumed that cp.(V, w) $ 0; therefore (w,v) cannot be on such a cycle. Thus f'(w,v) < Zv:ef(v)>oef(v) and thus f'(v,w) > 0. But f (v, w) > 0 and cp.(v, w) $ 0 contradict the complementary slackness constraint. I The idea behind the strongly polynomial algorithms based on capacity-scaling is to contract the arcs satisfying the condition c,. Theorem 5.3.1. Let f and p be a pseudoflow and a price function, respectively, and let (v, w) be an arc satisfying the condition of Theorem 5.3.1. The reduced cost function cp defines an equivalent problem with cp(v,w) = 0. By the above theorem the optimal prices of the vertices v and w are the same. Define a transshipment problem on the network formed by contracting the arc (v,w), with cost function cp, capacity function u, and demand function d such that the demand of the new vertex is d(v) + d(w). Theorem 5.3.2 For any optimal price function p* for the new problem, the price function p'(r) = p(r) + p*(r) for r E V (with p'(v) = p*(w) defined to be the price of the new vertex) is optimal for the original problem. Proof: The optimal price function is the optimal solution of a linear program. the linear programming dual of the minimum-cost flow problem. The price function p' = p + p* is the optimal 5.3. STRONGLY POLYNOMIAL ALGORITHMS 59 Step 1 Run the Edmonds-Karp algorithm for the first [2log n] scaling phases with cost function cp, starting with A = max11v Id(v)I. Let f' be the pseudoflow and let p' be the price function found by the algorithm. Step 2 Contract every arc (v,w) with Jf(v,w)] > nA and update the price function by setting p(v) .- p(v) + p'(v) for all v E V (where all vertices v E V contracted into the same vertex v' of the current network have the same price p'(v) = p'(v')). Figure 5.3: The inner loop of the simple strongly polynomial algorithm. solution of same linear program with the additional constraint c(v, w) + p'(v) - p'(w) = 0. By Theorem 5.3.1, the two linear programs have the same set of optimal solutions. I We shall describe a simple strongly polynomial algorithm based on this idea. It is a variation of Fujishige's algorithm [30], though our description is simplified by using ideas from [82]. One iteration of the algorithm consists of running the Edmonds-Karp algorithm for the first 2lognz scaling phases starting with A = maxVEV Id(v)I (we cannot find the power of two used by the Edmonds-Karp algorithm in strongly polynomial time), and then contracting all arcs that satisfy the condition of Theorem 5.3.1. The next iteration considers the contracted network. We shall prove that -ach iteration will contract at least one arc. The algorithm terminates when the current pseudoflow is a feasible solution. Theorem 5.3.2 guarantees that the price function found at this point is optimal. Given an optimal price function, the optimal flow can be found by a single maximum flow computation. (See Figure 5.3 for a more detailed description.) The following theorem implies that at least one arc is contracted at each iteraion of the inner loop. Theorem 5.3.3 [82] Let f and A be the pseudoflow and the scaling parameter at the end of an iteration of the algorithm of Figure 5.3. Then EvEv Ief (v)l < nA. Furthermore, if f is not feasible, then there exists an arc (v, w) with If(v, w)I > nA. Proof: After the A-scaling phase either Sf(A/2) or Tf(A/2) is empty. For a pseudoflow f, the excesses sum to zero. Therefore FvV Ief(v)I < nA. Also, every excess is less than nA. After [2logn] scaling phases, A < (maxWEV Id(w)I)/n 2 . For a vertex v whose demand has the maximum absolute value this implies that IE f(v,w)l > Id(v)l - nA > n(n - 1)A. (5.2) wEV Consequently, at least one arc incident to v carries at least nA flow. 3 One iteration takes O(n(n log n + m)log n) time. Each iteration contracts at least one arc, and therefore there are at most n iterations. 60 CHAPTER 5. CAPACITY-SCALING ALGORITHMS Theorem 5.3.4 The algorithm of Figure 5.3 solves the transshipment problem in O(n 2 log n(n log n + m)) time, and the minimum-cost circulation problem in O(m'logn(nlogn + m)) time. This simple algorithm wastes a lot of time by throwing away the current flow after each contraction. If there are several vertices with demand close to the maximum then it might make sense to run somewhat more than 2 log n scaling phases. The proof of Theorem 5.3.3 would then apply to more than one vertex, and more than one arc could be contracted in an iteration. This idea does not help much if there are not enough vertices with close-to-maximum demand. On the other hand, by Lemma 5.2.2 the number of shortest path computations during a phase can be bounded in terms of the number of vertices with relatively large demand (it is bounded either by IS(A/2)I or by IT(A/2)I). One could try to find a point for stopping the scaling algorithm that balances the number of shortest path computations done and the number of arcs contracted. Gail and Tardos [37] developed an O(n 2 (n log + m)log n)-time minimum-cost circulation algorithm based on this idea. In fact, the balancing technique of Galil and Tardos [37], together with Orlin's [82] proof of Theorem 5.3.3, can be used to give an O(n(n log n + m) log n)-time algorithm for the transshipment problem. Instead of pursuing this approach further, however, we shall present a variant of a much simpler algorithm due to Orlin [82] that has same efficiency and does not rely on delicate balancing. This algorithm contracts arcs during the scaling phases, without starting the scaling from scratch after each contraction; the algorithm finds an optimal flow directly, without using an optimal price function. The algorithm is a modification of the Edmonds-Karp algorithm, but runs in O(n log n) iterations instead of the O(n log D) iterations of the Edmonds-Karp algorithm. Orlin'-, algorithm runs the Edmonds-Karp scaling algorithm starting with A = max,Ev Id(v)l and ccntracts all arcs that carry at least 4nA flow. The scaling algorithm continues as long as the current pseudoflow is non-zero on some uncontracted arc. When the pseudoflow is zero on every arc, the scaling is restarted by setting A to be the maximum absolute value of a demand. Roughly speaking, a vertex v is the start of at most 2 log n shortest path computations during the algorithm. For a vertex v of the original graph, this follows from the fact that v V Sf(A) U Tj(A) unless Id(v)l > A, and therefore an arc incident to v will be contracted at most 2logn scaling phases after v first served as a starting vertex. (See the proof of Theorem 5.3.3.) This may rot be true for the contracted vertices, however. A vertex created by a contraction may have very small demand relative to its current excess. We say that a vertex v is active if v E S1 (A/2) U T1 (A/2). A vertex v is activated by a Ascaling phase if v is not active at the end of the 2A-scaling phase, and it becomes active during the A-scaling phase (i.e., either it becomes active at the beginning of the phase or the vertex is created by contraction during the phase). We can prove the following lemma using a proof similar to that of Lemma 5.2.2. Lemma 5.3.5 The number of shortest path computations during a phase is bounded by the number of vertices activated by the phase. 5.3. STRONGLY POLYNOMIAL ALGORITHMS 61 Id(v)I. If A = 0 then STOP. Step 2. while S1 (A/2) 6 0 and Tj(A/2) 0 0 do begin Step 1. Let A = maxEv Choose 8 E S1(A/2) and t E T 1 (A/2). Perform augment(s,1), sending A units of flow along the augmenting path selected. Contract every arc (v, w) with f(v, w) > 4nA. end. Step 3. If f is zero on all uncontracted arcs go to Step 1, otherwise let A = A/2 and go to Step 2. Figure 5.4: Orlin's Algorithm. Theorem 5.3.6 [82] A vertex can be activated at most r2lognl + 1 times before it is contracted. Proof: A vertex can be activated once due to contraction. When a vertex v is activated for the second time, it must already exist at the end of the previous scaling phase, when it was not active. Let A be the scaling parameter in the phase when v is activated for the second time. At the beginning of this phase, A/2 < Jej(v)l < A. But d(v) - ej(v) is an integer multiple of 2A. This implies that A/2 < lef(v)l < Id(v)I. After O(logn) more scaling phases the scaling parameter 2 will be less than td(v)i/4n . Using ai, argument analogous to the proof of Theorem 5.3.3 we can conclude that some arc incident to v is contracted. I The previous lemma and theorem bound the number of shortest path computations during the algorithm. All other work done is linear per scaling phase. At least one arc is contracted in each group of O(log n) scaling phases. Therefore, there are at most O(n log n) scaling phases. Theoreiii 5.3.7 [82] Orlin's algorithm solves the transshipment problem in O(n logmin{n, D}(n log n + in)) time, and the minimum-cost circulation problem in O(mlog min{n, U}(n log n + m)) time. Remark: In contrast to the simpler strongly polynomial algorithm discussed earlier, Orlin's algorithm constructs an optimal pseudoflow directly, without first constructing an optimal price function. To see this, consider the amount of flow moved during the A-scaling phase for some value A. Lemma 5.3.5 gives a 2nA bound. Suppose an arc is contracted during the A-scaling phase. The overall amount of flow that is moved after this contraction can be bounded by 4nA. Therefore the pseudoflow constructed by the algorithm is feasible in the uncontracted network. Orlin [80] has observed that capacity scaling ideas can be used to guide pivot selection in the dual network simplex algorithm for the transshipment problem. More recently Orlin, Plotkin and Tardos (personal communication, 1989) have obtained an O(nn log min(n, D)) bound on the number of pivots based on such a pivot selection strategy. 62 CHAPTER 5. CAPACITY-SCALING ALGORITHMS Chapter 6 The Generalized Flow Problem 6.1 Introduction The generalized flow problem models the following situation in financial analysis. An investor wants to take advantage of the discrepancies in the prices of securities on different stock exchanges and of currency conversion rates. His objective is to maximize his profit by trading on different exchanges and by converting currencies. The generalized flow problem, considered in this chapter, models the above situation, assuming that a bounded amount of money is available to the investor and that bounded amounts of securities can be traded without affecting the prices. Vertices of the network correspond to different currencies and securities, and arcs correspond to possible transactions. The generalU7ed flow prhbrn was first considered by Jewell [61]. This problem is very similar to the minimum-cost circulation problem, and several of the early minimum-cost circulation algorithms have been adapted to this problem. The first simple combinatorial algorithms were developed by Jewell [61] and Onaga [79]. These algorithms are not even pseudopolynomial, however, and tor real-valued data they need iot terminate. Several variations suggested in the early 70's result in algorithms running in finite (but exponential) time. The paper by Truemper [101] contains an excellent summary of these results. The generalized flow problem is a special case of linear programming, and therefore it can be solved in polynomial time using any general-purpose linear programming algorithm, such as the ellipsoid method [69] and the interior-point algorithms [65, 88]. Kapoor and Vaidya [64] developed techniques to use the structure of the matrices that arise in linear programming formulations of network flow problems to speed up interior-point algorithms. The first two polynomial-time combinatorialalgorithms were developed by Goldberg, Plotkin, and Tardos [42]. We shall review one of their algorithms in detail and sketch the other one. The current fastest algorithm for the generali d flow problem, due to Vaidya, is based on an interior-point method for linear programming and runs in O(n 2 m'1 5 log(nB)) time [102]. This algorithm combines the ideas from the paper of Kapoor and Vaidya [64] with the current fastest 63 64 CHAPTER 6. THE GENERALIZED FLOW PROBLEM linear programming algorithm (which uses fast matrix multiplication). The two combinatorial 2 algorithms due to Goldberg, Plotkin and Tardos [42] run in O(mn (M + n log n) log n log B) and O(m 2 n 2 log n log 2 B) time, respectively. More recently Maggs (personal communication, 1989) improved the latter bound through better balancing to O (n2m2log2B og((n2 logn)/m)+Iog ogB) 6.2 Simple Combinatorial Algorithms Jewell [61] and Onaga [78] suggested solving the generalized flow problem by using adaptations of augmenting path algorithms for the maximum flow and minimum-cost circulation problems. Theorem 1.6.2 is the basis of Onaga's augmenting path algorithm for the restricted problem. Starting from the zero flow, this algorithm iteratively augments the flow along a flow-generating cycle in the residual graph, thereby increasing the excess at the source. The structure of the restricted problem makes it possible to find the most efficient flow-generating cycle, i.e., the residual flow-generating cycle with highest-gain. The highest-gain cycle consists of an arc (r,s) and a highest-gain simple path from s to r for some vertex r. The highest-gain simple path to r is the shortest path, if the length of each arc (v,w) is defined as l(v,w) = -log7(v,w). Such a shortest path exists since deleting the arcs entering the source from the residual graph of the initial zero flow yields a graph with no negative cycles. The main observation of Onaga is that if the augmentation is done using the most-efficient flow-generating cycle, then all flow-generating cycles in the residual graph of the resulting generalized flow pass through the source. Onaga's algorithm iteratively augments the generalized flow along the most efficient flow-generating cycle in the residual graph. This algorithm maintains the invariant that every vertex is reachable from s in the current residual graph., a property that can be verified by induction on the number of augmentations. Consider the special case of a network with two diktingui-hed vertices s,t E V and with all gains other than y(t,s) and y(s, t) equal to one. The generalized flow problem in this network is equivalent to the maximum flow problem, and Onaga's algorithm specializes to the Ford-Fulkerson maximum flow algorithm. Consequently, the algorithm does not run in finite time. The next algorithm we describe uses a maximum flow computation as a subroutine. To describe this algorithm, we need to introduce a relabeling technique of Glover and Klingman [38]. This technique can be motivated as follows. Recall the financial analysis interpretation of the generalized flow problem, in which vertices correspond to different securities or currencies, and arcs correspond to possible transactions. Suppose one country decides to change the unit of currency. (For example. Great Britain could decide to introduce the penny as the basic currency unit, instead of the pound. or Italy could decide to erase a couple of O's at the end of its bills.) This causes an appropriate update of the exchange rates. Some of the capacities change as well (a million £ limit on the £ - DM exchange would now read as a limit of 100 million pennies). It is easy to see that such a relabeling defines an equivalent problem. Such a relabeling can be used, for example, to normalize 6.2. SIMPLE COMBINATORIAL ALGORITHMS 65 local units of currency to current market conditions. To formally define the relabeled problem, let /I(v) E R+ denote the number of old units corresponding to each new unit at vertex v E V. Given a function y, we shall refer to p(v) as the label of V. Definition: For a function : V R+ and a network N = (V,E,y,u), the relabeled network is Nm = (V, E,-y., u,), where the relabcled capacities and the relabeled gains are defined by ' uA(v, w) = U(VW)/A(V) -Y(V,W) = -Y(v'W)P(v)/P(w). For a generalized psc'-doflow g and a labeling i, the generalized pseudoflow corresponding to g in the relabeled network is g,(v,w) = g(v,w)/p(v), the relabeled residual capacity is defined by ug,,(v,w) = (u(v,w) - g(v,w))/p(v), and the relabeled excess by e9,(v) =eg(v)/;4L'). Thc corresponding pseudoflows have the same residual graph. Now we present a canonical relabeling. The residual graph of a generalized pseudoflow g can be canonically relabeled if every vertex v E V is reachable from the source and every flow-generating cycle in the residual graph contains the source. For a vertex v E V, the canonical label P(t,) is defined to be the gain of a highest-gain simple path from s to v in the residual graph. That is. one new unit corresponds to the amount of flow that can reach the vertex v if one old unit of flow is pushed along a most-efficient simple path in the residual graph from s to v, ignoring capacity restrictions along the path. Observe that in a restricted network, the highest-gain paths from s to each other vertex can be found using any single-source shortest path algorithm. Theorem 6.2.1 After a canonical relabeling, the following properties hold: Every arc (V,w) G E, such that w y s has -y,(v,w) _<1; there exists a path from s to every other vertex r in the residual graph with y(v, w) = 1 for all arcs (v, w) on the path; the most efficient flow-generating cycles each consist of an arc (r,s) E Eq and an (s,r)-path for some r E V, such that -Y(t,W) = 1 along the path and -y.(r,s) = max(Q,(V,w): (v,w) E Eg). Next we describe a simple finite algorithm, due to Truemper [101]. The algorithm, described in Figure 6.1, is a refinement of Onaga's algorithm in which augmentation along all of the maximum gain cycles is done simultaneously. The algorithm maintains a generalized pseudoflow g whose residual graph has every vertex reachable from the source and has all flow-generating cycles passing through the source. The algorithm first canonically relabels the residual graph. Now the highestgain residual cycles have maxymum relabeled gain on arcs entering the source and have a relabeled gain of one on all other arcs. Consider the subgraph induced by arcs with relabeled gain of one and arcs of maximum relabeled gain entering the source. A circulation that maximizes the sum of the flow on the latter arcs can be found by a maximum flow computation. This circulation gives an augmentation of the current generalized flow g. After such an augmentation, all gain cycles in 66 CHAPTER 6. THE GENERALIZED FLOW PROBLEM Step 1 Find pu, a canonical labeling from the source. If - (v,w) < 1 on every arc of the residual graph of the current generalized flow g, then STOP (the current flow is optimal). Step 2 Let a = max(V,w)EE, 7Y(v, w). Consider the network G' consisting of all arcs (v, w) with yp(v,w) = 1, and all arcs in A = {(v,s)I7u(v,s) = a} and their opposites. Find a (standard) circulation f' in G' that maximizes 2 "(.,,)EA f'(VS), the flow into s. Step 3 Let g' be the generalized flow corresponding to f', i.e, g'(v, w) = f'(v, w) if v 0s and g'(s, v) = --y,(v, s)f'(v, s). Update the current solution by setting g(v, w) = g(v, w) + g'(v, w)p(v) V(v, w) E E, Figure 6.1: The Inner loop of Truemper's algorithm. the residual graph pass through the source, and the maximum gain of a flow-generating cycle in the residual graph is decreased. Therefore, this algorithm runs in finite time. Truemper's algorithm is in some sense an analog of Jewell's [61] minimum-cost flow algorithm. which augments the flow along all of the ch=,cst .- gmenting paths at once using a maximum flow subroutine. Using the same network as Zadeh [105] used for the minimum-cost flow problem. Ruhe [90] gave an example on which Truemper's algorithm takes exponential time. Goldberg, Plotkin, and Tardos [42] gave two algorithms, one that uses a minimum-cost flow computation as a subroutine, and one that builds on ideas from several maximum and minimumcost flow algorithms. In the next section we describe the first algorithm in detail and very briefly outline the second one. 6.3 Polynomial-Time Combinatorial Algorithms The main idea of the first polynomial-time algorithm of Goldberg-Plot kin-Tardos is best described by contrasting the algorithm with Truemper's. In each iteration, both algorithms solve a simpler flow problem, and interpret the result as an augmentation in the generalized flow network. Truemper's algorithm is slow because at each iteration it considers only arcs with unit relabeled gain and some of the arcs adjacent to the source, disregarding the rest of the graph completely. The Goldberg-Plotkin-Tardos algorithm, which we shall call algorithm MCF, considers all arcs. It assigns a cost c(v, w) = - log -y,(v, w) to each arc (v, w) and sol,, the resulting minimum-cost circulation problem (disregarding gains). The interpretationof a pseudoflow f is a generalized pseudoflow g, such that g(v, w) = f(v, w) if f(v, w) >_0 and g(v, w) = --y,(w, v)f(w, v) otherwise. In Truemper's algorithm, the interpretation of a feasible circulation on the subnetwork G' is a feasible generalized flow on the original network. In Algorithm MCF, however, the interpretation of a minimum-cost circulation is a generalized pseudoflo":. Ar!: of the flow that have a relabeled gain of less than 1 produce vertices with deficits in the interpretation. A connection between a pseudoflow f and its interpretation is given by the 6.3. POLYNOMIAL-TIME COMBINATORIAL ALGotk1THMS 67 Step 1 Find j, a canonical labeling of the residual graph of the generalized pseudoflow g. If t,(v,w) < 1 on every arc of the re,"dual graph and Vv E (V - {8}) : e,,, (v) = 0, then HALT (the current flow is optimal). Step 2 Introduce costs c(v,w) = -log7-(v,w) on the arcs of the network. Solve the capacitated trawss.hipment problem in the residual relabeled network of G with demands d(v) = e,M,(v) for v E V - {s}. Step 3 Let 9' be the interpret.t;i-.n of f'. Update the current solution by setting g(v, w) = g(v, u) + g'(v, w)p (v) V(v, w) E E. Figure 6.2: Inner loop of Algorithm MCF. following lemma. Lemma 6.3.1 The residual graphs of a pseudoflow f and its interpretation g as a generalized pseudoflow are the same. The firs' iteration of the algorithm solves a minimum-cost circulation problem, and it creates excess at the source and deficits at various other vertices. Each subsequent iteration uses some of the excess at the source to balance the deficits created by the previous iterations by solving a capacitated transshipment problem. More precisely, Algorithm MCF, shown in Figure 6.2, maintains a generalized pseudoflow g in the original (non-relabeled) network, such that the excess at every vertex other than the source is non-positive. The algorithm proc-cds in iterations. At ei iteration it canonically relabels the residual graph, solves the corresponding capacitated transshipment problem in thc :elabdled network, and interprets the result as a generalized augmentation. The most important property of a minimum-cost pseudoflow, for this application, is that its residual graph has no negative cycles. Th"s implies (by Lemma 6.3.1) that the residual graphs of the generalized pseudoflows produced by the algorithm in Step 3 have no flow-generating cycles. The following lemma (analogous to Lemma 1.6.1) implies that this is enovgh to ensure that the new residual graph can be canonically relabeled. Lemma 6.3.2 Let g be a generalized pseudoflow in a restricted network. If the residua! graph of GQ contains no flow-generating cycles, and the excess at every vertex other than s is non-positive, then every vertex is reachable from s in G. The relabeled gain factors are at most 1 (except for arcs entering the source in the first iteratioi:): therefore the flow computation creates deficits, but no excesses at vertices other than the source. Using the decomposition of pseudoflows (Theorem 1.7.3), one can prove that in each relabeled network there exists a flow satisfying all the deficits. These properties are summarized in the following lemma. Lemma 6.3.3 For a generalized pseudoflow g that is constructed by Step 3, the following properties hold: The residual graph of g has no flow-generating cycle; canonical relabeling applies to the residual 68 CHAPTER 6. THE GENERALIZED FLOW PROBLEM graph of g; all excesses, except the excess of the source, are non-positive; the pseudoflow required in Step 2 of the next iteration exists. The algorithm terminates once a generalized flow has been found. Throughout the computation the residual graph contains no flow-generating cycles; therefore (by Theorem 1.6.2) if the algorithm ever finds a generalized flow, then this fiow is optimal. Next we bound the number of iterations of the algorithm. Consider the generalized pseudoflow g existing at the beginning of an iteration. Because there are deficits but no flow-generating cycles in the residual graph, the current excess at the source is an overestimate of the value of the maximum possible excess. It is easy to see that, after a canonical relabeling, the sum of the deficits at all the vertices other than the source is a lower bound on the amount of the Uverestimate. We use this value, DeI(g,p) = Zt,g(-eg,(r)), as a measure of the proximity of a generalized pseudoflow to an optimal generalized flow. It is not too hard to show that if Def(g,lt) is smaller than L then the algorithm terminates in one more iteration. (Recall that L is the product of denominators of arc gains.) Theorem 6.3.4 [42] If Def(g,p) < L before Step 2 of an iteration, then the algorithm produces a generalized flow in Step 3. An important observation is that the labels u are monotonically decreasing during the algorithm. The decrease in the labels is related to the price function associated with the minimum-cost pseudofiow computation. Let p' denote the optimal price function associated with the pseudoflow f found in Step 2. Assume, without loss of generality, that p'(s) = 0. Lenima 6.3.5 For each vertex v, the canonical relabeling in Step 1 of the next iteration decreases the label ji(v) by at least a factor of 2-'.' The key idea of the analysis is to distinguish two cases: Case 1, in which the pseudoflow f' is along "cheap" paths, (i.e., p'(v) is small, say p'f v) < log 1.5, for every v,E V); and Case 2, in which there exists a vertex v such that p'(v) is "large" (_ log 1.5). In the first case, Def(fp) decreases significantly; in the second case, Lemma 6.3.5 implies that at least one of the labels decreases significantly. The label of a vertex v is the gain of the current most efficient path from s to r. This limits the number of times Case 2 can occur. Theorem 6.3.4 can be used to guarantee that Case 1 cannot occur too often. The following lemma provides a tool for estimating the total deficit created when interpreting a minimum-cost pseudoflow as a generalized pseudoflow. Lemma 6.3.6 Let f' be a pseudoflow along a simple path P from s to some other vertex r that satisfies one unit of deficit at v. Let g' be the interpretation of f' as a generalized pseudoflow. Assume that all relabeled gains along the path p are at most 1, and denote them by 11. ...% Then after 6.3. POLYNOMIAL-TIME COMBINATORIAL ALGORITHMS 69 augmenting by g', the unit of deficit at v is replaced by deficits on vertices along P that sum up to at most (Ili<k< 71)' - 1. Proof: The deficit created at the ith vertex of the path is (1 - ji), for i = 1,... ,k. Using the assumption that the gain factors along the path are at most 1, the sum of the deficits can be bounded by k E(-Tly/ k - ... 1 1. This lemma can be used to bound the value of Def(gp) after an application of Step 3. Let p' be an optimal price function associated with the pseudoflow f', such that p'(s) = 0. Let 3 = max .EV p'(v) be the maximum price. Using the pseudoflow decomposition theorem (Theorem 1.7.2). one can show that a minimum-cost pseudoflow f' can be decomposed into flows along paths from the source s to the other vertices and cycles in the residual graph of g, such that the cost ol the cycles is zero and the cost of each path ending at a vertex v is at most p'(v). To bound the amount of deficit created by the interpretation, we shall consider the interpretation of the flow on the cycles and the paths one-by-one. Since the cycles consists of zero-cost arcs only, the interpretation of the flow along these cycles does not create deficits. When interpreting flow along a path from s to 7'. the deficit at v is replaced by deficits that sum to at most 2P'(v) -- I < 2- - 1 times the deficit at v satisfied by this path. This proves the following lemma. Lemma 6.3.7 The value of Def(g,y) after an application of Step 3 can be bounded by 23 - 1 times its value before the step. In particular, if p'(v) < log 1.5 for every vertex v, then Def(g. p) decreases by a factor of 2. The remair.ng difficulty is the fact that the function Def(g,p) can increase, both when Case 2 applies in Step 3 and due to the relabeling in Step 1 (changing the currency unit from £ to penny increases a 5 £ deficit to 500 pennies). The increase in Def(g,p) can be related to the decrease in the labels, however. The deficit at a vertex increases in Step 1 by a factor of a if and only if the label of this vertex decreases by the same factor. The increase in Def(gq) during Step 3 is bounded by 213, where 3 = maxp'(v), by Lemma 6.3.7. By Lemma 6.3.5, this means that p(c) for some vertex v dec'reases by at least 03 luring Step 1 of the next iteration. lience, in both cases, if Def(gq) increases by a during a step, then there exists a vertex v for which p(v) decreases by at least q during the next execution of Step 1. Let T < B' denote an upper bound on the gain of any simple path. The iabel of a vertex v is the gain of a simple path from s to v in the residual graph. Therefore the labels are at least T - 1 and at most T. and the label of a vertex cannot decrease by more than a factor of 72 during the algorithm. This gives the following lemma. CHAPTER 6. THE GENERALIZED FLOW PROBLEM 70 Lemma 6.3.8 The growth of the function Def(g,Iu) throughout an execution of the algorithm is bounded by a factor of T ( ) n = BO(n 2 ) , Theorem 6.3.9 (42] The algorithm terminates in O(n 2 log B) iterations. Proof: The label of any vertex v is the gain of a simple path from s to v, and it is monotonicallv decreasing during the algorithm. Therefore, it cannot decrease by a constant factor more than O(logT) times for any given vertex, and Case 2 cannot occur more than O(nlogT) times. When Case 1 applies, the value of Def(g,p) decreases by a factor of 2. The value of Def(g, 1) is at most O(nBT) after the first iteration; and, by Theorem 6.3.4, the algorithm terminates when this value decreases below L 1 . Lemma 6.3.8 limits the increase of Def(g,p) to TOM during the algorithm. Hence, Case 1 cannot occur more than O(log(nBT) + nlogT +]ogL) = O(n 2 logB) times. I To get a bound on the running time, one has to decide which algorithm to use for computing the minimum-cost pseudoflow in Step 2. The bcst choice is Orlin's O(,m(m + -!ogn)log n) time algorithm, discussed in Chapter 5. Theorem 6.3.10 [42] Algorithm MCF solves the generalized flow problem using at most O(n 2 r(m + n log n)log n log B) arithmetic operations on numbers whose size is bounded by O(rn log B). We conclude this section with a brief discussion of the other algorithm of [42]. The algorithm is based on ideas from two flow algorithms: the cycle-canceling algorithm of Goldberg and Tarjan [45] described in Chapter 4, and the fat-path maximum flow algorithm of Edmonds aird Karp [21], which finds the maximum flow in a network by repeatedly augmenting along a highest-capacity residual path. Each iteration of the algorithm starts with cycle-canceling. Canceling a flow-generating cycle in a generalized flow network creates an excess at some vertex of the cycle. If we cancel cycles other than the most efficient ones, the residual graph of the resulting generalized pseudoflow will have flow-generating cycles that do not contain the source. The algorithm cancels flow-generating cycles by an adaptation of the Goldberg-Tarjan algorithm. Using the dynamic tree data structure, this phase can be implemented to run in O(n 2 rnlognlogB) time. The resulting excesses at vertices other than the source are moved to the source along augmenting paths in the second phase of the algorithm. Consider a generalized pseudoflow that has non-negative excesses and whose residual graph has no flow-generating cycles. The key idea of the second phase is to search for augmenting paths from vertices with excess to the source that result in a significant increase in the excess of the source. The algorithm maintains a scaling parameter A such that the maximum excess at the source is at most 21nA more than the current excess. It looks for an augmenting path that can increase the 71 6.3. POLYNOMIAL-TIME COMBINATORIAL ALGORITHMS excess of the source by at least A. If the residual graph of the current generalized pseudoflow has no flow-generating cycles, then one can find a sequence of such paths such that, after augmenting the generalized flow along these paths, the maximum excess at the source is at most mA moje than the current excess. This phase can be implemented in 0(m(m + nlogn)) time by a sequence of 0(m) single-source shortest path computations on graphs with arcs of non-negative cost. Now A is divided by two and a new iteration is started. After O(m log B) iterations, when the excess at the source is very close to the optimum value, Truemper's algorithm can be applied to bring all of the remaining excesses to the source. Maggs (personal communication, 1989) has observed that this algorithm can be improved through better balancing of the two phases. Dividing A by a factor of 2' after each iteration decreases the number of iterations by a factor of a and increases the time required for the second phase by a factor of 2'. The parameter a is chosen so that the time required for the two phases is balanced. The resulting algorithm runs in 0 (n2m2log2 Bog(( 2 1))+Iogo) time. 72 CHAPTER 6. THE GENERALIZED FLOW PROBLEM Bibliography [1] R. K. Ahuja, A. V. Goldberg, J. B. Orlin, and R. E. Tarjan. Finding Minimum-Cost Flows by Double Scaling. Technical Report CS-TR-164-88, Department of Computer Science, Princeton University, 1988. [2] R. K. Ahuja, K. Mehlhorn, J. B. Orlin, and R. E. Tarian. Faster Algorithms for the Shortest Path Problem. Technical Report CS-TR-154-88, Department of Computer Science, Princeton University, 1987. [3] R. K. Ahuja and J. B. Orlin. A Fast and Simple Algorithm for the Maximum Flow Problem. Sloan Working Paper 1905-87, Sloan School of Management, M.I.T., 1987. [4] R. K. Ahuja, J. B. Orlin, and R. E. Tarja!, Improved Time Bounds for the Maximum Flow Problem. Technical Report CS-TR-118-s , Department of Computer Science, Princeton University, 1987. (SIAM J. Comput., to appear). [5] B. Awerbuch. Personal Communication. Laboratory for Computer Science, M.I.T., 1985. [6] M. L. Balinski. Signature Methods for the Assignment Problem. Oper. Res., 33:527-536, 1985. [7] F. Barahona and E. Tardos. Note on Weintraub's Minimum Cost Circulation Algorithm. SIAM J. Comput., to appear. [8] C. Berge. Graphs and Hypergraphs. North Holland, 1973. [9] D. P. Bertsekas. A Distributed Algorithm for the Assignment Problem. Unpublished Manuscript, Lab. for Decision Systems, M.I.T., 1979. [10] D. P. Bertsekas. A Distributed Asynchronous Relaxation Algorithm for the Assignment Problem. In Proc. 24th IEEE Conference on Decision and Control, Ft. Lauderdale, FL. 1985. [11] D. P. Bertsekas. Distributed Asynchronous Relaxation Methods for Linear Network Flow Problems. Technical Report LIDS-P-1986, Lab. for Decision Systems, M.I.T., 1986. [12] D. P. Bertsekas and J. Eckstein. Dual Coordinate Step Methods for Linear Network Flow Problems. Math. Programming,42:203-243, 1988. [13] R. G. Bland and D. L. Jensen. On the Computational Behavior of a Polynomia]-Tini Network Flow Algorithm. Technical Report 661, School of Operations Research and Industrial Engineering, Cornell University, 1985. 73 74 BIBLIOGRAPHY [14] R. G. Busacker and P. J. Gowen. A Procedure for Determinimg a Family of Minimal-Cost Network Flow Patterns. Technical Report 15, O.R.O., 1961. [15] R. G. Busacker and T. L. Saaty. Finite Graphs and Networks: An Introduction with Applications. McGraw-Hill, New York, NY., 1965. [16] J. Cheriyan and S.N. Maheshwari. Analysis of a Preflow Push Algorithm for Maximum Netwrok Flow. Unpublished manuscript, Department of Computer Science and Engineering, Indian Institute of Technology, New Deli, India, 1987. [17] R. V. Cherkasky. Algorithm of Construction of Maximal Flow in Networks with Complexity of O(V 2'v'E) Operations. Mathematical Methods of Solution of Economical Problems, 7:112125, 1977. (In Russian). [18] W. H. Cunningham. A Network Simplex Method. Math. Programming, 11:105-116, 1976. [19] G. B. Dantzig. Linear Programming and Eztensions. Princeton Univ. Press, Princeton, NJ. 1962. [20] E. A. Dinic. Algorithm for Solution of a Problem of Maximum Flow in Networks with Power Estimation. Soviet Math. Dokl., 11:1277-1280, 1970. [21) J. Edmonds and R. M. Karp. Theoretical Improvements in Algorithmic Efficiency for Network Flow Problems. J. Assoc. Comput. Mach., 19:248-264, 1972. [22] G. M. Engel and H. Schneider. Diagonal Similarity and EquivaJence for Matrices Over Groups with 0. Czech. Math. J., 25:389-403, 1975. [23] S. Even. Graph Algorithms. Computer Science Press, Potomac, MID, 1979. [24] S. Even and R. E. Tarjan. Network Flow and Testing Graph Connectivity. SLAM J. Comput., 4:507-518, 1975. [25] L. R. Ford, Jr. and D. R. Fulkerson. Maximal Flow Through a Network. Canadian Journal of Math., 8:399-404, 1956. [26] L. R. Ford, Jr. and D. R. Fulkerson. Flows in Networks. Princeton Univ. Press, Princeton, NJ., 1962. [27] S. Fortune and J. Wyllie. Parallelism in Random Access Machines. In Proc. 10th ACM Symp. on Theory of Computing, pages 114-118, 1978. [28] G. Frederickson. Fast Algorithms for Shortest Paths in Planar Graphs, with Applications. SIAM J. Comput., 16:1004-11q22, 1987. [29] M. L. Fredman and R. E. Tarjan. Fibonacci Heaps and Their Uses in Improved Network Optimization Algorithms. J. Assoc. Comput. Mach., 34:596-615, 1987. [30] S. Fujishige. A Capacity-rou -ding Algorithm for the Minimum-Cost Circulation Problem: a Dual Framework of the Tardos Algorithm. Math. Prog., 35:298-308, 1986. [31] D. R. Fulkerson. An Out-of-Kilter Method for Minimal Cost Flow Problems. SL4M J. Appl. Math, 9:18-27, 1961. BIBLIOGRAPHY 75 [32] H. N. Gabow. Scaling Algorithms for Network Problems. J. of Comp. and Sys. Sci., 31:148168, 1985. [331 H. N. Gabow and R. E. Tarjan. Almost-Optimal Speed-ups of Algorithms for Matching and Related Problems. In Proc. 20th ACM Symp. on Theory of Computing, pages 514-527, 1988. [34] H. N. Gabow and R. E. Tarjan. Faster Scaling Algorithms for Network Problems. SIAM J. Comput., to appear. [35] Z. Galil. An O(V 5 / 3 E 2/ 3 ) Algorithm for the Maximal Flow Problem. Acta Informatica, 14:221-242, 1980. 136] Z. Galil and A. Naamad. An O(EV log: V) Algurithm for the Maximal Flow Problem. J. Comput. System Sci., 21:203-217, 1980. [37] Z. Gil and L. Tardos. An O(n2 (m + nlogn)logn) Minimum Cost Flow Algorithm. J. Assoc. Comput. Maci., 35:374-386, 1988. [3S] F. Glover and D. Klingman. On the Equivalence of Some Generaliz, d Network Problems to Pure Network Problems. Math. Programming,4:269-278, 1973. [39] A. V. Goldberg. A New Max-Flow Algorithm. Technical Report MIT/LCS/TM-291, Laboratory for Computer Science, M.I.T., 1985. [40] A. V. Goldberg. Efficient Graph Algorithms for Sequential and Parallcl Computers. PhD thesis, M.I.T., January 1987. (Also available as Technical Report TR-374, Lab. for Computer Science, M.I.T., 1987). [41] A. V. Goldberg, M. D. Grigoriadis, and R. E. Tarjan. Efficiency of the Network Simple). Algorithm for the Maximum Flow Problem. Technical Report LCSR-TR-117, Laboratory for Computer Science Research, Department of Computer Science, Rutgers University, 1988. [42] A. V. Goldberg, S. A. Plotkin, and E. Tardos. Combinatorial Algorithms for the Generalized Circulation Problem. Technical Report STAN-CS-88-1209, Stanford University, 1988. (Also available as Technical Memorandum MIT/LCS/TM-358, Laboratory for Computer Science, M.I.T., 1988.). [43] A. V. Goldberg, S. A. Plotkin, and P. M. Vaidya. Sublinear-Time Parallel Algorithms for Matching and Related Problems. In Proc. 29th IEEE Symp. on Found. of Comp. Sci., pages 174-185, 1988. [44] A. V. Goldberg and R. E. Tarjan. A New Approach to the Maximum Flow Problem. In Proc. 18th ACM Symp. on Theory of Computing, pages 136-146, 1986. [45] A. V. Goldberg and R. E. Tarjan. Finding Minimum-Cost Circulations by Canceling Negative Cycles. Technical Report MIT/LCS/TM-334, Laboratory for Computer Science, M.I.T.. 1987. (Also available as Technical Report CS-TR 107-87, Departmcnt of Computer Science, Princeton University). 76 BIBLIOGRAPHY [46] A. V. Goldberg and R. E. Tarjan. Finding Minimum-Cost Circulations by Successive Approximation. Technical Report MIT/LCS/TM-333, Laboratory for Computer Science, M.I.T., 1987. (Math. of Oper. Res., to appear). [47] A. V. Goldberg and R. E. Tarjan. Solving Minimum-Cost Flow Problems by Successive Approximation. In Proc. 19th ACAl Symp. on Theory of Computing, pages 7-18, 1987. [48] A. V. Goldberg and R. E. Tarjan. A New Approach to the Maximum Flow Problem. J. Assoc. Comput. Mach., 35:921-940, 1988. [49] A. V. Goldberg and R. E. Tarjan. A Parallel Algorithm for Finding a Blocking Flow in an Acyclic Network. Technical Report STAN-CS-88-1228, Department of Computer Science, Stanford University, 1988. [50] D. Goldfarb and M. D. Grigoriadis. A Computational Comparison of the Dinic and Network Simplex Methods for Maximum Flow. Annals of Oper. Res., 13:83-123, 1988. [51] D. Goldfarb and J. Hao. A Primal Simplex Algorithm that Solves the Maximum Flow Problem in at Most nm Pivots and O(n 2 rn) Time. Technical report, Department of IE and OR, Columbia University, 1988. [52) M. Gondran and M. Minoux. Graphs and Algorithms. Wiley, 1984. [53] M. D. Grigoriadis. An Efficient Implementation of the Network Simplex Method. Proy. Study, 26:83-111, 1986. Math. [54] M. Gr6tschel, L. Lovisz, and A. Schrijver. Geo~aetrc Algorithms and Combinatorial Optimization. Springer Verlag, 1988. [55] R. Hassin. Maximum Flows in (s,t) Planar Networks. Information Processing Let., 13:107108, 1981. [56] R. Hassin and D. B. Johnson. An O(nlog 2 n) Algorithm for Maximum Flow in Undirected Planar Network. SIAM J. Comput., 14:612-624, 1985. [57] J. E. Ilopcroft and R. M. Karp. An n5/2 Algorithm for Maximum Matching in Bipartite Graphs. SIAM J. Comput., 2:225-231, 1973. [58] M. Iri. A New Method of Solving Transportation-Network Problems. J. Op. Res. Soc. Japan. 3:27-87, 1960. [59] A. Itai and Y. Shiloach. Maximum Flow in Planar Network. SIAM J. Comput., 8:135-150. 1979. [60] P. A. Jensen and J. W. Barnes. Network Flow Programming. J. Wiley & Sons, 1980. [61] W. S. Jewell. Optimal Flow through Networks. Technical Report 8, M.I.T., 1958. [62] W. S. Jewell. Optimal Flow through Networks with Gains. Oper. Res., 10:476-499, 1962. [63] D. B. Johnson and S. Venkatesan. Using Divide and Conquer to Find Flows in Directed Planar Netwroks in O(n'- log n) Time. In Proc. 20th Annual Allerton Conf. on Communication, Control, and Computing, pages 898-905, 1982. BIBLIOGRAPHY 77 [641 S. Kapoor and P. M. Vaidya. Fast Algorithms for Convex Quadratic Programming and Multicommodity Flows. In Proc. 18th ACM Symp. on Theory of Computing, pages 147-159, 1986. [65] N. Karmarkar. A New Polynomial-Time Algorithm for Linear Programming. Combinatorica, 4:373-395, 1984. [66] R. M. Karp. A Characterization of the Minimum Cycle Mean in a Digraph. Discrete Math., 23:309-311, 1978. [67] R. M. Karp, E. Upfal, and A. Wigderson. Constructing a Maximum Matching is in Random NC. Combinatorica,6:35-48, 1986. [68] A. V. Karzanov. Determining the Maximal Flow in a Network by the Method of Preflows. Soviet Math. Dok., 15:434-437, 1974. [69] L. G. Khachian. Polynomial Algorithms in Linear Programming. MAfatematiki i Matematicheskoi Fiziki, 20:53-72, 1980. Zhurnal Vychishtflnoi [70] M. Klein. A Primal Method for Minimal Cost Flows with Applications to the Assignment and Transportation Problems. Management Science, 14:205-220, 1967. [71] D. K6nig. Theorie der Endlichen und Unendlichen Graphen. Chelsea Publishing Co., New York, 1950. [72] H. W. Kuhn. The Htungarian Method for the Assignment Problem. Naval Res. Logist. Quart., 2:83-97, 1955. [73] E. L. Lawler. Combinatorial Optimization: Networks and Matroids. ltolt. Reinhart, and Winston, New York, NY., 1976. [74] V. M. Malhotra, M. Pramodh Kumar, and S. N. Maheshwari. An O(1I1"I) Algorithm for Finding Maximum Flows in Networks. Information Processing Let., 7:277-278, 1978. [75] G. L. Miller and J. Naor. Flow in Planar Graphs with Multiple Sources and Sinks. Unpublished Manuscript, Computer Science Department, Stanford University, Stanford, CA. 1989. [761 G. J. Minty. Monotone Networks. Proc. Roy. Soc. London, A(257):194-212, 1960. [77] K. Mulmuley, U. V. Vazirani, and V. V. Vazirani. Matching is as easy as matrix inversion. Combinatorica,pages 105-131, 1987. [781 K. Onaga. Dynamic Programming of Optimum Flows in Lossy Communication Nets. IEEE Trans. Circuit Theory, 13:308-327, 1966. [79] K. Onaga. Optimal Flows in General Communication Networks. J. Franklin Inst., 283:308327, 1967. [80] J. B. Orlin. Genuinely Polynomial Simplex and Non-Simplex Algorithms for the Minimum Cost Flow Problem. Technical Report No. 1615-84, Sloan School of Management, MIT. 1984. [81] J.B. Orlin. On the Simplex Algorithm for Networks and Generalized Networks. Math. Prog. Studies, 24:166-178, 1985. 78 BIBLIOGRAPHY 182] J. B. Orlin. A Faster Strongly Polynomial Minimum Cost Flow Algorithm. In Proc. 20th ACM Symp. on Theory of Computing, pages 377-387, 1988. [83 J. B. Orlin and R. K. Ahuja. New Distance-Directed Algorithms for Maximum-Flow and Parametric Maximum-Flow Problems. Sloan Working Paper 1908-87, Sloan School of Management, M.I.T., 1987. [84] J. B. Orlin and R. K. Ahuja. New Scaling Algorithms for Assignment and Minimum Cycle Mean Problems. Sloan Working Paper 2019-88, Sloan School of Management, M.I.T., 1988. [85] P. Van Emde Boas and R. Kaas and E. Zijlstra. Design and implementation of an efficient priority queue. Math. Systems Theory, 10:99-127, 1977. [86) C. H. Papadimitriou and K. Steiglitz. Combinatorial Optimization: Algorithms and Cornplexity. Prentice-Hall, Englewood Cliffs, NJ, 1982. [87] J. 1H. Reif. Minimum s-t Cut of a Planar Undirected Network in O(nlog 2 n) Time. SIAM J. Comput., 12:71-81, 1983. [88] J. Renegar. A Polynomial Time Algorithm, Based on Newton's Method, for Linear Programming. Mathematical Programming,40:59-94, 1988, [89] H. Rbck. Scaling Techniques for Minimal Cost Network Flows. In U. Pape, editor, Discrctc Structures and Algorithms, pages 181-191. Carl Hansen, Mfinich, 1980. [90] G. Ruhe. Parametric Maximal Flows in Generalized Networks - Complexity -A Algorithms. Optimization, 19:235-251, 1988. [91] Y. Shiloach. An O(nIlog2 I) Maximum-Flow Algorithm. Technical Report STAN-CS-78-802, Computer Science Department, Stanford University, Stanford, CA, 1978. [92] Y. Shiloach and U. Vishkin An O(n 2 logn) Parallel Max-Flow Algorithm. J. Algorithms, 3:128-146, 1982. [93] D. D. Sleator. An O(nmlogn) Algorithm for Maximum Network Flow. Technical Report STAN-CS-80-831, Computer Science Depaxtment, Stanford University, Stanford, CA, 1980. [94] D. D. Sleator and R. E. Tarjan. A Data Structure for Dynamic Trees. J. Comput. Systcm Sci., 26:362-391, 1983. [95] E. Tardos. A Strongly Polynomial Minimum Cost Circulation Algorithm. 5(3):247-255, 1985. Combinatorica, [96] R. E. Tarjan. Data Structures and Network Algorithms. Society for Industrial and Applied Mathematics, Philadelphia, PA, 1983. [97] R. E. Tarjan. A Simple Version of Karzanov's Blocking Flow Algorithm. OperationsResearch Letters, 2:265-268, 1984. [98] R. E. Tarjan. Amortized Computational Complexity. SIAM J. Alg. Disc. Math., 6:306-31R. 1985. [99] R. E. Tarjan. Efficiency of the Primal Network Simplex Algorihm for the Minimum-Cost Circulation Problem. Technical Report CS-TR-187-88, Department of Computer Science, Princeton University. 1988. BIBLIOGRAPHY 79 [100) N. Tomizawa. On Some Techniques Useful for Solution of Transportation Networks Problems. Networks, 1:173-194, 1972. [1011 K. Truemper. On Max Flows with Gains and Pure Min-Cost Flows. SIAM J. Appl. Math., 32:450-456, 1977. [102] P. M. Vaidya. Speeding up Linear Programming Using Fast Matrix Multiplication. Technical Memorandum, AT&T Bell Laboratories, Murray Hill, NJ, 1989. [103] A. Weintraub. A Primal Algorithm to Solve Network Flow Problems with Convex Costs. Management Science, 21:87-97, 1974. [104] M. A. Yakovleva. A Problem on Minimum Transportation Cost. In V. S. Nemchinov, editor, Applications of Mathematics in Economic Research, pages 390-399. Izdat. Social'no-Ekon. Lit., Moscow, 1959. [105] N. Zadeh. A Bad Network Problem for the Simplex Method and Other Minimum Cost Flow Problems. Mathematical Programming, 5:255-266, 1973. i 4! i