A Fast Approximation Scheme for Fractional Covering Problems with Variable Upper Bounds
Lisa Fleischer
Abstract
We present the first combinatorial approximation scheme for mixed positive packing and covering linear programs that yields a pure approximation guarantee. Our algorithm returns solutions that simultaneously satisfy general positive covering constraints and packing constraints that are variable upper bounds. The returned solution has positive linear ob jective function value at most 1 + times the optimal value. Our approximation scheme is based on Lagrangianrelaxation methods. Previous such approximation schemes for mixed packing and covering problems does not simultaneously satisfy packing and covering constraints exactly. We show how to exactly satisfy general positive covering constraints simultaneously with variable upper bounds. A natural set of problems that our work addresses are linear programs for various network design problems: generalized Steiner network, vertex connectivity, directed connectivity, capacitated network design, group Steiner forest. These are all NP-hard problems for which there are approximation algorithms that round the solution to the corresponding linear program. Solving the linear program is often the computational bottleneck in these problems, and thus a fast approximation scheme for the LP relaxation means faster approximation algorithms. For the special case of survivable network design, we introduce a new modification of the push-relabel maximum flow algorithm that allows us to perform each iteration in amortized O(m + n log n) time, instead of one maximum flow per iteration that is implied by the straight forward adaptation of our general algorithm. (m is the number of edges and n is the number of vertices in the network.) In conjunction with an observation that reduces the number of iterations to m log n for {0, 1} constraint matrices, the modification allows us to obtain an algorithm that is faster than existing exact or approximate algorithms by a factor of at least O(m) and by a factor of O(m log n) if the number of demand pairs is (n).

integer vectors b, c, and u. It has the form (1.1) (1.2) (1.3) (1.4) min cT x s.t. Ax  b Px  u x0

1

Intro duction Problem Statement. A mixed positive packing and covering program (MPCP) is a linear program described by nonnegative, integer n × m matrix A; nonnegative, integer n × m matrix P ; and nonnegative
Scho ol of Industrial Administration, Carnegie Mellon University, Pittsburgh, PA 15213, USA, Email: lkf@andrew.cmu.edu. Supp orted in part through NSF CAREER Award CCR-0049071, NSF award EIA-0049084, and the IBM T. J. Watson Research Center.
 Graduate

The problem without (1.3) is a fractional covering problem. The problem without (1.2) and with (1.1) changed to max cT x is a fractional packing problem. Mixed packing and covering programs can be solved exactly using standard LP solvers. However, these solvers are often slow for large problems. In response to this, approximation schemes for MPCP's were proposed in the early 1990's [29, 19, 20], generalizing ideas that appeared earlier in approximation schemes for fractional multicommodity flow problems [26, 27, 31]. Since this initial work, there has been substantial progress made on refining the algorithms. The result has been substantial reduction in complexity of both the run time and the analysis [10, 14, 24, 25, 30, 33, 34]. Empirically, approximate LP algorithms have proven to be effective in obtaining "warm starts" for exact LP solvers, and have thus enabled significant reductions in run times [1, 2, 3, 16] even when exact solutions are desired. All of the above-mentioned approximation schemes that consider mixed positive packing and covering problems find solutions that are only approximately feasible: these schemes find a solution of a fixed value that either satisfies Ax  b and P x  (1 + )u or satisfies Ax  b/(1 + ) and P x  u. Since x does not satisfy all constraints exactly, the value of x can be arbitrarily far from the value of a true solution. Our Contribution. In Section 2, we present the first Lagrangian based approximation scheme that simultaneously satisfies arbitrary covering constraints and packing constraints that are variable upper bounds. For the general positive covering problem with variable upper bounds, we present the first approximation scheme that finds exactly feasible solutions with value at most 1 + times the optimum. Our algorithm uses O( -2 m log(C m)) calls to an oracle that finds a most

Copyright © 2004 by the Association for Computing Machinery, Inc. and the Society for industrial and Applied Mathematics. All Rights reserved. Printed in The United States of America. No part of this book may be reproduced, stored, or transmitted in any manner without the written permission of the publisher. For information, write to the Association for Computing Machinery, 1515 Broadway, New York, NY 10036 and the Society for Industrial and Applied Mathematics, 3600 University City Science Center, Philadelphia, PA 19104-2688

1001


| violated constraint. (Here C := min|jcc||0 cj  ||c|| .) : j> Thus, our algorithm can approximately solve exponentially sized linear programs that are described implicitly by a polynomial time oracle. A natural set of problems that our work addresses are linear programs for various network design problems. For these problems, a provably good integer solution can be found by rounding the solution to an LP relaxation. Examples include generalized Steiner network [23] (also called survivable network design), element connectivity [12], vertex connectivity [5, 11], directed connectivity [28], capacitated network design [4], and group Steiner forest [15]. Our algorithm builds on the approximation schemes of Garg and K¨nemann [14] and Fleischer [10]. One o key innovation of our algorithm is to use scaling phases to control the increase in variable values. At the end of the -phase, the algorithm has a solution x that satisfies Ax  b and x  u. Thus, the solution x/ is a feasible solution to MPCP. Throughout the -phase, the solution x satisfies Ax  b/(1 + ) and x  u. Our algorithm works by using an oracle that given the current iterate x, finds a constraint (A(i), b(i)) such that A(i)x < b(i). It then increases x to increase the value of the left hand side of this constraint. This is iterated. One key lemma is that if the initial problem is feasible and A(i)x < b(i), then there is some index j such that x(j ) < u(j ) and A(i, j ) > 0. Thus our algorithm never gets stuck. In addition, with each increase, there is some index j for which we either add a substantial amount to x(j ), or x(j ) hits its phase upper bound of u(j ). We show that after at most O( 1 log(mC )) phases and O( 1 m log(mC )) iterations, 2 2 we encounter a feasible, -optimal1 solution and a proof of its approximate optimality. The second innovation is a new method of using dual updates to decrease the value of a dual solution to compensate for not increasing primal variables that are currently at their upper bound. The dual LP of our problem is not a mixed packing and covering LP, and has negative coefficients in the ob jective function. Thus, this requires a new method of handling dual variables. In Section 3 we show how to remove dependency of number of phases and iterations on log C when A is a {0, 1} matrix, by replacing log C m with log m . For the special case of survivable network design (also called generalized Steiner problem), we show how to implement all the iterations that occur in one phase

