# Approximation Algorithms for Metric Facility Location and k-Median Problems Using the Primal-Dual Schema and Lagrangian Relaxation.

Abstract. We present approximation algorithms for the metric uncapacitated facility location problem and the metric k-median problem achieving guarantees of 3 and 6 respectively. The distinguishing feature of our algorithms is their low running time: O(m log m) and O(m log m(L + log(n))) respectively, where n and m are the total number of vertices and edges in the underlying complete bipartite graph on cities and facilities. The main algorithmic ideas are a new extension of the primal-dual schema and the use of Lagrangian relaxation to derive approximation algorithms.

Categories and Subject Descriptors: F.2.2 [Analysis of Algorithms and Problem Complexity]: Nonnumerical Algorithms and Problems

General Terms: Algorithms, Theory

Additional Key Words and Phrases: Approximation algorithms, facility location problem, k-median problem, Lagrangian relaxation, linear programming

1. Introduction

Given costs for opening facilities and costs for connecting cities to facilities, the uncapacitated facility location problem seeks a minimum cost solution that connects each city to an open facility. Clearly, this problem is applicable to a number of industrial situations. For a modern-day application, consider the problem of locating proxy servers on the web. For this reason, it has occupied a central place in operations research since the early 60s,(1) and has been studied from the perspectives of worst-case analysis, probabilistic analysis, polyhedral combinatorics and empirical heuristics (see Cornuejols et al. [1990] and Nemhauser and Wolsey [1990]). In the last few years, there has been renewed interest in tackling this problem, this time from the perspective of approximation algorithms.(2) In this paper, we carry this further by developing an approximation algorithm based on the primal-dual schema. We further use this algorithm as a subroutine to solve a related problem, the k-median problem. The latter problem differs in that there are no costs for opening facilities, instead a number k is specified, which is an upper bound on the number of facilities that can be opened. The two algorithms achieve approximation guarantees of 3 and 6 respectively.

Both of our algorithms work only for the metric case, that is, when the connecting costs satisfy the triangle inequality; both problems are NP-hard for this case as well. If the connection costs are unrestricted, approximating either problem is as hard as approximating set cover, and therefore cannot be done better than O(log n) factor, unless NP [subset or equal to] P. For the first problem, this is straightforward to see, and for the second, this is established by Lin and Vitter [1992].

The distinguishing feature of our algorithms is their low running time: O(m log m) and O(m log m(L + log(n))) respectively, where n and m are the total number of vertices and edges in the underlying complete bipartite graph on cities and facilities (n = [n.sub.c] + [n.sub.f] and m = [n.sub.c] x [n.sub.f], where [n.sub.c] and [n.sub.f] are the number of cities and facilities) and L is the number of bits needed to represent a connecting cost. In particular, the running time of the first algorithm is dominated by the time taken to sort the connecting costs of edges. It is worth pointing out that our facility location algorithm is also suitable for distributed computation.

The first constant factor algorithm for the metric uncapacitated facility location problem was given by Shmoys et al. [1997] improving on Hochbaum's [1982] bound of O(log n) (see Lin and Vitter [1992] for another O(log n) factor algorithm). Their approximation guarantee was 3.16. After some improvements [Guha and Khuller 1998; Chudak 1998], the current best factor is (1 + 2/e), due to Chudak and Shmoys [1998]. The drawback of these algorithms, based on LP-rounding, is that they need to solve large linear programs, and so have prohibitive running times for most applications. A different approach was recently used by Korupolu et al. [1998]. They showed that a well-known local search heuristic achieves an approximation guarantee of (5 + [Epsilon]), for any [Epsilon] [is greater than] 0. However, even this algorithm has a high running time of O([n.sup.6] log n/[Epsilon]). Regarding hardness results, the work of Guha and Khuller [1998] and S. Sviridenko (personal communication) establishes that a better factor than 1.463 is not possible, unless NP [subset or equal to] P.

Researchers have felt that the primal-dual schema should be adaptable in interesting ways to the combinatorial structure of individual problems, and that its full potential has not yet been realized in the area of approximation algorithms. Our work substantiates this belief. We extend the scope of this schema in the following way: All primal-dual approximation algorithms obtained so far(3) work with a pair of covering and packing linear programs, that is, a primal-dual pair of LP's such that all components of the constraint matrix, objective function vector and right hand side vector are nonnegative. This includes, for instance, Williamson et al. [1995] and Goemans et al. [1994], in which the overall LP-relaxation does have negative coefficients; however, the problem is decomposed into phases, and the relaxation used in each phase is a covering program. On the other hand, our algorithm works with primal and dual programs that do have negative coefficients.

Despite this added complexity, our algorithm has a simple description: Each city j keeps raising its dual variable, [[Alpha].sub.j], until it gets connected to an open facility. All other primal and dual variables simply respond to this change, trying to maintain feasibility or satisfying complementary slackness conditions. For the latter, we give a new mechanism as well.

Until the work of Rajagopalan and Vazirani [1999] (which relaxed the dual program itself), all approximation algorithms based on the primal-dual schema used the mechanism formalized in Williamson et al. [1995]. In the first phase, an integral primal solution is found, satisfying the primal complementary slackness conditions; however, this solution may have redundancies. In the second phase, a minimal solution is extracted, typically via a reverse delete procedure, and in the process, dual complementary slackness conditions get satisfied with a relaxation factor. The final algorithm has this factor as its approximation guarantee.

Our first phase is similar. In the second phase, we first ensure that all complementary slackness conditions are satisfied; however, the primal solution may be infeasible. The solution is augmented--this time the primal conditions need to be relaxed by a factor of 3, which is also the approximation guarantee of the algorithm.

The k-median problem also has numerous applications, especially in the context of clustering, and has also been extensively studied. In recent years, the problem has found new clustering applications in the area of data mining (see Bradley et al. [1998]).

A nontrivial approximation algorithm for this problem eluded researchers for many years. The breakthrough was made by Bartal, who gave a factor O(log n log log n) algorithm using a probabilistic approximation of metric spaces by tree metrics. After a slight improvement to a factor of O(log k log log k) [Charikar et al. 1998], a constant factor algorithm was recently obtained by Charikar et al. [1999] using a technique of Lin and Vitter [1992]. Their algorithm has an approximation guarantee of 6 2/3, however, it has the same drawback since it uses LP-rounding. Their algorithm uses several ideas from the constant factor algorithms obtained for the facility location problem, thus making one wonder if there is a deeper connection between the two problems.

In this paper, we establish such a connection between the two problems: that a Lagrangian relaxation of the k-median problem is the facility location problem. This enables us to use our algorithm for the facility location problem as a subroutine to solve the k-median problem. The Lagrangian relaxation technique has been used implicitly in the past by Garg [1996] to obtain a factor 3 algorithm for the k-MST problem. In this paper, we make its use transparent. We also abstract our ideas into a general method for deriving approximation algorithms using this technique.

These ideas also help solve a common generalization of the two problems--in which facilities have costs, and in addition, there is an upper bound on the number of facilities that can be opened. We give a factor 6 approximation algorithm for this problem as well; the previous bound was 9.8 [Charikar et al. 1999].