in O(mn min{n, k } log n) time, plus linear work per iteration, where m is the number of edges, n is the number of vertices in the network, and k is the number of pairs of vertices with positive connectivity requirement. Without modification to our general algorithm, the work per iteration is min{n, k } maximum flow computations. This can be reduced to one maximum flow per iteration using ideas from [10]. In Section 4, we show how to reduce this further to amortized m + n log n time per iteration. Our key insight involves a new modification to the push-relabel maximum flow algorithm to show that the work per phase can be performed in the same asymptotic time as one maximum flow computation. Thus we obtain an algorithm which is faster than existing algorithms by a factor of at least O(min{m, nk }), and faster by a factor of O(m log n) if the number of demand pairs is (n). Related Work. The symmetric problem, the positive packing problem with lower bounds on variables, is easily transformed into a pure positive packing problem by a change of variables. With an oracle that given x finds a constraint P x/u  , an -optimal solution can be found using earlier work. For instance, the algorithm of Young [34] solves this problem by solving the dual pure covering problem, where the (perhaps exponential number of ) variables are handled implicitly via the oracle. Independently and concurrently with our own work, Garg and Khandekar [13] develop a very different approximation scheme for linear program relaxations of problems on clutters. These yield exact polynomial time approximation schemes for problems when A is restricted to be a 0-1 matrix. This captures both the survivable network design problem, and the group Steiner problem, but not capacitated network design problems. Their algorithm does not run in polynomial time for general matrices A. In addition, their algorithm specialized to the case of SND requires O(m min{n, k } log2 k ) iterations, where each iteration requires a minimum cost flow computation. In contrast, for the special case of SND our algorithm requires only O(m log n) iterations, where the amortized time per iteration is O(n min{n, k } log n), and an additional min{n, k } minimum cost flow computations. 1.1 Preliminaries and Notation. For vector c, c(j ) denotes the j -th component of c. For a matrix A, A(i) denotes the i-th row of A and A(i, j ) denotes the j th component of the ith row of A. For a specified A, b, u, and c, OPT denotes the value of the optimal solution to (1.1)-(1.4) when P is the identity matrix. Denote by 1 the vector of all 1's. For two vectors v and w, v w and v T w both denote the dot product of v and

1 A feasible solution to a minimization problem is -optimal if its value is at most (1 + ) times the optimal value. A feasible solution to a maximization problem is -optimal if its value is at least the optimal value divided by (1 + ).

1002


w. Without loss of generality, we assume u = 1, since we can scale A and c correspondingly: multiply A(i, j ) and c(j ) by u(j ), and replace u(j ) by 1. We also assume that minj c(j ) > 0. If not, all variables in the set {j |c(j ) = 0} are included in the solution, and the problem is reduced by removing these variables and modifying A and b accordingly. Thus, we also assume OPT > 0, since otherwise the problem is trivially solved. Without loss of generality, we assume minj c(j ) = 1 by dividing c by the minimum (nonzero) component. 2 Fractional Covering with Variable Upp er Bounds In this section, we show how to obtain a solution that exactly satisfies (1.2)-(1.4) and has value in (1.1) at most (1+ )OPT, when P x  u are variable upper bounds and A may have an exponential number of rows and may be given implicitly by an approximate constraint oracle : If mini A(i)x/bi = , then an approximate constraint oracle returns a constraint p with A(p)x/bp  (1 + ). The general form of our problem is (P) minimize sub ject to P (x) := 1  i  n : 1  j  m : 1  j  m : The corresponding dual LP: (D) maximize D(y , z ) := sub ject to 1  j  m : 1  i  n : 1  j  m : in A(i, j )y (i) y (i) -z (j )  c(j ) 0 z (j )  0. in b(i)y (i) - jm z (j ) jm c(j )x(j ) A(i, j )x(j )  b(i) x(j )  1 x(j )  0
 (1+

An -phase starts with x  /(1 + ) and A(ii)x  b( ) ) for all 1  i  n. It iteratively picks a constraint p

=1 j

=1

=1

=1

2.1 The Approximation Scheme. Algorithm FractionalCoverUB(c, A, b, ) (appearing in Figure 1) operates in -phases, during which  is fixed. This algorithm differs from existing packing and covering schemes in that x   throughout the -phase. This restriction is necessary so that division by  yields 0 a feasible solution. We show that by proceeding in (y , z )  0. Denote by (x1 this initial value of x, and A i)   0 -phases, this restriction does not increase the number define (0) as C mini b(i) . Then (0)  C = x (j ), so 0 that x /(0) is feasible for (P). of iterations by too much.

with A(p)x <  (lines (2) and (16)) and increases A(p)x b(p) while maintaining x  . The phase ends when A( i ) x   for all 1  i  n. Thus, at the end of the b ( i) phase, x/ is a feasible primal solution. The next phase starts by setting  = (1 + ). To maintain x   while increasing A(p)x, FractionalCoverUB increases only variables in the set Q(p) = {1  j  m|A(p, j ) > 0 and x(j ) < } (line (6)). The amount of increase (specified in line (10)) is determined so that at least one variable x(j ) either increases by a (1 + ) factor or hits the phase upper bound of . If Q(p) is empty, we show below in Corollary 2.1 that this implies that all variables x(j ) with A(p, j ) > 0 must be set to 1 in every feasible solution. Thus, the subroutine Reduce (lines (7)-(8)) removes j from the problem and reduces b accordingly, by setting b(i) = b(i) - A(i, j ) for all constraints i. To prove -optimality of the final solution x , the algorithm also maintains dual solutions (y , z ). The algorithm returns primal and dual solutions that are feasible, and have ob jective function value within a (1+ ) factor of each other, proving their -optimality via LP duality theory. Each increase to the primal solution x is balanced with an increases to (y , z ) so that the ob jective function value of one can be related to the other. Since (P) is not a pure covering problem, the dual problem is not a pure packing problem. The algorithm incorporates a new dual update procedure to maintain this balance. These lines (13)-(15) and (20), are not necessary to find an -optimal solution to (P), but they are helpful to prove -optimality. Given (y , z ), let  be the minimum quantity such that (y / , z / ) is feasible for (D). Thus,  := i maxj [ A(i, j )y (i) - z (j )) ] / c(j ). In each iteration, one y -component is increased (line (14)). This increases the left hand side value of some dual constraints. Each time the left hand side of dual constraint j is increased by c(j ) (implying that if j determines  then the new  is increased by 1), this increase is balanced in (P) by multiplying x(j ) by at least (1 + ). Thus, if y (p) is increased and A(p, j ) > 0, but x(j ) is already at its upper bound, FractionalCoverUB increases z (j ) by A(p, j ) times the increase to y (p) (lines (14)-(15)) so that these two increases cancel each other and the left hand side of dual constraint j remains unchanged. |c | Initially x   /C where C = minj:|c(j|)0 c(j ) , and >

1003


The algorithm terminates when P (x)  1. (Denote this last phase by t.) The actual value of P is not important; what is important is that this implies that there have been a sufficient number of iterations, i.e. at least log1+ 1 = O( log(mC ) ) iterations.  2

Substituting p in for i in this last line, it follows that there is at least one index in the set {j | A(p, j ) > 0} that satisfies xt (j )  1. Since minj c(j )  1, we have that P (xt )  1.

We use the preceding lemmas to help bound the number of iterations. The following key lemma is useful 2.2 Algorithm Analysis. In this section, we prove to show that either the value of at least one variable x(j ) the following theorem, which is one of the main results strictly increases in each iteration, or we can reduce the of this paper: size of our problem by fixing variables and removing Theorem 2.1. Algorithm FractionalCoverUB returns constraints. feasible, -optimal solutions to (P) and (D), and re- Lemma 2.4. If (P) is feasible, then b(p)  j quires at most O( -2 m log(mC )) oracle cal ls. A(p, j ) inside the while loop of Let and denote the vectors x, y , and z at the end of the lth iteration in the k th phase. Let xk , y k , and z k be the corresponding vectors at the end of the k th phase. Let (k ) be the value of  in the k th phase. The first lemma shows that x stays bounded by  in the -phase. Lemma 2.1. Throughout FractionalCoverUB, x(j )   for al l 1  j  m. Proof. This is true at the start of the algorithm by the initial setting of x and . By definition of u in line (10) of FractionalCoverUB, after incrementing x(j ), we have that x(j )  x(j )(1 + c(j )/A(p,j ) )  x(j )(1 + u (c(j )/A(p,j ))(-x(j )) ) = x(j ) +  - x(j ) = . (c(j)/A(p,j))x(j)
k xk , yl , l k zl

FractionalCoverUB.

Q(p) /

Proof. Suppose to the contrary that b(p) < j A We have Q(P ) j (p, j ) inside the while lo op. / j that A(p, j )x(j )  = Q(p) A(p, j )x(j ) / j  A(p, j ) > b(p) > A(p)x, where the Q(p) / last inequality follows by the criteria to enter the while loop. This is a contradiction. Corollary 2.1. If (P) is feasible, then either Q(p) = j  or b(p) = A(p, j ). In the latter case, for every j with A(p, j ) > 0, variable x(j ) = 1 in every feasible solution. j Proof. If the problem is feasible, then b(i)  A(i, j ) for all constraij ts i. Thus, by Lemma 2.4, if Q(p) = , n then b(p) = A(p, j ) and every feasible solution must set x(j ) = 1 if A(p, j ) > 0.