The capacitated facility location problem, in which each facility i can serve at most [u.sub.i] cities, has no nontrivial approximation algorithms. Part of the problem is that all LP-relaxations known for this problem have unbounded integrality gap (see Shmoys et al. [1997]). In Section 5, we give a factor 4 approximation algorithm for the variant in which each facility can be opened an unbounded number of times; if facility i is opened [y.sub.i] times, it can serve at most [u.sub.i] [y.sub.i] cities. A special case of this version, in which the capacities of all the facilities are assumed to be equal, is solved with factor 3 in Chudak and Shmoys [1999], again using LP-rounding. The special case of uniform capacities is solved within a factor of 5, using at most one extra copy of a facility at each location in Chudak and Williamson [1999].

Building on ideas presented in this paper, Charikar and Guha [1999] have obtained the following improved results: a factor 1.853 algorithm for the facility location problem and a factor 4 algorithm for the k-median problem, both with running times of O([n.sup.3]).

2. The Metric Uncapacitated Facility Location Problem

The uncapacitated facility location problem seeks a minimum cost way of connecting cities to open facilities. It can be stated formally as follows: Let G be a bipartite graph with bipartition (F, C), where F is the set of facilities and C is the set of cities. Let [f.sub.i] be the cost of opening facility i, and [c.sub.ij] be the cost of connecting city j to (opened) facility i. The problem is to find a subset I [subset or equal to] F of facilities that should be opened, and a function [Phi]: C [right arrow] I assigning cities to open facilities in such a way that the total cost of opening facilities and connecting cities to open facilities is minimized. We will consider the metric version of this problem, that is, the [c.sub.ij]'s satisfy the triangle inequality.

We will adopt the following notation: [n.sub.c] = |C| and [n.sub.f] = |F|. The total number of vertices [n.sub.c] + [n.sub.f] = n and the total number of edges [n.sub.c] x [n.sub.f] = m.

Consider the following integer program for this problem, due to Balinski [1966]. In this program, [y.sub.i] is an indicator variable denoting whether facility i is open, and [x.sub.ij] is an indicator variable denoting whether city j is connected to the facility i. The first constraint ensures that each city is connected to at least one facility, and the second ensures that this facility must be open.

(1) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The LP-relaxation of this program is:

(2) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The dual program is:

(3) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

2.1. RELAXING PRIMAL COMPLEMENTARY SLACKNESS CONDITIONS. Our algorithm is based on the primal-dual schema. As stated in the introduction, instead of the usual mechanism of relaxing dual complementary slackness conditions, we relax the primal conditions. Before showing how this is done, let us give the reader some feel for how the dual variables "pay" for a primal solution by considering the following simple setting: suppose LP (2) has an optimal solution that is integral, say I [subset or equal to] F and [Phi]: C [right arrow] I. Thus, under this solution, [y.sub.i] = 1 iff i [element of] I, and [x.sub.ij] = 1 iff i = [Phi](j).

Let ([Alpha], [Beta]) denote an optimal dual solution. The reader can verify that primal and dual complementary slackness conditions imply the following facts:

--Each open facility is fully paid for, that is, if i [element of] I, then

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