In the next two lemmas, we establish a bound on the number of -phases by proving a significant increase Under the conditions of Corollary 2.1, we may in  between phases. reduce the problem size by removing all variables with A(p, j ) > 0 and updating b appropriately (i.e. For all i, Lemma 2.2. If there are t phases, then for al l 1  k < b(i)  b(i) - j :A(p,j )>0 A(i, j ).) This is the role of the t, (k + 1)  (1 + )(k). subroutine Reduce in line (8). Proof. All phases k except possibly the last phase t Lemma 2.5. The algorithm terminates after at most (1+ ) terminate when when min A(p)x/b(p)  (k ). At 2m log iterations. 1+ C  the start of phase k + 1, (k + 1) is set to (1 + Proof. Each time p is selected such that Q(p) = , ) min A(p)x/b(p)  (1 + )(k) in line (4). at least one variable is removed from the problem Lemma 2.3. After the completion of log1+ C phases, (Corollary 2.1). This can happen at most m times. In  P (x)  1. each iteration that Q(p) = , the increment amount is non-zero. With each increment, the variable j that T  (p) Proof. At start, (0) = C Ab(p) 1 for some constraint defines u in line (10) has an x-value that either increases T (p)T p satisfying Ab(p) 1 = mini A(i) ) 1 . At phase t := by a factor of 1 + , or hits . Variable x(j ) starts b( i equal to  /C , and never exceeds (1 + )/c(j), by the log1+ C , we have that  termination condition. Since x is nondecreasing, x(j ) can be multiplied by (1 + ) at most log1+ C (1+ ) times. A(p)T 1  = (0) C b(p)  Since  increases at most log1+ C times, the number   log1+ C / = (0)(1 + ) of times any variable can hit  is bounded by this same (t) by Lemma 2.2 quantity. Summing this over all 1  j  m yields the Tt for all i.  A(bi()i)x lemma.

1004


FractionalCoverUB(c, A, b, ) Input: cost vector c  Rm , n constraints Ax  b given by oracle. Output: -optimal primal and dual solutions x and (y , z ). (1) Initialize x   /C, (y , z )  0, x  x,   argmini A(i)x/b(i). (2) p  argmini A(i)x/b(i) (3) while P (x) < 1, (4)   (1 + )A(p)x/b(p) (5) while A(p)x/b(p) <  and P (x) < 1, (6) Q(p) = {1  j  m|A(p, j ) > 0 and x(j ) < }. (7) If Q(p) = , (8) Reduce(c, A, b, x, y , z ). (9) Else, ,j ))( (10) u = minj Q(p) min{c(j )/A(p, j ), (c(j )/A(px(j) -x(j )) } (11) For j  Q(p), (12) x(j )  x(j )(1 + c(j )/A(p,j ) ) u (13) yp  y p + u (14) For j  Q(p), / (15) z (j )  z (j ) + uA(p, j ) (16) p  argmini A(i)x/b(i) (17) end while (18) If cT x/ < cT x / , then x  x,   . (19) end while (20)   ((C (1 + ))1- m)-1/ . (21) return x / , ln 1 (y , z ).
m

­ oracle cal l ­ ­ -phase ­ ­ iteration ­

­ increment amount ­

­ oracle cal l ­ ­ keep best solution ­ ­ end -phase ­

Figure 1: An FPTAS for general positive covering problem with variable upper bounds. We prove in Lemma 2.7 the feasibility of the solution to (P) returned by FractionalCoverUB. For this we require the following technical lemma. Applying Lemma 2.6 with {gr }r = {xk (j )}k,l , {fr }r = l
i

{

Lemma 2.6. Given f0 = 0 and g0 > 0 and sequence {a1 , a2 , . . . , at } satisfying 0  ar  1 for al l r, suppose fr  fr-1 + ar and gr  gr-1 (1 + ar ). Then gr /g0  We are now ready to prove Theorem 2.1. The proof (1 + )fr for al l r. can be found in the Appendix. Proof. Initially g0 /g0 = 1 = (1 + )f0 . Since 1 + a  (1 + )a fr r all 0  a  1, the product r (1 + ar )  3 Zero-One Constraint Matrices o (1 + ) ar , when 0  ar  r 1 for all r. Thus, When A is a 0-1 matrix, and we can bound the value of OPT from above and below within a polynomial gr /g0  r=1 (1 + ai )  (1 + ) i=1 ai  (1 + )fr . i factor of accuracy, then we can replace the dependence Lemma 2.7. The final dual solution (y t , z t ) divided by in number of phases and iterations on log(C m) with log1+ C (1+ ) is a feasible solution for (D).  log( m ). This is done via two observations. Let x be a Proof. By our update of i (j ), y (p) and z (j ) in one iterx A(i,j )y (i)-z (j ) ation, when the ratio increases by an solution of value OPT. c(j ) 1) If c(j ) > OPT, then x (j ) = 0 (see Lemma 3.1 u additive term of c(j )/A(p,j )  1, the quantity x(j ) is below); u multiplied by (1 + c(j )/A(p,j ) ). By the termination cri- 2) if c(j ) < OPT/m, then setting x(j ) = 1 for all such teria of P (x) > 1, we have that xt (j )/x0 (j ) 
C (1+ ) . 

k k A(i,j )yl (i)-zl (j ) uk l }k,l , and {ar }r = { c(j )/A,(p,j ) }k,l , c(j ) i A(i,j )y k (i)-z k (j ) we have that is bounded by c(j ) (1+ ) log1+ C  for all k , which implies the lemma.

variables j increases the value of the solution by at most

1005