--Suppose city j is connected to facility i, that is, [Phi](j) = i. Then, j does not contribute for opening any facility besides i, that is, [[Beta].sub.i'j] = 0 if i' [is not equal to] i. Furthermore, [[Alpha].sub.j] - [[Beta].sub.ij] = [c.sub.ij]. So, we can think of [[Alpha].sub.j] as the total price paid by city j; of this, [c.sub.ij] goes towards the use of edge (i, j), and [[Beta].sub.ij] is the contribution of j towards opening facility i.

Suppose the primal complementary slackness conditions were relaxed as follows, while maintaining the dual conditions:

[inverted] Aj [element of] C: (1/3)[c.sub.[Phi](j)j] [is less than or equal to] [[Alpha].sub.j] - [[Beta].sub.[Phi](j)j] [is less than or equal to] [c.sub.[Phi](j)j],

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Then, the cost of the (integral) solution found would be within thrice the dual found, thus leading to a factor 3 approximation algorithm. However, we would like to obtain the stronger inequality stated in Theorem 7, in which the dual pays at least one-third the connection cost, but must pay completely for opening facilities. This stronger inequality will be needed in order to use this algorithm to solve the k-median problem.

For this reason, we will relax the primal conditions as follows: The cities are partitioned into two sets, directly connected and indirectly connected. Only directly connected cities will pay for opening facilities, that is, [[Beta].sub.ij] can be non-zero only if j is a directly connected city and i = [Phi](j). For an indirectly connected city j, the primal condition is relaxed as follows:

(1/3)[c.sub.[Phi](j)j] [is less than or equal to] [[Alpha].sub.j] [is less than or equal to] [c.sub.[Phi](j)j].

All other primal conditions are maintained, that is for a directly connected city j,

[[Alpha].sub.j] - [[Beta].sub.[Phi](j)j] = [c.sub.[Phi](j)j],

and for each open facility i,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

2.2. THE ALGORITHM. Our algorithm consists of two phases. In Phase 1, the algorithm operates in a primal-dual fashion. It finds a dual feasible solution, and also determines a set of tight edges and temporarily open facilities, [F.sub.t]. Phase 2 consists of choosing a subset I of [F.sub.t] to open, and finding a mapping, [Phi], from cities to I.

Algorithm 1

Phase 1. We would like to find as large a dual solution as possible. This motivates the following underlying process for dealing with the non-covering-packing pair of LP's. Each city j keeps raising its dual variable, [[Alpha].sub.j], until it gets connected to an open facility. All other primal and dual variables simply respond to this change, trying to maintain feasibility or satisfying complementary slackness conditions.

A notion of time is defined in this phase, so that each event can be associated with the time at which it happened; the phase starts at time 0. Initially, each city is defined to be unconnected. Throughout this phase, the algorithm raises the dual variable [[Alpha].sub.j] for each unconnected city j uniformly at unit rate, that is, [[Alpha].sub.j] will grow by 1 in unit time. When [[Alpha].sub.j] = [c.sub.ij] for some edge (i, j), the algorithm will declare this edge to be tight. Henceforth, dual variable [[Beta].sub.ij] will be raised uniformly, thus ensuring that the first constraint in LP (3) is not violated, [[Beta].sub.ij] goes towards paying for facility i. Each edge (i, j) such that [[Beta].sub.ij] [is greater than] 0 is declared special.

Facility i is said to be paid for if [[Sigma].sub.j][[Beta].sub.ij] = [f.sub.i]. If so, the algorithm declares this facility temporarily open. Furthermore, all unconnected cities having tight edges to this facility are declared connected and facility i is declared the connecting witness for each of these cities. (Notice that the dual variables [[Alpha].sub.j] of these cities are not raised anymore.) In the future, as soon as an unconnected city j gets a tight edge to i, j will also be declared connected and i will be declared the connecting witness for j (notice that [[Beta].sub.ij] = 0, and so edge (i, j) is not special). When all cities are connected, the first phase terminates. If several events happen simultaneously, the algorithm executes them in arbitrary order.

Remark 2. At the end of Phase 1, a city may have paid towards temporarily opening several facilities. However, we want to ensure that a city pay for only the facility that it is eventually connected to. This is ensured in Phase 2, which chooses a subset of temporarily open facilities for opening permanently.

Phase 2. Let [F.sub.t] denote the set of temporarily open facilities and T denote the subgraph of G consisting of all special edges. Let [T.sup.2] denote the graph that has edge (u, v) iff there is a path of length at most 2 between u and v in T, and let H be the subgraph of [T.sup.2] induced on [F.sub.t]. Find any maximal independent set in H, say I. All facilities in the set I are declared open.

For city j, define [F.sub.j] = {i [element of] [F.sub.t] | (i, j) is special}. Since I is an independent set, at most one of the facilities in [F.sub.t] is opened. If there is a facility i [element of] [F.sub.j] that is opened, then set [Phi](j) = i, and declare city j directly connected. Otherwise, consider tight edge (i', j) such that i' was the connecting witness for j. If i' [element of] I, again set [Phi](j) = i' and declare city j directly connected (notice that in this case [[Beta].sub.i'j] = 0). In the remaining case that [is not an element of] I, let i be any neighbor of i' in graph H such that i [element of] I. Set [Phi](j) = i and declare city j indirectly connected.

I and [Phi] define a primal integral solution: [x.sub.ij] = 1 iff [Phi](j) = i, and [y.sub.i] = 1 iff i [element of] I. The values of [[Alpha].sub.j] and [[Beta].sub.ij] obtained at the end of Phase 1 form a dual feasible solution.

2.3. ANALYSIS. We will show how the dual variables [[Alpha].sub.j]'s pay for the primal costs of opening facilities and connecting cities to facilities. Denote by [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] the contributions of city j to these two costs respectively; [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. If j is indirectly connected, then [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. If j is directly connected, then the following must hold:

[[Alpha].sub.j] = [c.sub.ij] + [[Beta].sub.ij],

where i = [Phi](j). Now, let [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

LEMMA 3. Let i [element of] I. Then,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

PROOF. Since i is temporarily open at the end of Phase 1, it is completely paid for, that is,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

The critical observation is that each city j that has contributed to [f.sub.i] must be directly connected to i. For each such city, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Any other city j' that is connected to facility i must satisfy [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. The lemma follows. []

COROLLARY 4. [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Recall that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] was defined to be 0 for indirectly connected cities. So, only the directly connected cities pay for the cost of opening facilities.

LEMMA 5. For an indirectly connected city [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], where i = [Phi](j).

PROOF. Let i' be the connecting witness for city j. Since j is indirectly connected to i, (i, i') must be an edge in H. In turn, there must be a city, say j', such that (i, j') and (i', j') are both special edges. Let [t.sub.1] and [t.sub.2] be the times at which i and i' were declared temporarily open during Phase 1.

Since edge (i', j) is tight, [[Alpha].sub.j'] [is greater than or equal to] [c.sub.i'j']. We will show that [[Alpha].sub.j] [is greater than or equal to] [c.sub.ij'] and [[Alpha].sub.j] [is greater than or equal to] [c. sub.i'j']. Then, the lemma will follow by using the triangle inequality.

Since edges (i', j') and (i, j') are tight, [[Alpha].sub.j'] [is greater than or equal to] [c.cub.ij'] , and [[Alpha].sub.j'] [is greater than or equal to] [c.sub.i'j']. During Phase 1, [[Alpha].sub.j'] stops growing as soon as one of the facilities that j' has a tight edge to opens. Therefore, [[Alpha].sub.j'] [is less than or equal to] min ([t.sub.1], [t.sub.2]). Finally, since i' is the connecting witness for j, [[Alpha].sub.j] [is greater than or equal to] [t.sub.2]. Therefore, [[Alpha].sub.j] [is greater than or equal to] [[Alpha].sub.j'], and the required inequalities follow. []

[ILLUSTRATION OMITTED]

REMARK 6. If instead of picking all special edges in T, all tight edges were picked, then Lemma 5 does not hold. However, if the facilities in H are ordered in the order in which they were temporarily opened, and I is picked to be the lexicographically first maximal independent set, then Lemma 5 holds again.

THEOREM 7. The primal and dual solutions constructed by the algorithm satisfy:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

PROOF. For a directly connected city [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], where [Phi](j) = i. Combining with Lemma 5 we get

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Adding this to the following inequality obtained from Corollary 4 gives the theorem:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. []

2.4. RUNNING TIME. Sort all the edges by increasing cost--this gives the order and the times at which edges go tight. For each facility, i, we maintain the number of cities that are currently contributing towards it, and the anticipated time, [t.sub.i], at which it would be completely paid for if no other event happens on the way. Initially all [t.sub.i]'s are infinite, and each facility has 0 cities contributing to it. The [t.sub.i]'s are maintained in a binary heap so we can update each one and find the current minimum in O (log [n.sub.f]) time. Two types of events happen, and they lead to the following updates:

--An edge (i, j) goes tight.

--If facility i is not temporarily open, then it gets one more city contributing towards its cost. The amount contributed towards its cost at the current time can be easily computed. Therefore, the anticipated time for facility i to go be paid for can be recomputed in constant time. The heap can be updated in O (log [n.sub.f]) time.

--If facility i is already temporarily open, city j is declared connected, and [[Alpa].sub.j] is not raised anymore. For each facility i' that was counting j as a contributor, we need to decrease the number of contributors by 1, and recompute the anticipated time at which it gets paid for.

--Facility i is completely paid for. In this event, i will be declared temporarily open, and all cities contributing to i will be declared connected. For each of these city, we will execute the second case of the previous event, that is, update facilities that they were contributing towards.

The next theorem follows by observing that each edge (i, j) will be considered at most twice. First, when it goes tight. Second, when city j is declared connected. For each consideration of this edge, we will do O (log [n.sub.f]) work.

THEOREM 8. Algorithm 1 achieves an approximation factor of 3 for the facility location problem, and has a running time of O (m log m).

2.5. TIGHT EXAMPLE. The following infinite family of examples shows that the analysis of our algorithm is tight: The graph has n cities, 1, 2, ..., n and two facilities 1 and 2. For each city j, [c.sub.2j] = 1. [c.sub.11] = 1 and all other [c.sub.ij]'s follows from the tight triangle inequalities, [f.sub.1] and [f.sub.2] are [element of] and (n + 1) [Epsilon] respectively, for a small number [Epsilon].

The optimal solution is to open facility 2 and connect all cities to it, at a total cost of (n + 1) [Epsilon] + n. Algorithm 1 will however open facility 1 and connect all cities to it, at a total cost of [Epsilon] + 1 + 3(n - 1).

2.6. EXTENSION TO ARBITRARY DEMANDS. A small extension to Algorithm 1 enables it to handle the following generalization to arbitrary demands. For each city j, a nonnegative demand [d.sub.j] is specified; any open facility can serve this demand. The cost of serving this demand via facility i is [c.sub.ij][d.sub.i].

The only change to IP (1) and LP (2) is that in the objective function, [c.sub.ij][x.sub.ij] is replaced by [c.sub.ij][d.sub.j][x.sub.ij]. This changes the first constraint in the dual (3) to

[inverted] Ai [element of] F, j [element of] C: [[Alpha].sub.j] - [[Beta].sub.ij] [is less than or equal to] [c.sub.ij][d.sub.j].

The only change to Algorithm 1 is that for each city j, [[Alpha].sub.j] is raised at rate [d.sub.j]. Notice that because of the change in the first constraint in the dual, edge (i, j) still goes tight at time [c.sub.ij]. However, once (i, j) goes tight, [[Beta].sub.ij] will be increasing at rate [d.sub.j], and so facility i may get opened earlier than in the unit demands case.

An easy way to see that this modification works is to reduce to the unit demands case by making [d.sub.j] copies of city j. The change proposed above to Algorithm 1 is more general, since it works even if [d.sub.j] is nonintegral, and even if it is exponentially large.

3. The Metric k-Median Problem

The k-median problem differs from the facility location problem in two respects: there is no cost for opening facilities, and there is an upper bound, k, on the number of facilities that can be opened; k is not fixed, it is supplied as part of the input. Once again, we will assume that the edge costs satisfy the triangle inequality.

The power of primal-dual algorithms lies in efficiently making "judicious" local improvements. On the other hand, the constraint that at most k facilities be opened is a global constraint--one that is not easy to satisfy through such an algorithm. We observe that the Lagrangian relaxation of the k-median problem is the facility location problem. This enables us to replace this global constraint by a penalty for opening each facility.

Following is an integer program for the k-median problem. The indicator variables [y.sub.i] and [x.sub.ij] play the same role as in (1).

(4) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The LP-relaxation of this program is:

(5)[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The dual program is:

(6)[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

3.1. THE HIGH LEVEL IDEA. The similarity in the linear programs of the two problems is exploited as follows: Take an instance of the k-median problem, assign a cost of z for opening each facility, and find optimal solutions to LP (2) and LP (3), say (x,y) and ([Alpha], [Beta]), respectively. By the strong duality theorem,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Now, suppose that the primal solution (x,y) happens to open exactly k facilities (fractionally), that is, [[Sigma].sub.i] [y.sub.i] = k. Then, we claim that (x,y) and ([Alpha], [Beta], z) are optimal solutions to LP (5) and LP (6) respectively. Feasibility is easy to check. Optimality follows by substituting [[Sigma].sub.i][y.sub.i] = k in the above equality, and rearranging terms to show that the primal and dual solutions achieve the same objective function value:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Let's use this idea, together with Algorithm 1 and Theorem 7, to obtain a "good" integral solution to LP (5). Suppose with a cost of z for opening each facility, Algorithm 1 happens to find solutions (x, y) and ([Alpha], [Beta]), where the primal solution opens exactly k facilities. By Theorem 7,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Now, observe that (x, y) and ([Alpha], [Beta], z) are primal (integral) and dual feasible solutions to the k-median problem satisfying

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Therefore, (x, y) is a solution to the k-median problem within thrice the optimal.

Notice that proof of factor 3 given above would not work if less than k facilities were opened; if more than k facilities are opened, the solution is infeasible for the k-median problem. The remaining problem is to find a value of z so that exactly k facilities are opened. Several ideas are required for this. The first is the following principle from economics: taxation is an effective way of controlling the amount of goods coming across the border--raising tariffs will reduce in-flow and vice versa. In a similar manner, raising z should reduce the number of facilities opened and vice versa.

It is natural now to seek a modification to Algorithm 1 that can find a value of z so that exactly k facilities get opened. This would lead to a factor 3 approximation algorithm. We don't know if this is possible. Instead, we present the following strategy which leads to a factor 6 algorithm. For the rest of the discussion, assume that we never encountered a run of the algorithm which resulted in exactly k facilities being opened.

Clearly, when z = 0, the algorithm will open all facilities, and when z is very large, it will open only one facility. The later value of z can be picked to be [nc.sub.max], where [c.sub.max] is the length of the longest edge. We will conduct a binary search on the interval [0, [nc.sub.max]] to find [z.sub.2] and [z.sub.1] for which the algorithm opens [k.sub.2] [is less than] k and [k.sub.1] [is less than] k facilities respectively, and furthermore, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], where [c.sub.min] is the length of the shortest non-zero edge. Let ([x.sup.s], [y.sup.s]) and ([x.sup.l], [y.sup.l]) be the two primal solutions found, with [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (the superscripts s and l denote "small" and "large" respectively). Further, let ([[Alpha].sup.s], [[Beta.sup.s]) and ([[Alpha].sup.l], [[Beta].sup.l]) be the corresponding dual solutions found.

Let (x, y) = a ([x.sup.s], [y.sup.s]) + b ([x.sup.l], [y.sup.l]) be a convex combination of these two solutions, with [ak.sub.1] + [bk.sub.2] = k; under these conditions, a = ([k.sub.2] - k) / ([k.sub.2] - [k.sub.1]) and b = (k - [k.sub.1])/([k.sub.2] - [k.sub.1]). Since (x, y) is a feasible (fractional) solution to the facility location problem that opens exactly k facilities, it is also a feasible (fractional) solution to the k-median problem. In this solution, each city is connected to at most two facilities.

LEMMA 9. The cost of (x, y) is within a factor of (3 + 1/[n.sub.c]) of the cost of an optimal fractional solution to the k-median problem.

PROOF. By Theorem 7, we have:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Since [z.sub.1] [is greater than] [z.sub.2], ([[Alpha].sup.l], [[Beta].sup.l]) is a feasible dual solution to the facility location problem even if the cost of facilities is [z.sup.1]. We would like to replace [z.sub.2] by [z.sub.1] in the second inequality, at the expense of the increased factor. This is achieved using the upper bound on [z.sub.1] - [z.sub.2], and the fact that [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] We get:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Multiplying this inequality by b and the first inequality by a and adding, we get

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

where [Alpha] = a[[Alpha].sup.s] + b[[Alpha].sup.l] Let [Beta] = a[[Beta].sup.s] + b[[Beta].sup.l]. Observe that ([Alpha], [Beta], [z.sub.1]) is a feasible solution to the dual of the k-median problem. The lemma follows. []

In Section 3.2, we give a randomized rounding procedure that obtains an integral solution to the k-median problem from (x, y) with an increase in cost by a small factor. In Section 3.3, we derandomize this procedure.

3.2. RANDOMIZED ROUNDING. We give a randomized rounding procedure that produces an integral solution to the k-median problem from (x, y). In the process, it increases the cost by a multiplicative factor of 1 + max(a, b).

Let A and B be the sets of facilities opened in the two solutions, |A| = [k.sub.1] and |B| = [k.sub.2]. For each facility in A, find the closest facility in B--these facilities are not required to be distinct. Let B' [subset] B be these facilities. If |B'| [is less than] [k.sub.1], arbitrarily include additional facilities from B - B' into B' until |B'| = [k.sub.1].

With probability a, open all facilities in A, and with probability b = 1 - a, open all facilities in B' In addition, a set of cardinality k - [k.sub.1] is picked randomly from B - B' and facilities in this set are opened. Notice that each facility in B - B' has a probability of b of being opened. Let I be the set of facilities opened, |I| = k.

The function [Phi]: C [right arrow] I is defined as follows: Consider city j, and suppose that it is connected to [i.sub.1] [element of] A and [i.sub.2] [element of] B in the two solutions. If [i.sub.2] [is element of] B', then one of [i.sub.1] and [i.sub.2] is opened by the procedure given above, [i.sub.1] with probability a and [i.sub.2] with probability b. City j is connected to the open facility.

If [i.sub.2] [element of] B - B', let [i.sub.3] [element of] B' be the facility in B that is closest to [i.sub.1]. City j is connected to [i.sub.2], if it is open. Else, it is connected to [i.sub.1], if it is open. If neither [i.sub.2] or [i.sub.1] is open, then j is connected to [i.sub.3].

[ILLUSTRATION OMITTED]

Denote by cost(j) the connection cost for city j in the fractional solution (x, y); [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

LEMMA 10. The expected connection cost for city j in the integral solution, E[[c.sub.[Phi](j)j] [is less than or equal to] (1 + max(a, b))cost(j). Moreover, E[[c.sub.[Phi](j)j]] can be efficiently computed.

PROOF. If [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Consider the second case, that [i.sub.2] [is not an element of] B'. Now, [i.sub.2] is open with probability b. The probability that [i.sub.2] is not open and [i.sub.1] is open is (1 - b)a = [a.sup.2], and the probability that both [i.sub.2] and [i.sub.1] are not open is (1 - b)(1 - a) = ab. This gives

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Since [i.sub.3] is the facility in B that is closest to [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII], where the second inequality follows from the triangle inequality. Again, by the triangle inequality, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Therefore,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Now, [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. Therefore,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Clearly, in both cases, E[[c.sub.[Phi](j)j]] is easy to compute.[]

Let ([x.sup.k], [y.sup.k]) denote the integral solution obtained to the k-median problem by this randomized rounding procedure. Then,

LEMMA 11

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

and moreover, the expected cost of the solution found can be computed efficiently.

3.3. DERANDOMIZATION. Derandomization follows in a straightforward manner using the method of conditional expectation. First, the algorithm opens the set A with probability a, and the set B' with probability b = 1 - a. Pick A, and compute the expected value if k - [k,sub.1] facilities are randomly chosen from B - B'. Next, do the same by picking B' instead of A. Choose to open the set that gives the smaller expectation.

Second, the algorithm opens a random subset of k - [k.sub.1] facilities from B - B'. For a choice D [subset] B - B', |D| [is less than or equal to] k - [k.sub.1], denote by E[D, B - (B' [union] D)] the expected cost of the solution if all facilities in D and additionally k - [k.sub.1] - |D| facilities are randomly opened from B - (B' [union] D). Since each facility of B - (B' [union] D) is equally likely to be opened, we get

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

This implies that there is an i such that

E[D [union] {i}, B - (B' [union] D [union] {i})] [is less than or equal to] E[B', B - (B' [union] D)].

Choose such an i and replace D by D [union] {i}. Notice that the computation of E[D [union] {i}, B - (B' [union] D [union] {i})] can be done as in Lemma 11.

3.4. RUNNING TIME. It is easy to see that a [is less than or equal to] 1 - 1/[n.sub.c] (this happens for [k.sub.1] = k - 1 and [k.sub.2] = [n.sub.c]) and b [is less than or equal to] 1 - 1/k (this happens for [k.sub.1] = 1 and [k.sub.2] = k + 1). Therefore, 1 + max(a, b) [is less than or equal to] 2 - 1/[n.sub.c]. Altogether, the approximation guarantee is (2 - 1/[n.sub.c])(3 + 1/[n.sub.c]) [is less than] 6. Using the method of conditional probabilities, this procedure can be derandomized, as in Section 3.3. The binary search will make O([log.sub.2]([n.sup.3][c.sub.max]/[c.sub.min])) = O(L + log n) probes. The running time for each probe is dominated by the time taken to run Algorithm 1; randomized rounding takes O(n) time and derandomization takes O(m) time. Hence, we get

THEOREM 12. The algorithm given above achieves an approximation factor of 6 for the k-median problem, and has a running time of O(m log m(L + log(n))).

The running time of the algorithm can also be made strongly polynomial by standard method of discretizing the costs to integers of magnitude O(poly).

3.5. TIGHT EXAMPLE. We do not have a tight example of factor 6 for the complete k-median algorithm. However, we give below an infinite family of instances which show that the analysis of the randomized rounding procedure cannot be improved.

The two solutions ([x.sup.s], [y.sup.s]) and ([x.sup.l], [y.sup.l]) open one facility, [f.sub.0], and k + 1 facilities, [f.sub.1], ..., [f.sub.k+1] respectively. The distance between [f.sub.0] and any other [f.sub.i] is 1, and that between two facilities in the second set is 2. All n cities are at a distance of 1 from [f.sub.0], and at a distance of [Epsilon] from [f.sub.k+1]. The rest of the distances are given by the triangle inequality. The convex combination is constructed with a = 1/k and b = 1 - 1/k.

Now, the cost of the convex combination is an + b[Epsilon]n. Suppose the algorithm picks [f.sub.1] as the closest neighbor of [f.sub.0]. Now, the expected cost of the solutions produced by the randomized rounding procedure is n(b[Epsilon] + [a.sup.2] + ab(2 + [Epsilon])). Letting [Epsilon] tend to 0, the cost of the convex combination is essentially na, and that of the rounded solution is na(1 + b).

3.6. A LAGRANGIAN RELAXATION TECHNIQUE FOR APPROXIMATION ALGORITHMS. Lagrangian relaxation is a fundamental technique in combinatorial optimization. In this section, we will abstract the ideas developed above to give one method of using this technique to derive approximation algorithms. This method does not require the constraints of the problem to be linear, and in fact we will present it in a very general setting.

Let [P.sub.1] be the following optimization problem:

(7) minimize f(x) subject to P(x)(7) g(x) = k

where f and g are arbitrary real valued functions, P is an arbitrary predicate, and k is a constant. Let [OPT.sub.1] denote the optimal value of this problem, and let a be the value of x at which the optimum is attained.

The Lagrangian relaxation technique consists of relaxing certain constraints by moving them into the objective function, together with associated Lagrange multipliers. We will use this technique to relax the constraint g(x) = k. Let z be the Lagrange multiplier. Now, for any value of z,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

is a lower bound on [OPT.sub.1]. To see this notice that substituting x = a in the above expression gives [OPT.sub.1]. Therefore,

(8) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

is also a lower bound on [OPT.sub.1]. Let L be the value of this expression. Let us rewrite this expression as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Now, for each value of z, consider the following optimization problem, which we will call [P.sub.2](z):

(9) minimize f(x) + zg(x)

(10) subject to P(x)

Let [OPT.sub.2](z) denote the optimum value of this problem. We will show how to derive an approximation algorithm for problem [P.sub.1] using an approximation algorithm for problem [P.sub.2].

Let A be an approximation algorithm which, for each z, finds a solution x satisfying

f(x) + [Alpha]zg(x) [is less than or equal to] [Alpha][OPT.sub.2](z)

for some constant [Alpha] [is greater than or equal to] 1. Notice that we have multiplied one term on the left-hand side by [Alpha] as well, and so this is stronger than an a factor approximation algorithm for problem [P.sub.2]. It must pick a solution so zg(x) is completely paid for by [OPT.sub.2](z).

THEOREM 13. Suppose there exists approximation algorithm A defined above. Suppose further that there is a polynomial time procedure R that uses A as a subroutine and finds a value of z for which the solution found by A satisfies g(x) = k. Then, there is an a factor approximation algorithm for problem [P.sub.1].

PROOF. By the premise, we can find in polynomial time a value of z and a solution x such that

f(x) + [Alpha]zg(x) [is less than or equal to] [Alpha][OPT.sub.2](z) and g(x) = k.

Substituting, we get

f(x) [is less than or equal to] [Alpha]([OPT.sub.2](z) - zk).

The important observation is that for any value of z,

[OPT.sub.2](z) - zk [is less than or equal to] L,

since L was defined to be the optimal value of expression (8). Therefore, f(x) [is less than or equal to] [Alpha]L. Since L is a lower bound on [OPT.sub.1], we get f(x) [is less than or equal to] [Alpha][OPT.sub.1]. Since g(x) = k, x is a feasible solution to problem [P.sub.1]; moreover, it comes within an [Alpha] factor of the optimal. []

Procedure R will be problem dependent. For instance, for the k-MST problem, after getting the two solutions, Garg (personal communication) uses additional structural properties to obtain a tree containing exactly k vertices.

For the k-median problem presented above, this involved doing a binary search to find two very close values of z for which g attains values [k.sub.1] and [k.sub.2] with [k.sub.1] [is less than] k [is less than] [k.sub.2], taking a convex combination of these solutions and doing a randomized rounding to get an integral solution with a further slight loss in the approximation factor. Other than the last step of randomized rounding, the remaining steps apply to any problem with linear constraints.

For instance, consider the following variant of the k-median problem. Instead of being specified a bound k on the number of facilities to be opened, we are specified the cost of opening each facility and an upper bound allowed for opening facilities. Subject to this constraint, the problem is to minimize the total connection cost. For this problem, we do not know how to carry out the last step of randomized rounding, and leave this as an open problem.

Another interesting phenomenon, which we call decoupling, can happen when we take Lagrangian relaxation. Suppose we have two kind of facilities, hospital and school. Suppose the total number of hospitals and schools we can open is at most k (in practice, this might be the result of a budget constraint) so that each city is connected to one hospital and one school. This problem can be thought of as two facility location problems coupled with a k-median kind of constraint. If we take its Lagrangian relaxation, we get rid of k-median kind of constraint and get two independent instances of the facility location problem, which can be solved separately.

4. A Common Generalization of the Two Problems

Consider the uncapacitated facility location problem with the additional constraint that at most k facilities can be opened. This is a common generalization of the two problems solved in this paper--if k is made [n.sub.f], we get the first problem and if the facility costs are set to zero, we get the second problem.

The techniques of this paper yield a factor 6 algorithm for this generalization as well. The high-level idea is as follows: Using the Lagrangian relaxation technique, we will first remove the restriction that at most k facilities be opened, and instead set the cost of opening each facility i to [f.sub.i] + z. Now, binary search on z will yield two values of z, close to each other, for which Algorithm 1 opens [k.sub.1] [is less than] k and [k.sub.2] [is greater than] k facilities respectively. An appropriate convex combination of these two solutions gives a fractional solution that opens exactly k facilities, with the additional property that each city is connected to at most two facilities. The cost of this solution is within thrice the cost of an optimal fractional solution. Notice that the randomized rounding procedure it ensures that the expected cost of opening facilities in the rounded solution is the same as the cost of opening facilities in the convex combination. Finally, the derandomization procedure can also be carried out in this setting.

THEOREM 14. There is a factor 6 approximation algorithm for common generalization of uncapacitated facility location and k-median problems in which facilities have costs and at most k of them can be opened.

5. Dealing with Capacities

We consider the following variant of the capacitated metric facility location problem. Each facility can be opened an unbounded number of times; if facility i is opened [y.sub.i] times, it can serve at most [u.sub.i][y.sub.i] cities. The LP-relaxation of this problem has the following extra constraint:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Let the dual variable corresponding to this constraint be [[Gamma].sub.i]. Then, the dual program is:

(11) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

For each facility i, let us fix [[Gamma].sub.i] = 3[f.sub.i]/4[u.sub.i]. This step enables us to get rid of the variables [[Gamma].sub.i] from LP (11), and the resulting linear program is again the dual of an uncapacitated facility location problem. The primal program for this modified dual is:

(12) [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

It is easy to see that [c.sub.ij] + (3[f.sub.i]/4[u.sub.i]) still satisfies the triangle inequality. Using Algorithm 1, we can now find a 0/1 integral solution to this LP satisfying

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

by Theorem 7. Now, our solution to the capacitated problem is: [x.sub.ij]'s are as in this solution, and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

This gives the following relationship between [y.sub.i] and [Y.sub.i]:

[y.sub.i] [is less than or equal to] [Y.sub.i] + [[Sigma].sub.j [element of] c][x.sub.ij]/[u.sub.i].

Using this relationship and the above inequality, we get:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

This implies

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

thereby giving an approximation guarantee of factor 4.

Remark 15. Generalizations of the problems considered in Sections 4 and 5 to the case of arbitrary demands for cities can also be solved within the factors given above, using ideas from Section 2.6.

6. [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] Clustering

Our k-median algorithm extends, in a fairly straightforward manner, to obtaining a constant factor algorithm for the problem of [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] clustering. This holds even for the case that the number of clusters and the dimension of the space are arbitrary--a case for which constant factor algorithms were not observed before. However, we note that such a result follows from previous constant factor k-median algorithms as well. The result below should be considered preliminary --the factor obtained is too high. See Drineas et al. [1999] for a factor 2 algorithm for the case that k is fixed.

Given a set of n points S = {[v.sub.1], ..., [v.sub.n]} in d-dimensional space and a number k, the problem is to find a minimum cost k-clustering, that is, to find k points, called centers, [f.sub.1], ..., [f.sub.k], so as to minimize the sum of squares of distances from each point [v.sub.i] to its closest center. This naturally defines a partitioning of the n points into k clusters.

Suppose points [v.sub.1], ... [v.sub.t] form one of these clusters with center [f.sub.1]. Define the centroid of [v.sub.1], ... [v.sub.t] to be c = ([v.sub.1] + ... + [v.sub.t])/t. It is well known that

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII],

where [parallel]u - v[parallel] denotes the square of the Euclidean distance between points u and v. So, each center must be the centroid of its cluster. Therefore, this problem can be stated as a k-median problem. The cities are the n given points, [v.sub.1], ..., [v.sub.n], and the facilities are the centroid of each subset of points. The cost of connecting a city to a facility is the square of the Euclidean distance between them. Since there are exponentially many facilities, the corresponding LP is exponential sized, and we do not know how to deal efficiently with it.

One way of getting around this difficulty is to choose centers from the given points only. Suppose the closest point from c to [v.sub.1], ..., [v.sub.t] is [v.sub.1], say. Then, using the above equality, we get

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII].

Therefore, the cost of the optimal clustering with the given points as centers is within a factor of 2 of the optimal clustering on centroids. The former problem can be expressed as a polynomial sized k-median LP, and its Lagrangian relaxation as a polynomial sized facility location LP. Our facility location algorithm solves the Lagrangian relaxation with a factor of 9. The reason for the larger factor is that edge costs do not satisfy the triangle inequality. Instead, the statement of Lemma 5 needs to be modified to [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]. For the same reason, the factor for randomized rounding also increases to 6. This gives an overall factor of 2 x 9 x 6 = 108 for [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]-clustering.

7. Discussion

A large fraction of the theory of approximation algorithms, as we know it today, is built around linear programming, which provides two main algorithm design techniques: rounding and the primal-dual schema. Both techniques have yielded algorithms with good approximation guarantees, often achieving the integrality gap of the relaxation being used. However, with respect to the running times of the algorithms derived, the two methods differ widely. Rounding resorts to the "big hammer" approach of solving the linear program and therefore leads to inefficient algorithms. On the other hand, the primal-dual schema leaves enough room to exploit the special combinatorial structure of individual problems and has therefore lead to efficient algorithms. Once the algorithm is obtained, typically the scaffolding of linear programming can be completely dispensed with to obtain a purely combinatorial algorithm. As was done in this paper, it seems worthwhile examining various algorithms derived using rounding, to see if efficient combinatorial algorithms achieving the same factors can be obtained.

Besides the objective measure of running time, another aspect in which primal-dual algorithms are superior to rounding based algorithms is the ease with which the core algorithmic idea can be modified, generalized and adapted to special circumstances or variants of the original problem. In this respect, our algorithm has met special success. Besides the various generalizations covered in this paper, the core idea has been used to solve:

--A fault tolerant version of the facility location problem, in which we are given a connectivity requirement ri with each city j, specifying the number of open facilities city j should be connected to Jain and Vazirani [2000].

--The prize collecting version of both facility location and k-median problems. In this version, we are not required to connect each city to an open facility; however, there is a specified penalty which we have to pay if a city is not connected [Charikar et al. 2001].

--The outlier version of the facility location problem, in which we are specified a number T, and are required to connect only T cities to open facilities [Charikar et al. 2001]. Charikar et al. [2001] reduce this problem to the prize collecting version by using the Lagrangian relaxation technique.

--The on-line median problem, in which k is not prespecified and is chosen on-line [Mettu and Plaxton 2000].

It is instructive to compare the current status of primal-dual approximation algorithms with the (mature) status of exact primal-dual algorithms. In the latter setting, only one underlying mechanism is used: iteratively ensuring all complementary slackness conditions. On termination, an optimal (integral) solution to the LP is obtained. In the former setting, we are not seeking an optimal solution to the LP (since the LP may not have any optimal integral solutions), and so there is a need to introduce a further relaxation. Relaxing complementary slackness conditions (which itself can be carried out in more than one way) is only one of the possibilities (see Rajagopalan and Vazirani [1999] for an alternative mechanism). Another point of difference is that in the exact setting, more sophisticated dual growth algorithms have been given, for example, Edmonds [1965]. In the approximation setting, other than Rajagopalan and Vazirani [1999], all primal-dual algorithms use a simple greedy dual growth algorithm.

So far, the primal-dual schema has been used for obtaining good integral solutions to an LP-relaxation. However, it seems powerful enough for the following more general scenario: when the NP-hard problem is captured not through an integer program, but in some other manner, and there is an LP that provides a relaxation of the problem. In this setting, the primal-dual schema will try to find solutions that are feasible for the original NP-hard problem, and are near-optimal in quality. This open problem was first mentioned in Vazirani [1995].

In Section 3.6, we have stated our Lagrangian relaxation technique in a very general setting in which the constraints of the problem are provided by arbitrary predicates. This includes, for instance, the possibility of nonlinear constraints. It will be interesting to see if this technique finds applications in non-linear settings. It will also be interesting to derive an approximation algorithm for a problem in which there are two global constraints, via the Lagrangian relaxation technique, for instance, the outlier k-median problem, in which we are specified the number of facilities that can be opened and the number of cities that need to be connected.

At a more detailed level, the issue of modifying Algorithm 1 so it opens exactly k facilities deserves some thought--this is a possible avenue for improving the factor for the k-median problem. It would be nice to improve the running time of the facility location algorithm in case the metric is specified as the closure of a sparse graph, rather than a complete bipartite graph. Another question is to obtain a non-trivial approximation algorithm for the capacitated facility location problem.

ACKNOWLEDGMENT. We wish to thank Pete Veinott for interesting discussions. Thanks to David Williamson sharing his insights about the Lagrangian relaxation technique and for pointing out that Garg's k-MST algorithm can be viewed as the use of the Lagrangian relaxation technique. Thanks also to Moses Charikar for suggesting a simplification to our facility location algorithm.

(1) See, for example, Balinski [1996], Kuehn and Hamburger [1963], and Stollsteimer [1961; 1963].

(2) See, for example, Chudak and Shmoys [1998], Guha and Khuller [1998], Korupolu et al. [1998], Lin and Vitter [1992], and Shmoys et al. [1997].

(3) See, for example, Bar-Yehuda and Even [1981], Goemans and Williamson [1995; 1997], Williamson et al. [1995], Goemans et al. [1994], Rajagopalan and Vazirani [1999], and Jain et al. [1999].

REFERENCES

AGRAWAL, A., KLEIN, P., AND RAVI, R. 1995. When trees collide: An approximation algorithm for the generalized Steiner problem on networks. SIAM J. Comput. 24, 440-456.

ARORA, S., RAGHAVAN, P., AND RAO, S. 1998. Approximation schemes for Euclidean k-medians and related problems. In Proceedings of the 30th Annual ACM Symposium on Theory of Computing (Dallas, Tex., May 23-26). ACM, New York, pp. 106-113.

BALINSKI, M. L. 1966. On finding integer solutions to linear programs. In Proceedings of the IBM Scientific Computing Symposium on Combinatorial Problems. IBM, New York, pp. 225-248.

BARTAL, Y. 1996. Probabilistic approximation of metric spaces and its algorithmic applications. In Proceedings of the 37th IEEE Symposium on Foundations of Computer Science. IEEE Computer Society Press, Los Alamitos, Calif., pp. 184-193.

BAR-YEHUDA, R., AND EVEN, S. 1981. A linear time approximation algorithm for the weighted vertex cover problem. J. Algorithms 2, 198-203.

BRADLEY, P. S., FAYYAD, U. M., AND MANGASARIAN, O. L. 1998. Mathematical programming for data mining: Formulations and challenges. Microsoft Technical Report, January.

CHARIKAR, M., CHEKURI, C., GOEL, A., AND GUHA, S. 1998. Rounding via trees: Deterministic approximation algorithms for group Steiner trees and k-median. In Proceedings of the 30th ACM Symposium on Theory of Computing (Dallas, Tex., May 23-26). ACM, New York, pp. 114-123.

CHARIKAR, M., AND GUHA, S. 1999. Improved combinatorial algorithms for the facility location and k-median problems. In Proceedings of the 40th Annual Symposium on Foundations of Computer Science. IEEE Computer Society Press, Los Alamitos, Calif., pp. 378-388.

CHARIKAR, M., GUHA, S., TARDOS, E., AND SHMOYS, D. B. 1999. A constant-factor approximation algorithm for the k-median problem. In Proceedings of the 31st Annual ACM Symposium on Theory of Computing (Atlanta, Ga., May 1-4). ACM, New York, pp. 1-10.

CHARIKAR, M., KHULLER, S., MOUNT, D. M., AND NARSHIMHAN, G. 2001. Algorithms for facility location problems with outliers. In Proceedings of the 12th Annual ACM-SIAM Symposium on Discrete Algorithms. ACM, New York, pp. 642-651.

CHUDAK, F. 1998. Improved approximation algorithms for uncapacitated facility location. In Integer Programming and Combinatorial Optimization, R. E. Bixby, E. A. Boyd, and R. Z. Rios-Mercado, eds. Lecture Notes in Computer Science; vol. 1412. Springer-Verlag, New York, pp. 180-194.

CHUDAK, F., AND SHMOYS, D. 1998. Improved approximation algorithms for the uncapacitated facility location problem. Unpublished manuscript.

CHUDAK, F., AND SHMOYS, D. 1999. Improved approximation algorithms for the capacitated facility location problem. In Proceedings of the 10th Annual ACM-SIAM Symposium on Discrete Algorithms (Baltimore, Md., Jan. 17-19). ACM, New York, pp. S875-S876.

CHUDAK, F. A., AND WILLIAMSON, D. P. 1999. Improved approximation algorithms for capacitated facility location problems. In Proceedings of the Integer Programming and Combinatorial Optimization.

CORNUEJOLS, G., NEMHAUSER, G. L., AND WOLSEY, L. A. 1990. The uncapacitated facility location problem. In Discrete Location Theory. P. Mirchandani and R. Francis, eds. Wiley, New York, pp. 119-171.

DRINEAS, P., FRIEZE, A., KANNAN, R., VEMPALA, S., AND VINAY, V. 1999. Clustering in large graphs and matrices. In Proceedings of the 10th Annual ACM-SIAM Symposium on Discrete Algorithms (Baltimore, Md., Jan. 17-19). ACM, New York, pp. 291-299.

EDMONDS, J. 1965. Maximum matching and a polyhedron with 0,1-vertices. J. Res. NBS B 69B, 125-130.

GARG, N. 1996. A 3-approximation for the minimum tree spanning k-vertices. In Proceedings of the 37th Annual IEEE Symposium on Foundation of Computer Science. IEEE Computer Society Press, Los Alamitos, Calif., pp. 302-309.

GARG, N., VAZIRANI, V., AND YANNAKAKIS, M. 1993. Primal-dual approximation algorithms for integral flow in multicut in trees, with application to matching and set cover. In Proceedings of the 20th International Colloquium on Automata, Languages and Programming.

GOEMANS, M. X., GOLDBERG, A. V., PLOTKIN, S., SHMOYS, D., TARDOS, E., AND WILLIAMSON, D. P. 1994. Improved approximation algorithms for network design problems. In Proceedings of the 5th Annual ACM-SIAM Symposium on Discrete Algorithms (Arlington, Va., Jan. 23-25). ACM, New York, pp. 223-232.

GOEMANS, M. X., AND WILLIAMSON, D. P. 1995. A general approximation technique for constrained forest problems. SIAM J. Comput. 24, 296-317.

GOEMANS, M. X., AND WILLIAMSON, D. P. 1997. The primal-dual method for approximation algorithms and its application to network design problems. In Approximation Algorithms for NP-hard Problems, D. Hochbaum, ed. PWS, pp. 144-191.

GUHA, S., AND KHULLER, S. 1998. Greedy strikes back: Improved facility location algorithms. In Proceedings of the 9th Annual ACM-SIAM Symposium on Discrete Algorithms (San Francisco, Calif., Jan.). ACM, New York, pp. 649-657.

HOCHBAUM, D. S. 1982. Heuristics for the fixed cost median problem. Math. Prog. 22, 148-162.

JAIN, K., MANDOIU, I., VAZIRANI, V. V., AND WILLIAMSON, D. P. 1999. A primal-dual schema based approximation algorithm for the element connectivity problem. In Proceedings of the 10th Annual ACM-SIAM Symposium on Discrete Algorithms (Baltimore, Md., Jan. 17-19). ACM, New York, pp. 484-489.

JAIN, K., AND VAZIRANI, V. V. 2000. An approximation algorithm for the fault tolerant metric facility location problem. In Proceedings of the 3rd Annual APPROX Conference. Lecture Notes in Computer Science, vol. 1671. Springer-Verlag, New York.

KAUFMAN, L., VANDEN EEDE, M., AND HANSEN, P. 1977. A plant and warehouse location problem. Oper. Res. Quart. 28, 547-557.

KORUPOLU, M. R., PLAXTON, C. G., AND RAJARAMAN, R. 1998. Analysis of a local search heuristic for facility location problems. In Proceedings of the 9th Annual ACM-SIAM Symposium on Discrete Algorithms (San Francisco, Calif., Jan.) ACM, New York, pp. 1-10.

KUEHN, A. A., AND HAMBURGER, M. J. 1963. A heuristic program for locating warehouses. Manage. Sci. 9, 643-666.

LIN, J.-H., AND VITTER, J. S. 1992. Approximation algorithms for geometric median problems. Inf. Proc. Lett. 44, 245-249.

LIN, J.-H., AND VITTER, J. S. 1992. E-approximation with minimum packing constraint violation. In Proceedings of the 24th Annual ACM Symposium on Theory of Computing (Victoria, B.C., Canada, May 4-6). ACM, New York, pp. 771-782.

METTU, R. R., AND PLAXTON, C. G. 2000. The online median problem. In Proceedings of the 41st Annual IEEE Symposium on Foundations of Computer Science. IEEE Computer Society Press, Los Alamitos, Calif., pp. 339-348.

NEMHAUSER, G. L., AND WOLSEY, L. A. 1990. Integer and Combinatorial Optimization. Wiley, New York.

RAJAGOPALAN, S., AND VAZIRANI, V. V. 1999. On the bidirected cut relaxation for the metric Steiner tree problem. In Proceedings of the 10th Annual ACM-SIAM Symposium on Discrete Algorithms (Baltimore, Md., Jan. 17-19). ACM, New York, pp. 742-751.

SHMOYS, D. B., TARDOS, E., AND AARDAL, K. 1997. Approximation algorithms for facility location problems. In Proceedings of the 29th Annual ACM Symposium on Theory of Computing (El Paso, Tex., May 4-6). ACM, New York, pp. 265-274.

STOLLSTEIMER, J. F. 1961. The effect of technical change and output expansion on the optimum number, size and location of pear marketing facilities in a California pear producing region. Ph.D. thesis, Univ. California at Berkeley, Berkeley, Calif.

STOLLSTEIMER, J. F. 1963. A working model for plant numbers and locations. J. Farm Econom. 45, 631-645.

VAZIRANI, V. V. 1995. Primal-dual schema based approximation algorithms. In Proceedings of the 1st Annual International Conference, COCOON. 650-652.

WILLIAMSON, D. P., GOEMANS, M. X., MIHAIL, M., AND VAZIRANI, V. V. 1995. A primal-dual approximation algorithm for generalized Steiner network problems. Combinatorica 15 (Dec.), 435-454.

RECEIVED MARCH 1999; REVISED AUGUST 2000; ACCEPTED AUGUST 2000

The research was supported by National Science Foundation (NSF) grant CCR 98-20896.

Authors' address: College of Computing, Georgia Institute of Technology, Atlanta, GA 30332-0280, e-mail: {kjain,vazirani}@cc.gatech.edu.
COPYRIGHT 2001 Association for Computing Machinery, Inc.
No portion of this article can be reproduced without the express written permission from the copyright holder.