OPT, and hence choosing appropriately maintains our Lemma 3.1 also holds in the setting when bounded, approximation guarantee. multiple copies of edges are allowed. Garg and Of course, we can only use these observations if we can Khandekar [13] give a version of Lemma 3.1 when there bound OPT. are no upper bounds on the value of x(e). Theorem 3.1. If we have an estimate EST of OPT 4 Fractional Network Design such that The generalized Steiner network (GSN) problem is (3.5) OPT  EST   OPT defined by a graph G = (V , E ) with m = |E | and for  bounded by a polynomial in m, and al l entries n = |V |, an edge-cost vector c and a connectivity oracle of A are 0 or 1, then the number of phases and iter- f . Given a subset of vertices S  V , f (S ) specifies the ations required by FractionalCoverUB are O( -2 log m ) minimum number of edges with exactly one end point in S that must be included in any feasible solution. and O( -2 m log m ), respectively. The ob jective is to find a minimum cost set of edges Proof. Given EST, remove all variables j with c(j ) > satisfying f on all subsets S  V .  EST using Lemma 3.1. For all j with c(j ) < EST/m, A natural special case of this problem is the undiset x(j ) = 1. Then, after scaling the costs of the rected edge connectivity problem. We will refer to this remaining variables so that minj c(j ) = 1, we have that case as survivable network design (SND). In this probC := ||c||  m / . Thus, log(C m) = log m . lem, each pair (i, j ) of vertices has a connectivity re quirement r(i, j ) that is the number of (i, j ) edge disLemma 3.1. If c(j ) > OPT, then x (j ) = 0. joint paths required in any feasible solution. Then Proof. Suppose variable j0 has c(j0 ) > OPT and fSND (S ) = maxiS,j S r(i, j ). The integer version of / x (j0 ) = µ > 0. Let F be the set of variables with SND is NP-hard. Jain [23] describes an approximation strictly fractional x value, i.e. F := {j | 0 < x (j ) < 1}. algorithm that finds an integer solution with value at We modify x to obtain x as follows: most twice the optimal value. His algorithm iteratively ( rounds the solution to the linear program formulation, x j) = 0 if j = j0 and thus its efficiency depends on the efficiency of the  1 if x (j ) > 1 - µ LP solver.  x (j ) if x (j )  1 - µ. The linear program for GSN is given below. For 1-µ each edge e there is a variable x(e). In an integer i The cost of x s less than OPT: solution x, x(e) = 1 means that e is in the final solution  OPT - OPTµ  cx - c(j0 )µ and x(e) = 0 means that e is not in the solution. Define <  OPT. cx 1-µ 1-µ (S ) for S  V to be the set of edges with exactly one We now show that x is a feasible solution to (P), end point in S . Note that for the problem to be feasible contradicting the optimality of x . Consider constraint f () = f (V ) = 0, and f (S )  |(S )|. p with A(p, j0 ) = 1 and µp := b(p) - [A(p)T x - e x (j )] > 0. Define F> := {j  F - {j0 } : x (j ) > minimize E c(e)x(e) e  1 - µ} and F := {j  F - {j0 } : x (j )  1 - µ}. S  V : x(e)  f (S ) (S ) We show that the increase of values of variables in (GSN) e : x(e)  1 eitj er F> or F compensates for the loss of x (j ). If h j e : x(e)  0 [ 1 - A(p, j )x (j )]  µp , then A(p, j )x (j )  F> j  < b(p). Otherwise, F> [ 1 - A(p, j )x (j )] = µ We denote the optimal solution value to (GSN) µp j  µ. Since b(p) is integer and x is feasible, by OPT. To implement FractionalCoverUB to solve   1 - µ for some F A(p, j )x (j ) = N - µp + µ (GSN), we require an oracle that finds a constraint e positive integer N . This implies that the variables in (S ) x(e) < f (S ) if one exists. F are increased sufficiently to compensate for the loss Williamson [32] describes efficient oracles for gen of value x (j ): eral f . For SND, there is a well-known oracle for j fSND (S ) that uses a Gomory-Hu tree. We introduce A(p, j ) [ x (j ) - x (j ) ] = a new, ee cient procedure for finding inequalities of the ffi F j form 1 (S ) x(e) < fSND (S ) if they exist. Define k A(p, j ) x (j ) [ - 1 ]  µ. as the cardinality of the set {(i, j )  V × V | rij > 0}, 1-µ F so that k = O(n2 ). Naively using a Gomory-Hu tree, each search requires min{k , n} maximum flow compu-

1006


tations in a network with capacities determined by x. A simplification based on an idea introduced in [10] for multicommodity flows, cycles through each edge of the tree, and reduces the computation per search to one maximum flow. We show how to perform all searches in an -phase in the same asymptotic time as min{n, k } maximum flow computations. Since the average number of iterations per phase is O(m), this implies a further reduction in complexity. 4.1 Gomory-Hu Trees and an oracle for SND. A Gomory-Hu tree of an undirected, weighted graph G is an edge-weighted tree TG that concisely represents connectivity information about G. The weight of the minimum weight cut separating vertices i and j in G equals the weight of the minimum weight edge on the unique path from i to j in TG . Gomory and Hu [18] show that every graph G has a corresponding TG . Gusfield [21] describes an efficient procedure to construct TG using n - 1 maximum flow computations in G. Checking Feasibility. Let H = (V , EH ) be a graph on V that contains an edge for every pair (i, j ) with rij > 0. Note that k = |EH |. Interpret rij as the weight of edge (i, j ). In the Gomory-Hu tree TH of H , the (implied) connectivity requirement for (i, j ) is the weight of the minimum weight edge in the path from i to j in TH . Let Fe  E be the cutset of G defined by the removal of edge e in TH : The removal of e creates 1 2 subtrees TH and TH . Define Fe  E to be the set of 1 edges in G with one end point in TH and one end point 2 in TH . Let x be a vector of edges capacities for G. Let  (i, j ) be the value of the maximum flow from i to j (= value of minimum cut separating i and j ) in G with edge capacities x. Lemma 4.1. x satisfies al l connectivity requirements r if and only if  (i, j )  rij for al l (i, j )  TH .

Lemma 4.2. An estimate EST of OPT for SND that satisfies (3.5) with  = n can be found with min{n, k } minimum cost flow computations in a unit capacity graph. Proof. k  n: For edge (i, j )  TH , the minimum cost flow of value rij from i to j in the graph G with unit capacities and costs c may be assumed to be integral. The set of edges that carry flow thus define a minimum cost network that allows for connectivity rij between i and j . Let lij be the cost of this flow. Then clearly lij  OPT. Repeating this for each edge in in TH yields a collection of edges. This collection of edges is sufficient to meet all connectivity requirements by Lee ma 4.1. Thus we have that maxeTH le  OPT  m TH le  (n - 1) maxeTH le . Hence, the estimate maxeTH le satisfies (3.5) for  = n. If k < n, then we don't need to use TH and can simply compute le for each e  EH to obtain maxeEH le  OPT  k maxeEH le . Corollary 4.1. An -optimal solution to SND can be found using FractionalCoverUB with at most O( -2 log n ) phases and at most O( -2 m log n ) iterations. Proof. This follows from Theorem 3.1 and Lemma 4.2

We now explain the applicability of Lemma 4.1 to finding violated inequalities of SND. Lemma 4.3. Given TH , a set S that satisfies x(S ) < fSND (S ), if one exists, can be found with at most min{k , n} s-t minimum cut computations.

Proof. If k  n, then for each edge (s, t)  TH , compute a minimum s-t x-capacity cut Sst in G. If x(Sst )  fSND (S ) for all (s, t)  TH , then by Lemma 4.1, Proof. The "only if " direction is clear. Let u and v x(S )  fSND (S ) for all S  V . If k < n, then compute a minimum s-t x-capacity be nonadjacent vertices in TH with ru,v > 0. Let (u = w0 , w1 , . . . , wq = v ) be the adjacent vertices on cut for each (s, t)  EH . the unique path P from u to v in TH . Each cut in G Lemma 4.3 implies that min{k , n} flow computaseparating u from v must separate at least one pair of tions per iteration are sufficient to find a most vioconsecutive vertices in P . Since  (i, j )  rij for each (i, j )  TH , each cut in G separating u from v has value lated inequality for SND. Using a technique introduced at least the minimum weight edge in P , which is at least in [10], this can be reduced to a single flow computation ruv . Hence, x satisfies the connectivity requirement for per iteration, plus an additional min{k , n} iterations per phase: The results of Section 2 still hold if we replace a (u, v ). most violated inequality with any inequality that satisLemma 4.1 is useful both to obtain a good estimate fies A(p)T x < b(p). Finding some violated inequality of OPT for SND, and also for finding violated inequali- instead of the most violated one may be done more efties of SND. We first explain its applicability to obtain- ficiently. To do this, we fix an ordering of the edges ing good estimates of OPT. in TH . In the -phase, we go through this ordering

1007


exactly once. For edge (i, j )  TH , we run a maximum flow computation for this edge. If  (i, j )  rij , then we are done with the edge for the phase. Otherwise, we find an inequality  (i, j ) < rij that we can use in this iteration. We continue with the iteration in FractionalCoverUB, increasing x. In the next iteration of FractionalCoverUB, we start the search for a violated inequality with edge (i, j ). We have established the following lemma: Lemma 4.4. The number of maximum flow computations required to obtain an -optimal solution to SND via algorithm FractionalCoverUB is O( -2 m log(n/ )).

4.2 Faster implementation of FractionalCoverUB for SND. In this section, we prove Theorem 4.2, which implies the second main result of this paper: Theorem 4.1. An -optimal solution to SND can be found with min{n, k } cal ls to a minimum cost flow algorithm for uniform capacity graphs, plus at most O( -2 m(m + n min{n, k } log n) log n) additional work. Theorem 4.1 represents a run time improvement over existing algorithms by a factor of at least O(min{m, nk }), and by a factor of O(m log n) if k = (n). Theorem 4.2. The number of maximum flow computations required to obtain an -approximate solution to SND via algorithm FractionalCoverUB is at most min{n, k } per -phase. To prove Theorem 4.2, we modify the variant of push-relabel algorithm of Goldberg and Tarjan [17] that uses the gap heuristic [8]. Due to space restrictions, we sketch the main ideas here. Push-Relab el algorithm. The push relabel algorithm for maximum s-t flow starts with a preflow : a function g that satisfies capacity constraints (g (v , w)  u(v , w) for given capacity vector u) and nonnegativity constraints (g (v , w)  0), but relaxes flow conservation so that nodes other than the source s may have excess : more flow may enter the node than leave it. It also uses vertex labels d that satisfy d(s) = n, and d(v )  d(w) + 1 if g (v , w) < u(v , w). Initially, d(v ) = 0 for all v = s. The push-relabel algorithm iteratively selects a vertex v with excess(v ) > 0 and performs either a push operation or a relabel operation. Push(v , w) applies if excess(v ) > 0, g (v , w) < u(v , w), and d(v ) = d(w) + 1. It increases g (v , w) and decreases excess(v ) by min{excess(v ), u(v , w) - g (v , w)}. Relabel(v ) applies if excess(v ) > 0 and g (v , w) = u(v , w) for all w with d(v ) = d(w) + 1. It increases d(v ) by 1.

The run time of this algorithm depends on the number of relabels and push operations. Goldberg and Tarjan show that with appropriate choice of v and w in each iteration, the number of push operations is O(mn) and the number of relabel operations is O(n2 ). They show that with the use of dynamic trees, each operation 2 can be done in amortized O(log n ) time. m Gap heuristic. The gap heuristic is due independently to Cherkassky [6, 7]; and Derigs and Meier [9]. Its practical use is demonstrated in Cherkassky and Goldberg [8]. It is based on the observation that if a vertex v with label d(v ) is eligible for a relabel operation, and there are no other vertices with this label, then the set of vertices with label d(v ) or higher will be on the source side of any minimum s-t cut. The heuristic says that all vertices with label at least d(v ) can be relabeled to n without affecting the correctness or asymptotic complexity of the algorithm. We call this set of vertices the cut implied by the gap. Mo dified push-relab el with gap heuristic for finding violated cuts. At the start of the -phase, the minimum x-capacity i-j cut in G has value at least rij /(1 + ). We use the push-relabel algorithm to increase this to rij . Pro of of Theorem 4.2: In one phase, we fix an ordering of edges in TH and cycle through these edges once. For (i, j )  TH , we start using the push-relabel algorithm to compute a minimum i-j cut in G with the current iterate x as the capacity vector. If the set of arcs leaving i has capacity less than rij we use {i} as the first cut considered by algorithm FractionalCoverUB in line (2) or line (16). Otherwise, we run the push-relabel algorithm until the first gap is encountered. If the gap implies a cut S that has capacity at least rij , we simply contract all v  S into the source, and continue with the modified push-relabel. If the cut S implied by the gap has x((S )) < rij , we use S as the determining inequality in line (16), freeze the current state of the modified push-relabel algorithm, and continue with the steps in FractionalCoverUB until it returns to line (16). Inside this iteration, FractionalCoverUB increases the xcapacity of edges only in (S ) (lines (6),(11),(12)). On return to line (16), FractionalCoverUB checks if S still has x((S )) < rij , and if so, continues on with it. Otherwise, the modified push-relabel is unfrozen and resumed until a new gap is found. This repeats until this modified push-relabel algorithm finds a flow of value at least rij , indicating that the connectivity between i and j is at least rij . The total number of pushes and relabels performed by this modified algorithm is no more than the number performed by the unmodified push-relabel algorithm. A phase ends after

1008


this modified push relabel is performed once for each edge in TH . Hence, the work devoted to finding cuts in one phase is at most min{n, k } push-relabel maximum flow computations. Pro of of Theorem 4.1: The algorithmic steps can be partitioned into three groups: time spent in maximum flow computations, time spent per violated cut, and time spent to bound OPT. By Theorem 4.2, the number of push-relabel maximum flow computations over the course of the algorithm is min{n, k } times the number of phases. Thus, by Corollary 4.1 and the run time complexity of the push-relabel algorithm, the total work to execute the maximum flow computations is O( -2 mn min{n, k } log2 n). There is at most linear work per cut, and the number of cuts is bounded by the number of iterations, so the number of steps is bounded by O( -2 m2 log n). Finally, to bound OPT we use min{n, k } minimum cost flow computations in uniform capacity graphs (Lemma 4.2). Acknowledgment Thanks to Neal Young for thoughtful discussion and feedback. References
[1] D. Bienstock. Experiments with a network design algorithm using -approximate linear programs. Submitted for publication, 1996. [2] D. Bienstock. Approximately solving large-scale linear programs. i. strengthening lower bounds and accelerating convergence. Technical report, CORC Report 1999-1, Columbia University, 1999. Extended abstract: Proc. 11th Ann. ACM-SIAM Symposium on Discrete Algorithms, January 2000. [3] D. Bienstock. Potential function algorithms for approximately solving linear programs: theory and practice. CORE Lecture Series monograph. 2001. Available from Center for Operations Research and Econometrics (CORE), U. Catholique de Louvain, Belgium. [4] R. D. Carr, L. Fleischer, V. J. Leung, and C. A. Phillips. Strengthening integrality gaps for capacitated network design and covering problems. In Proceedings of the 11th Annual ACM-SIAM Symposium on Discrete Algorithms, pages 106­115, 2000. [5] J. Cheriyan, S. Venpala, and A. Vetta. Approximation algorithms for minimum-cost k-vertex connected subgraphs. In Proceedings of the 34th Annual ACM Symposium on Theory of Computing, 2002. [6] B. V. Cherkassky. A fast algorithm for computing maximum flow in a network. In A. V. Karzanov, editor, Col lected Papers: Combinatorial Methods for Flow Problems, volume 3, pages 90­96. The Institute for Systems Studies, Moscow, 1979. In Russian. [7] B. V. Cherkassky. A fast algorithm for constructing a maximum flow through a network. In Selected topics

[8]

[9]

[10]

[11]

[12]

[13]

[14]

[15]

[16]

[17]

[18] [19]

[20]

[21] [22] [23]

in discrete mathematics (Moscow, 1972­1990), number 158 in Amer. Math. Soc. Transl. Ser. 2, pages 23­30. Amer. Math. Soc., Providence, RI, 1994. B. V. Cherkassky and A. V. Goldberg. On implementing the push-relabel method for the maximum flow problem. Algorithmica, 19:390­410, 1997. U. Derigs and W. Meier. Implementing goldberg's maxflow-algorithm--a computational investigation. Z. Oper. Res., 33(6):383­403, 1989. L. K. Fleischer. Approximating fractional multicommodity flow independent of the number of commodities. SIAM J. Discrete Math., 13(4):505­520, 2000. L. K. Fleischer. A 2-approximation for minimum cost {0, 1, 2} vertex connectivity. In 8th International Integer Programming and Combinatorial Optimization Conference, pages 115­129, 2001. L. K. Fleischer, K. Jain, and D. P. Williamson. An iterative rounding 2-approximation algorithm for the element connectivity problem. In IEEE [22], pages 339­ 347. N. Garg and R. Khandekar. Fast approximation algorithms for fractional steiner forest and related problems. In 43rd Annual IEEE Symposium on Foundations of Computer Science, 2002. To appear. N. Garg and J. K¨nemann. Faster and simpler algoo rithms for multicommodity flow and other fractional packing problems. In 39th Annual IEEE Symposium on Foundations of Computer Science, pages 300­309, 1998. N. Garg, G. Konjevod, and R. Ravi. A polylogarithmic approximation algorithm for the group steiner tree problem. J. Algorithms, 37(1):66­84, 2000. A. V. Goldberg, J. D. Oldham, S. Plotkin, and C. Stein. An implementation of a combinatorial approximation algorithm for minimum-cost multicommodity flow. In R. E. Bixby, E. A. Boyd, and R. Z. R´os-Mercado, i editors, Integer Programming and Combinatorial Optimization, volume 1412 of Lecture Notes in Computer Science, pages 338­352, Berlin, 1998. Springer. A. V. Goldberg and R. E. Tarjan. A new approach to the maximum flow problem. Journal of the ACM, 35:921­940, 1988. R. E. Gomory and T. C. Hu. Multi-terminal network flows. J. Soc. Indust. Appl. Math, 9:551­570, 1961. M. D. Grigoriadis and L. G. Khachiyan. Fast approximation schemes for convex programs with many blocks and coupling constraints. SIAM Journal on Optimization, 4:86­107, 1994. M. D. Grigoriadis and L. G. Khachiyan. Coordination complexity of parallel price-directive decomposition. Mathematics of Operations Research, 21:321­340, 1996. D. Gusfield. Very simple methods for all pair network flow analysis. SIAM J. Comput., 19:143­155, 1990. 42nd Annual Symposium on Foundations of Computer Science, 2001. K. Jain. A factor 2 approximation algorithm for the generalized Steiner network problem. Combinatorica, 21(1):39­60, 2001.

1009


[24] G. Karakostas. Faster approximation schemes for fractional multicommodity flow. In Proceedings of the 13th Annual ACM-SIAM Symposium on Discrete Algorithms, 2002. [25] D. Karger and S. Plotkin. Adding multiple cost constraints to combinatorial optimization problems, with applications to multicommodity flows. In Proceedings of the 27th Annual ACM Symposium on Theory of Computing, pages 18­25, 1995. ´ [26] P. Klein, S. Plotkin, C. Stein, and E. Tardos. Faster approximation algorithms for the unit capacity concurrent flow problem with applications to routing and finding sparse cuts. SIAM Journal on Computing, 23:466­ 487, 1994. ´ [27] T. Leighton, F. Makedon, S. Plotkin, C. Stein, E. Tardos, and S. Tragoudas. Fast approximation algorithms for multicommodity flow problems. Journal of Computer and System Sciences, 50:228­243, 1995. [28] V. Melkonian and E. Tardos. Approximation algorithms for a directed network design problem. In 7th International Integer Programming and Combinatorial Optimization Conference, pages 345­360, 1999. ´ [29] S. A. Plotkin, D. Shmoys, and E. Tardos. Fast approximation algorithms for fractional packing and covering problems. Mathematics of Operations Research, 20:257­301, 1995. [30] T. Radzik. Fast deterministic approximation for the multicommodity flow problem. Mathematical Programming, 78:43­58, 1997. [31] F. Shahrokhi and D. W. Matula. The maximum concurrent flow problem. Journal of the ACM, 37:318­ 334, 1990. [32] D. P. Williamson. On the Design of Approximation Algorithms for a Class of Graph Problems. PhD thesis, MIT Computer Science, September 1993. [33] N. Young. Randomized rounding without solving the linear program. In Proceedings of the 6th Annual ACMSIAM Symposium on Discrete Algorithms, pages 170­ 178, 1995. [34] N. Young. Sequential and parallel algorithms for mixed packing and covering. In IEEE [22]. 538-546.

P (xk ) l

= P (xk-1 ) + u l = P (xk-1 ) l
l

j

A(p, j )xk-1 (j ) l
Q(p)

+ j A(p, j )]
Q(p) /

u[ (4.6)

j

A(p, j )xk-1 (j ) - (k ) l
=1

 P (xk-1 ) + u[b(p) - l

j

A(p, j )](k )
Q(p) /

= P (xk-1 ) + (D(k, l) - D(k , l - 1))(k ) l  P (xk ) + (k) 0 jl (D(k , j ) - D(k , j - 1))

=1

= P (xk ) + (k)[D(k) - D(k - 1)] 0 Inequality (4.6) follows since (k ) > A(p)T xk /b(p) l for all iterations l in the (k )-phase. Thus, at the end of the k th phase, we have that P (k )  P (k - 1) + (k)(D(k) - D(k - 1))  P (0) + i k (D(i) - D(i - 1))(i).
=1

By definition of , (i - 1) = (i)/(1 + ). Also, at the end of phase k , solution xk satisfies the property that if we scale by (k ), then xk /(k ) is feasible for (P): it satisfies both the upper bound constraints (by Lemma 2.1), and the covering constraints. Thus the optimal solution value, OPT, satisfies OPT  P (x / )  P (k )/(k ) for all 1  k  t - 1. We (i can now rewrite an upper bound for P (t), using + )  1
P (i-1) P (x / ) . t (1 + ) i (D(i) - D(i - 1))P (i - 1) P (x / ) =1

P (t)  P (0) +

App endix Pro of of Theorem 2.1: We first show  that D(y t / log1+ C (1+ ) , z t / log1+ C (1+ ) )   P (x / )/(1 + ). By feasibility of x (Lemma 2.1 and the inner while loop in FractionalCoverUB) and (y t , z t )/ log1+ C (1+ ) (Lemma 2.7), and weak LP  duality, this implies -optimality. kk Let D(k , l) = D(yl , zl ), and let D(k ) = D(y k , z k ). The change in ob jective function value of (D) in onj iteration is D(k , l) - D(k , l - 1) = u[b(p) - e We Q(p) A(p, j )], as determined in lines (13)-(15). / now bound the change in P (x) in terms of change in D. Recall that Q(p) = {1  j  m|A(p, j ) > 0 and x(j ) < }.

The analysis in [10, 14] shows that this inequality, along with Lemma 2.3, and the fact that P (0) = m , imply the following inequalities: 1  P (t)  m e
(1+ )D (t)/P (x / )

.

These in turn imply the following bound on the ratio between the best feasible solution to (GSN) found and final (unscaled) dual solution. P (x / ) (1 + )  . D(t) ln(m )-1 Choosing  = O(((C (1 + ))1- m)-1/ ) implies optimality (see, for example, the analysis in [10, 14]). For this choice of  , the number of phases is O( -2 log(mC )) by Lemma 2.3, and by Lemma 2.5 the number of oracle calls is O( -2 m log(mC )).

1010