Approximation Algorithms: Clustering and Facility Location

This is the lecture notes from Chandra Chekuri's CS583 course on Approximation Algorithms. Chapter 9: Clustering and Facility Location. You can read Chapter 10: Introduction to Network Design, here. Chapter 8: Introduction to Local Search, here.

Chapter 9

Clustering and Facility Location are two widely studied topics with a vast literature. Facility location problems have been well-studied in Operations Research and logistics. Clustering is ubiquitious with many applications in data analysis and machine learning. We confine attention to a few central problems and provide some pointers as needed to other topics. These problems have also played an important role in approximation algorithms and their study has led to a variety of interesting techniques. Research on these topics is still quite active.

For both classes of problems a key assumption that we will make is that we are working with points in some underlying metric space. Recall that a space

(V, d)

where

d : V \times V \to R_{+}

is a metric space if the distance function

d

satisfies metric properties: (i)

d (u, v) = 0

iff

u = v

(reflexivity) (ii)

d (u, v) = d (v, u)

for all

u, v \in V

(symmetry) and (iii)

d (u, v) + d (v, w) \geq d (u, w)

for all

u, v, w \in V

(triangle inequality). We will abuse the notation and use

d (A, B)

for two sets

A, B \subseteq V

to denote the quantity

min_{p \in A, q \in B} d (p, q)

. Similarly

d (p, A)

for

p \in V

and

A \subseteq V

will denote

min_{q \in A} d (p, q)

Center based clustering: In center based clustering we are given

n

points

P = {p_{1}, p_{2}, \dots, p_{n}}

in a metric space

(V, d)

, and an integer

k

. The goal is to cluster/partition

P

into

k

clusters

C_{1}, C_{2}, \dots, C_{k}

which are induced by choosing

k

centers

c_{1}, c_{2}, \dots, c_{k}

from

V

. Each point

p_{i}

is assigned to its nearest center from

c_{1}, c_{2}, \dots, c_{k}

and this induces a clustering. The nature of the clustering is controlled by an objective function that measures the quality of the clusters. Typically we phrase the problem as choosing

c_{1}, c_{2}, \dots, c_{k}

to minimize the clustering objective

\sum_{i = 1}^{n} d {(p_{i}, {c_{1}, \dots, c_{k}})}^{q}

for some

q

. The three most wellstudied problems are special cases obtained by choosing an appropriate

q

$k$ -Center is the problem when $q = \infty$ which can be equivalently phrased as $min_{c_{1}, c_{2}, \dots, c_{k} \in V} {max}_{i = 1}^{n} d (p_{i}, {c_{1}, \dots, c_{k}})$ . In other words we want to minimize the maximum distance of the input points to the cluster centers.
$k$ -Median is problem when $q = 1. min_{c_{1}, c_{2}, \dots, c_{k} \in V} \sum_{i = 1}^{n} d (p_{i}, {c_{1}, \dots, c_{k}})$ .
$k$ -Means is the problem when $q = 2. min_{c_{1}, c_{2}, \dots, c_{k} \in V} \sum_{i = 1}^{n} d {(p_{i}, {c_{1}, \dots, c_{k}})}^{2}$ .

We will mainly focus on the discrete versions of the problems where

V =

{p_{1}, p_{2}, \dots, p_{n}}

which means that the centers are to be chosen from the input points themselves. However, in many data analysis applications the points lie in

R^{d}

for some

d

and the centers can be chosen anywhere in the ambient space. In fact this makes the problems more difficult in a certain sense since the center locations now come from an infinite set. One can argue that limiting centers to the input points does not lose more than a constant factor in the approximation and this may be reasonable from a first-cut point of view but perhaps not ideal from a practical point of view. In some settings there may a requirement or advantage in choosing the cluster centers from the input set.

Facility Location: In facility location we typically have two finite sets

F

and

C

where

F

represents a set of locations where facilities can be located and

D

is a set of client/demand locations. We will assume that

V = F ⊎ D

and that there is a metrid

d

over

V

. There are several variants but one of the simplest one is the Uncapacitated Facility Location (UCFL) problem. In UCFL we are given

(F ⊎ D, d)

as well auxiliarly information which specifies the cost

f_{i}

of opening a facility at location

i \in F

. The goal is to open a subset of facilities in

F

to minimize the sum of the cost of the opened facilities and the total distance traveled by the clients to their nearest open facility. In other words we want to solve

min_{F^{'} \subseteq F} (\sum_{i \in F^{'}} f_{i} + \sum_{j \in D} d (j, F^{'}))

. The problem has close connections to

k

-Median problem. The term "uncapacitated" refers to the fact that we do not limit the number of clients that can be assigned to an open facility. In several OR applications that motivate facility location (opening warehouses or distribution facilities), capacity constraints are likely to be important. For this reasons there are several capacitated versions.

9.1 $k$ -Center

Recall that in

k

-Center we are given

n

points

p_{1}, \dots, p_{n}

in a metric space and an integer

k

and we need to choose

k

cluster centers

C = {c_{1}, c_{2}, \dots, c_{k}}

such that we minimize

max_{i} d (p_{i}, C)

. An alternative view is that we wish to find the smallest radius

R

such that there are

k

balls of radius

R

that together cover all the input points. Given a fixed

R

this can be seen as a Set Cover problem. In fact there is an easy reduction from Dominating Set to

k

-Center establishing the NP-Hardness. Moreoever, as we saw already in Chapter

1, k

-Center has no

2 - ϵ

-approximation unless

P = N P

via a reduction from Dominating Set. Here we will see two 2-approximation algorithms that are quite different and have their own advantages. The key lemma for their analysis is common and is stated below.

Lemma 9.1. Suppose there are $k + 1$ points $q_{1}, q_{2}, \dots, q_{k + 1} \in P$ such that $d (q_{i}, q_{j}) >$ $2 R$ for all $i \neq j$ . Then $OPT > R$ .

Proof. Suppose

OPT \leq R

. Then there are

k

centers

C = {c_{1}, c_{2}, \dots, c_{k}}

which induces

k

clusters

C_{1}, \dots, C_{k}

such that for each cluster

C_{h}

and each

p \in C_{h}

we have

d (p, c_{h}) \leq R

. By the pigeon hole principle some

q_{i}, q_{j}, i \neq j

are in the same cluster

C_{h}

but this implies that

d (q_{i}, q_{j}) \leq d (q_{i}, c_{h}) + d (q_{j}, c_{h}) \leq 2 R

which contradicts the assumption of the lemma.

Note that the lemma holds even if the centers can be chosen from outside the given point set

P

9.1.1 Gonzalez's algorithm and nets in metric spaces

The algorithm starts with an empty set of centers, and in each iteration picks a new center which is the farthest point from the set of centers chosen so far.

\underline{Gonzalez- k -Center (P, k)}

C \leftarrow \emptyset

2. For

i = 1

k

A. Let

c_{i} = \underset{c \in P}{\arg max} d (c, C)

C \leftarrow C \cup {c_{i}}

3. Output

C

Note that

c_{1}

is chosen arbitrarily.

Theorem 9.1. Let $R = max_{p \in P} d (p, C)$ where $C$ is the set of centers chosen by Gonzalez's algorithm. Then $R \leq 2 R^{*}$ where $R^{*}$ is the optimum $k$ -Center radius for $P$ .

Proof. Suppose not. There is a point

p \in P

such that

d (p, C) > 2 R^{*}

which implies that

p \notin C

. Since the algorithm chose the farthest point in each iteration and could have chosen

p

in each of the

k

iteration but did not, we have the property that

d (c_{i}, {c_{1}, \dots, c_{i - 1}}) > 2 R^{*}

for

i = 2

k

. This implies that the distance between each pair of points in the set

{c_{1}, c_{2}, \dots, c_{k}, p}

is more than

2 R^{*}

. By Lemma 9.1, the optimum radius must be larger than

R^{*}

, a contradiction.

Exercise 9.1. Construct an instance to demonstrate that the algorithm's worstcase performance is

2

Remark 9.1. Gonzalez's algorithm can be extended in a simple way to create a permutation of the points

P

; we simply run the algorithm with

k = n

. It is easy to see from the proof above that for any

k \in [n]

, the prefix of the permutation consisting of the first

k

points provides a 2-approximation for that choice of

k

. Thus, one can compute the permutation once and reuse it for all

k

9.1.2 Hochbaum-Shmoys bottleneck approach

A second algorithmic approach for

k

-Center is due to Hochbaum and Shmoys and has played an influential role in variants of this problem.

For a point

v

and radius

r

let

B (v, r) = {u ∣ d (u, v) \leq r}

denote the ball of radius

r

around

v

\underline{HS- k -Center (P, k)}

1. Guess

R^{*}

the optimum radius

C \leftarrow \emptyset, S \leftarrow P

3. While

(S \neq \emptyset)

A. Let

c

be an arbitrary point in

S

C \leftarrow C \cup {c}

S \leftarrow S ∖ B (c, 2 R^{*})

4. Output

C

Theorem 9.2. Let $C$ be the output of the HS algorithm for a guess $R$ . Then for all $p \in P, d (p, C) \leq 2 R$ and moreover if $R \geq R^{*}$ then $| C | \leq k$ .

Proof. The first property is easy to see since we only remove a point

p

from

S

if we add a center

c

C

such that

p \in B (c, 2 R)

. Let

c_{1}, c_{2}, \dots, c_{h}

be the centers chosen by the algorithm. We observe that

d (c_{i}, {c_{1}, \dots, c_{i - 1}}) > 2 R

. Thus, if the algorithm outputs

h

points then the pairwise distance between any two of them is more than

2 R

. By Lemma 9.1, if

h \geq k + 1

the optimum radius is

> R

. Hence, if the guess

R \geq R^{*}

the algorithm outputs at most

k

centers.

The guessing of

R^{*}

can be implemented by binary search in various ways. We omit these routine details.

Exercise 9.2. Describe an example where the algorithm uses exactly

k

centers even with guess

R^{*}

. Describe an example where the algorithm outputs less than

k

centers with a guess of

R^{*}

The

k

-Center problem is natural in geometric settings and one can see from the proof that the 2-approximation for the two algorithms holds even when allowing for centers to be chosen outside. A surprising result of Feder and Greene ^[1] shows that even in two dimensions (the Euclidean plane) one cannot improve the factor of 2 unless

P = N P

The

k

-Supplier problem is closely related to

k

-Center and is motivated by facility location considerations. Here we are given

F \cup D

which are in a metric space. We need to choose

k

centers

C \subseteq F

to minimize

max_{p \in D} d (p, C)

. Note that we don't have to cover the facilities. Hochbaum and Shmoys gave a variant of their algorithm that obtains a 3-approximation for k-Supplier and moreover showed that unless

P = N P

one cannot improve 3 ^[2]. Interestingly in Euclidean spaces 3 is not tight ^[3]. Several generalizations of

k

-Center which constrain the choice of centers have been considered - see ^[4] for a general framework that also considers outliers.

One can consider weighted version of

k

-Center or relabeled as priority version in subsequent work. We refer to the work of Plesnik ^[5] and a recent one ^[6] on this variant which has found applications in fair clustering.

We finally mention that

k

-Center clustering has nice connections to the notion of

r

-nets in metric spaces. Given a set of points

P

in a metric space and a radius

r

, an

r

-net is a set of centers

C

such that (i) for all

p \in P

we have

d (p, C) \leq r

(that is the points are covered by balls of radius

r

) and (ii) for any distinct

c, c^{'} \in C

we have

B (c, r / 2)

and

B (c^{'}, r / 2)

are disjoint (packing property or the property that no two centers are too close).

r

-nets provide a concise summary of distances in a metric space at scale

r

. One can use

r

-nets to obtain nearest-neighbor data structures and other applications, especially in low-dimensional settings. We refer the reader to ^[7],^[8].

LP Relaxation: The two

k

-Center algorithms we described are combinatorial. One can also consider an LP relaxation. Since it is a bottleneck problem, writing an LP relaxation involves issues similar to what we saw with unrelated machine scheduling. Given a guess

R

we can write an LP to check whether a radius

R

is feasible and then find the smallest

R

for which it is feasible. The feasibility LP can be written as follows. Let

x_{i}

be an indicator for whether we open a center at point

p_{i}

\begin{aligned} \sum_{i = 1}^{n} x_{i} & = k \\ \sum_{p_{i} \in B (p_{j}, R)} x_{i} & \geq 1 \forall p_{j} \in P \\ x_{i} & \geq 0 \forall p_{i} \in P \end{aligned}

Exercise 9.3. Prove that if

R

is feasible for the preceding LP then one can obtain a solution with

k

centers with max radius

2 R

Exercise 9.4. Generalize the LP for the

k

-Supplier problem and prove that one can obtain a 3-approximation with respect to lower bound provided via the LP approach.

9.2 Uncapacitated Facility Location

We now discuss UCFL. One can obtain a constant factor for UCFL via several techniques: LP rounding, primal-dual, local search and greedy. The best known approximation bound is

1.488

due to Li ^[9] while it is known that one cannot obtain a ratio better than

1.463

^[10]. We will describe the complicated algorithms and focus on the high-level approaches that yield some constant factor approximation.

It is common to use

i

to denote a facility in

F

and

j

to denote a demand/client.

9.2.1 LP Rounding

The first constant factor approximation for UCFL was via LP rounding by Aardal, Shmoys and Tardos using a filtering technique of Lin and Vitter. We start with the LP relaxation. We use a variable

y_{i}

for

i \in F

to indicate whether

i

is opened or not. We use a variable

x_{i, j}

to indicate whether

j

is connected to

i

(or assigned to

i

). One set of constraints are natural here: each client has to be assigned/connected to a facility. The other constraint requires that

j

is assigned to

i

only if

i

is open.

\begin{aligned} min \sum_{i \in F} f_{i} y_{i} + \sum_{j \in D} \sum_{i \in F} d (i, j) x_{i, j} \\ \sum_{i} x_{i, j} & = 1 \forall j \in D \\ x_{i, j} & \leq y_{i} i \in F, j \in D \\ x, y & \geq 0 \end{aligned}

Given a feasible solution

x, y

to the LP the question is how to round. We note that the LP relaxation does not "know" whether

d

is a metric or not. In fact when

d

is arbitrary (but non-negative) we obtain the non-metric facility location problem which is as hard as the Set Cover problem but not much harder - one can obtain an

O (\log n)

-approximation. However, we can obtain a constant factor when

d

is a metric and the rounding will exploit this property.

Given the fractional solution

x, y

for each

j

we define

α_{j}

to be the distance cost paid for by the LP: therefore

α_{j} = \sum_{i \in F} d (i, j) x_{i, j}

. Note that the LP cost is

\sum_{i} f_{i} y_{i} + \sum_{j} α_{j}

Lemma 9.2. For each $j$ and each $δ \in (0, 1)$ there is a total facility value of at least $(1 - δ)$ in $B (j, α_{j} / δ)$ . That is, $\sum_{i \in B (j, α_{j} / δ)} y_{i} \geq 1 - δ$ . In particular $\sum_{i \in B (j, 2 α_{j})} y_{i} \geq 1 / 2$ .

Proof. This essentially follows from Markov's inequality or averaging. Note that

α_{j} = \sum_{i} d (i, j) x_{i, j}

and

\sum_{i} x_{i, j} = 1

. Suppose

\sum_{i \in B (j, α_{j} / δ)} y_{i} < 1 - δ

. Since

x_{i, j} \leq y_{i}

for all

i

j

, we will have

α_{j} > δ α_{j} / δ

which is impossible.

We say that two clients

j

and

j^{'}

intersect if there is some

i \in F

such that

i \in B (j, 2 α_{j}) \cap B (j^{'}, 2 α_{j^{'}})

. The rounding algorithm is described below.

\underline{UCFL-primal-rounding}

1. Solve LP and let

(x, y)

be a feasible LP solution

2. For each

j

compute

α_{j} = \sum_{i} d (i, j) x_{i, j}

3. Renumber clients such that

α_{1} \leq α_{2} \leq \dots \leq α_{h}

where

h

is number of clients

4. For

j = 1

h

A. If

j

already assigned continue

B. Open cheapest facility

i

B (j, 2 α_{j})

and assign

j

i

C. For each remaining client

j^{'} > j

that intersects with

j

, assigne

j^{'}

i

5. Output the list of open facilities and the client assignment

It is not hard to see that every client is assigned to an open facility. The main issue is to bound the total cost. Let

F

be the total facility opening cost, and let

C

be the total connection cost. We will bound these separately.

Lemma 9.3.

F \leq 2 \sum_{i} f_{i} y_{i}

Proof. Note that a client

j

opens a new facility only if it has not been assigned when it is considered by the algorithm. Let

j_{1}, j_{2}, \dots, j_{k}

be the clients that open facilities. Let

A_{j} = F \cap B (j, 2 α_{j})

be the set of facilities in the ball of radius

2 α_{j}

around

j

. From the algorithm and the definition of intersection of clients, we see that the sets

A_{j_{1}}, A_{j_{2}}, \dots, A_{j_{k}}

are pairwise disjoint. The algorithm opens the cheapest facility in

A_{j_{ℓ}}

for

1 \leq ℓ \leq k

. Note that

\sum_{i \in A_{j ℓ}} y_{i} \geq 1 / 2

by Lemma 9.2. Hence the cost of the cheapest facility in

A_{j ℓ}

is at most

2 \sum_{i \in A_{j}} f_{i} y_{i}

(why?). By the disjointness of the sets

A_{j_{1}}, \dots, A_{j_{k}}

we see that the total cost of the facilities opened is at most

2 \sum_{i} f_{i} y_{i}

Lemma 9.4.

C \leq 6 \sum_{j} α_{j}

Proof. Consider a client

j

that opens a facility

i

B (j, 2 α_{j})

. In this case

j

is assigned to

i

and

d (i, j) \leq 2 α_{j}

. Now consider a client

j

that does not open a facility. This implies that there was a client

j^{'} < j

that opened a facility

i^{'}

and

j^{'}

and

j

intersect, and

j

was assigned to

i^{'}

. What is

d (i^{'}, j)

? We claim that

d (i^{'}, j) \leq 6 α_{j}

. To see this we note that

j^{'}

and

j

intersect because there is some facility

i \in B (j, 2 α_{j^{'}}) \cap B (j, 2 α_{j})

. By considering the path

j \to i \to j^{'} \to i^{'}

, via triangle inequality,

d (i^{'}, j) \leq 2 α_{j} + 2 α_{j^{'}} + 2 α_{j^{'}} \leq 6 α_{j}

since

α_{j^{'}} \leq α_{j}

. Thus the distance traveled by each client

j

to its assigned facility is at most

6 α_{j}

. The lemma follows.

The two preceding lemmas give us the following which implies that the algorithm yields a 6-approximation.

Theorem 9.3.

C + F \leq 6 {OPT}_{L P}

It should be clear to the reader that the algorithm and analysis are not optimized for the approximation ratio. The goal here was to simply outline the basic scheme that led to the first constnat factor approximation.

9.2.2 Primal-Dual

Jain and Vazirani developed an elegant and influential primal-dual algorithm for UCFL ^[11]. It was influential since it allowed new algorithms for

k

-median and several generalizations of UCFL in a clean way. Moreover the primal-dual algorithm is simple and efficient to implement. On the other hand we should mention that one advantage of having an LP solution is that it gives an explicit lower bound on the optimum value while a primal-dual method yields a lower bound via a feasible dual which may not be optimal. We need some background to describe the primal-dual method in approximation.

Complementary slackness: To understand primal-dual we need some basic background in complementary slackness. Suppose we have a primal

LP (P)

of the form

min c x

s.t

A x \leq b, x \geq 0

which we intentionally wrote in this standard form as a covering LP. It's

dual (D)

is a packing LP

max b^{t} y

s.t

y A^{t} \geq c, y \geq 0

. We will assume that both primal and dual are feasible and hence the optimum values are finite, and via strong duality we know that the optimum values are the same.

Definition 9.4. A feasible solution $x$ to $(P)$ and a feasible solution $y$ to $(D)$ satisfy the primal complementary slackness condition with respect to each other if the following is true: for each $i, x_{i} = 0$ or the corresponding dual constraint is tight, that is $\sum_{j} A_{i, j} y_{j} = c_{i}$ . They satisfy the dual complementary slackness condition if the following is true: for each $j, y_{j} = 0$ or the corresponding primal constraint is tight, that is $\sum_{j} A_{j, i} x_{i} = b_{j}$ .

One of the consequences of the duality theory of LP is the following.

Theorem 9.5. Suppose $(P)$ and $(D)$ are a primal and dual pair of LPs that both have finite optima. A feasible solution $x$ to $(P)$ and a feasible solution $y$ to $(D)$ satisfy the primal-dual complementary slackness properties with respect to each other if and only if they are both optimum solutions to the respective LPs.

We illustrate the use of complementary slackness in the context of approximation via Vertex Cover. Recall the primal covering LP and we also write the dual.

\begin{aligned} min \sum_{v \in V} w_{v} x_{v} \\ x_{u} + x_{v} & \geq 1 \forall e = (u, v) \in E \\ x_{v} & \geq 0 \forall v \in V \end{aligned}

\begin{aligned} max \sum_{e \in E} y_{e} \\ \sum_{e \in δ (v)} y_{e} & \leq w_{v} \forall v \in V \\ y_{e} & \geq 0, \forall e \in E \end{aligned}

Recall that we described a simple rounding algorithm. Given a feasible primal solution

x

. We let

S = {v ∣ x_{v} \geq 1 / 2}

and showed that (i)

S

is a vertex cover for the given graph

G

and (ii)

w (S) \leq 2 \sum_{v} w_{v} x_{v}

. Now suppose we have an optimum solution

x^{*}

to the primal rather than an arbitrary feasible solution. We can prove an interesting and stronger statement via complementary slackness.

Lemma 9.5. Let $x^{*}$ be an optimum solution to the primal covering LP. Then $S = {v ∣$ $x_{v}^{*} > 0}$ is a feasible vertex cover for $G$ and moreover $w (S) \leq 2 \sum_{v} w_{v} x_{v}^{*}$ .

Proof. It is easy to see that

S

is a vertex cover via the same argument that we have seen before. How do we bound the cost now that we may be rounding

x_{v}^{*}

to 1 even though

x_{v}^{*}

may be tiny? Let

y^{*}

be any optimum solution to the dual; one exists (why?). Via strong duality we have

\sum_{v} w_{v} x_{v}^{*} = \sum_{e} y_{e}^{*}

. Via primal-complementary slackness we have the property that if

x_{v}^{*} > 0

then

\sum_{e \in δ (v)} y_{e}^{*} = w_{v}

. Hence

w (S) = \sum_{v : x_{v}^{*} > 0} w_{v} = \sum_{v : x_{v}^{*} > 0} \sum_{e \in δ (v)} y_{e}^{*} .

Interchanging the order of summation we obtain that

w (S) = \sum_{v : x_{v}^{*} > 0} \sum_{e \in δ (v)} y_{e}^{*} \leq 2 \sum_{e \in E} y_{e}^{*} = 2 \sum_{v} w_{v} x_{v}^{*}

where we use the fact that an edge

e

has only two end points in the inequality.

Primal-dual for approximating Vertex Cover: We will first study a primaldual approximation algorithm for the Vertex Cover problem - this algorithm due to Bar-Yehuda and Even ^[12] can perhaps be considered the first primal-dual algorithm in approximation. Primal-dual is a classical method in optimization and is applicable in both continuous and discrete settings. The basic idea, in the context of LPs (the method applies more generally), is to obtain a solution

x

to a primal LP and a solution

y

to the dual LP together and certify the optimality of each solution via complementary slackness. It is beyond the scope of this notes to give a proper treatment. One could argue that understanding the setting of approximation is easier than in the classical setting where one seeks exact algorithms, since our goals are more modest.

Typically one starts with one of

x, y

being feasible and the other infeasible, and evolve them over time. In discrete optimization, this method is successfully applied to LPs that are known to be integer polytopes. Examples include shortest paths, matchings, and others. In approximation the LP relaxations are typically not integral. In such a setting the goal is to produce a primal and dual pair

x, y

where

x

is an integral feasible solution, and

y

is fractional feasible solution. The goal is to approximately bound the value of

x

via the dual

y

, and for this purpose we will enforce only the primal complementary slackness condition for the pair

(x, y)

. To make the algorithm and analysis manageable, the primal-dual algorithm is typically done in a simple but clever fashion – there have been several surprisingly strong and powerful approximation results via this approach.

We illustrate this in the context of Vertex Cover first. It is useful to interpret the dual as we did already in the context of the dual-fitting technique for Set Cover. We think of the edges

e \in E

as agents that wish to be covered by the vertices in

G

at minimum cost. The dual variable

y_{e}

can be thought of as the amount that edge

e

is willing to pay to be covered. The dual packing constraint

\sum_{e \in δ (v)} y_{e} \leq w_{v}

says that for any vertex

v

, the total payments of all edges incident to

v

cannot exceed its weight. This can be understood game theoretically as follows. The set of edges

δ (v)

can together pay

w_{v}

and get covered, and hence we cannot charge them more. The dual objective is to maximize the total payment that can be extracted from the edges subject to these natural and simple constraints. With this interpretation in mind we wish to produce a feasible dual (payments) and a corresponding feasible integral primal (vertex cover). The basic scheme is to start with an infeasible primal

x = 0

and a feasible dual

y = 0

and increase

y

while maintaining feasibility; during the process we will maintain primal complementary slackness which means that if a dual constraint for a vertex

v

becomes tight we set

x_{v} = 1

. Note that we are producing an integer primal solution in this process. How should we increase

y

values? We will do it in the naive fashion which is to uniformly increase

y_{e}

for all

e

that are not already part of a tight constraint (and hence not covered yet).

\underline{VC-Primal-Dual (G = (V, E), w : V \to R_{+})}

x \leftarrow 0, y \leftarrow 0 // initialization: primal infeasible, dual feasible

U \leftarrow E // uncovered edges that are active

3. While

(U \neq \emptyset)

A. Increase

y_{e}

uniformly for each

e \in U

until constraint

\sum_{e \in δ (v)} y_{e} = w_{v}

for

some vertex

a

B. Set

x_{a} = 1 // Maintain primal complementary slackness

U \leftarrow U ∖ δ (a) // Remove all edges covered by a

4. Output integer solution

x

and dual certificate

y

Remark 9.2. Note that when checking whether a vertex

v

is tight we count the payments from

e \in δ (v)

even though some of them are no longer active.

Lemma 9.6. At the end of the algorithm $x$ is a feasible vertex cover for $G$ and $\sum_{v} w_{v} x_{v} \leq 2 OPT$ .

Proof. By induction on the iterations one can prove that (i)

y

remains dual feasible throughout (ii)

a b \in U

at the start of an iteration iff

e = a b

is not covered yet (iii) each iteration adds at least one more vertex and hence the algorithm terminates in

\leq n

iterations and outputs a feasible vertex cover. The main issue is the cost of

x

For this we note that the algorithm maintains primal complementary slackness. That is,

x_{v} = 0

or if

x_{v} = 1

then

\sum_{e \in δ (v)} y_{e} = w_{v}

. Thus, we have

\sum_{v} w_{v} x_{v} = \sum_{v : x_{v} > 0} \sum_{e \in δ (v)} y_{e} \leq 2 \sum_{e} y_{e} \leq 2 {OPT}_{L P} .

We used the fact that

e

has at most two end points in the first inequality and the fact that

y

is dual feasible in the second inequality. In terms of payment what this says is that edge

u v

's payment of

y_{u v}

can be used to pay for opening

u

and

v

while the dual pays only once.

As the reader can see, the algorithm is very simple to implement.

Exercise 9.5. Describe an example to show that the primal-dual algorithm's worst case performance is 2. Describe an example to show that the dual value constructed by the algorithm is

≃ OPT / 2

. Are these two parts the same?

Remark 9.3. The algorithm generalizes to give an

f

-approximation for Set Cover where

f

is the maximum frequency of any element. There are examples showing that the performance of this algorithm in the worst-case can indeed be a factor of

f

. We saw earlier that the integrality gap of the LP is at most

1 + \ln d

where

d

is the maximum set size. There is no contradiction here since the specific primal-dual algorithm that we developed need not achieve the tight integrality gap.

Primal-Dual for UCFL: Now we consider a primal-dual algorithm for UCFL. Recall the primal LP that we discussed previously. Now we describe the dual LP written below. The dual has two types of variables. For each client

j

there is a variable

α_{j}

corresponding to the primal constraint that each client

j

needs to be connected to a facility. For each facility

i

and client

j

there is a variable

β_{i, j}

corresponding to the constraint that

y_{i} \geq x_{i, j}

\begin{aligned} max \sum_{j \in D} α_{j} \\ \sum_{i} β_{i, j} & \leq f_{i} \forall i \in F \\ α_{j} - β_{i, j} & \leq d (i, j) i \in F, j \in D \\ α, β & \geq 0 \end{aligned}

It is important to interpret the dual variables. There is a similarity to Set Cover since Non-Metric Facility Location is a generalization, and the LP formulation does not distinguish between metric and non-metric facility location - it is only in the rounding that we can take advantage of metric properties. The variable

α_{j}

can be interpreted as the amount of payment client

j

is willing to make. This comes in two parts - the payment to travel to a facility which it cannot share with any other clients, and the payment it is willing to make to open a facility which it can share with other clients. The variable

β_{i, j}

corresponds to the amount client

j

is willing to pay to facility

i

to open it. (In Set Cover there is no need to distinguish between

α_{j}

and

β_{i, j}

since there are no distances (or they can be assumed to be 0 or

\infty

).) The first set of constraints in the dual say that for any facility

i

, the total payments from all the clients

(\sum_{j} β_{i, j})

cannot exceed

cost f_{i}

. The second set of constraints specify that

α_{j} - β_{i, j}

is at most

d (i, j)

. One way to understand this is that if

α_{j} < d (i, j)

then client

j

will not even be able to reach

i

and hence will not contribute to opening

i

. Via this interpretation it is convenient to assume that

β_{i, j} = max {0, α_{j} - d (i, j)}

The primal-dual algorithm for UCFL will have a growth stage that is similar to what we saw for Vertex Cover. We increase

α_{j}

for each "uncovered" client

j

uniformly. A facility

i

will receive payment

β_{i, j} = max {0, α_{j} - d (i, j)}

from

j

. To maintain dual feasibility, as soon as the constraint for facility

i

becomes tight we open facility

i

(in the primal we set

y_{i} = 1

); note that any client

j

such that

β_{i, j} > 0

cannot increase its

α_{j}

any further and hence will stop growing since it is connected to an open facility. This process is very similar to that in Set Cover. The main issue is that we will get a weak approximation (a factor of

| F |

) since a client can contribute payments to a large number of facilities. Note that the process so far has not taken advantage of the fact that we have a metric facility location problem. Therefore, in the second phase we will close some facilities which means that a client may need to get connected to a facility that it did not contribute to - however we will use the metric property to show that a client does not need to travel too far to reach a facility that remains open.

With the above discussion in place, we describe the two phase primal-dual algorithm below. The algorithm also creates a bipartite graph

G

with vertex set

F \cup D

and initially it has no edges.

\underline{JV-primal-dual ((F \cup D, d), f_{i}, i \in F)}

α_{j} = 0

for all

j \in D, β_{i, j} \leftarrow 0

for all

i, j

initialize dual values to 0

A \leftarrow D

active clients that are unconnected

O \leftarrow \emptyset

temporarily opened facilities

4. While

(A \neq \emptyset)

Growth phase

A. Increase

α_{j}

uniformly for each

j \in A

while maintaining the invariant

max {0, α_{j} - d (i, j)} = β_{i, j}

until one of following conditions hold: (i)

α_{j} = d (i, j)

for some

j \in A

and

i \in O

or (ii)

\sum_{j} β_{i j} = f_{i}

for some

i \in F ∖ O

B. If (i) happens then add edge

(i, j)

G

and and

A \leftarrow A - {j}

j is

connected to already open facility i \in O

C. Else If (ii) happens then

O \leftarrow O \cup {i}

temporarily open i that became tight

2. for each

j \in D

such that

β_{i, j} > 0

add edge

(i, j)

G

note: clients

not in A may also get edges to i

A \leftarrow A ∖ {j : β (i, j) > 0}

make active clients connected to i inactive

5. Create graph

H

with vertex set

O

and edge set

Q

where

(i, i^{'}) \in Q

iff

there exists client

j

such that

β (i, j) > 0

and

β (i^{'}, j) > 0

O^{'}

is any maximal independent set in

H

7. Close facilities in

O ∖ O^{'}

Pruning phase

8. Assign each client

j

to nearest facility in

O^{'}

Example: The example in Fig 9.2.2 illustrates the need for the pruning phase. There are

2 n

clients and

n

facilities and the opening costs of the facilities are

n + 2

except for the first one which has an opening cost of

n + 1

. The first group of clients shown at the top of the figure have a connection cost of 1 to each facility. The second group of clients have the following property:

d (i_{h}, j_{h}^{'}) = 1

and

d (i_{ℓ}, j_{h}^{'}) = 3

when

ℓ \neq j

. The rest of the distances are induced by these. One can verify that in the growth phase all the facilities will be opened. However the total dual value after the growth phase is

5 n - 1

while the total cost of the opened facilitie is

Θ (n^{2})

and hence pruning is necessary to obtain a good approximation.

Figure 9.1: Example to illustrate the need for the pruning phase.

We will do a high-level analysis of the algorithm skipping a few fine details. A formal treatment with full details can be found in ^[13].

The algorithm maintains the property that

α_{j}

and

β_{i, j}

variables are dual feasible. Consider the graph

G

created by the algorithm. We call an edge

(i, j) \in G

a witness edge for

j

when

j

is removed from

A

(it happens exactly once for each

j

). We call an edge

(i, j) \in G

special if

β (i, j) > 0

which means that

j

paid to temporarily open

i

. We remark that a client

j

may have a special edge to

i

even though

(i, j)

is not its witness edge - this can happen if

i

is temporarily opened after

j

was already removed from

A

(due to other clients contributing to opening

i

later). One can associate a notion of time with the progression of the algorithm since dual variables are monotonically increased. This can be used to order events.

A basic property maintained by the algorithm is the following.

Claim 9.2.1. If $(i, j)$ is an edge in $G$ then $α_{j} - β_{i, j} = d (i, j)$ and hence $α_{j} \geq d (i, j)$ .

Proof. The algorithm adds an edge

(i, j)

only if

(i, j)

is a witness edge for

j

or if

β (i, j) > 0

(in which case

(i, j)

is special). The algorithm maintain the invariant

β_{i, j} = max {0, α_{j} - d (i, j)}

and hence if

β_{i, j} > 0

the claim is clear. If

β_{i, j} = 0

then

(i, j)

is a witness edge for

j

and case (i) happens when

j

is removed from

A

and in this case

α_{j} = d (i, j)

Analysis: We upper bound the cost of opening the facilities in

O^{'}

and the connection cost of the clients.

We leave the proof of the following lemma as an exercise.

Lemma 9.7. For each facility $i \in O$ we have $\sum_{(i, j) is special} β_{i, j} = f_{i}$ .

Since a client

j

can pay for multiple facilities which we cannot afford, the pruning phase removes facilities such that no client

j

is connected two facilities in

O^{'}

with special edges (otherwise those two facilities will have an edge in

H

). We say that a client

j

is directly connected to a facility

i \in O^{'}

(i, j)

is a special edge. We call all such clients directly connected clients and the rest of the clients are called indirectly connected. We let

D_{1} = \cup_{i \in O^{'}} Z_{i}

be the set of all directly connected clients and let

D_{2}

be the set of all indirectly connected clients.

For

i \in O^{'}

let

Z_{i}

be the directly connected clients. By the pruning rule we have the property that a client

j

is directly connected to at most one facility in

O^{'}

. We show that each facility in

O^{'}

can be paid for by its directly connected clients.

Lemma 9.8. For each

i \in O^{'}, \sum_{j \in Z_{i}} β_{i, j} = f_{i}

Proof. From Lemma 9.7 and the fact that if

i \in O^{'}

then every client

j

with special edge

(i, j)

must be directly connected to

i

From Claim 9.2.1 we see that if

j

is directly connected to

i

then

α_{j} - β_{i, j} = d (i, j)

, and hence

j

can pay its connection cost to

i

and its contribution to opening

i

. What about indirectly connected clients? The next lemma bounds their connection cost.

Lemma 9.9. Suppose $j \in D_{2}$ , that is, it is an indirectly connected client. Let $i$ be its closest facility in $O^{'}$ then $d (i, j) \leq 3 α_{j}$ .

Proof. Let

(i, j)

be the witness edge for

j

. Note that

i \in O

. Since

j

is an indirectly connected client there is no facility

i^{'} \in O^{'}

such that

(i^{'}, j)

is a special edge. Since

i \notin O^{'}

it must be because

i

was closed in the pruning phase and hence there must be a facility

i^{'} \in O^{'}

such that

(i, i^{'})

is an edge in

H

(otherwise

O^{'}

would not be a maximal independent set). Therefore there is some client

j^{'} \neq j

such that

(i^{'}, j^{'})

and

(i, j^{'})

are both special edges. We claim that

α_{j} \geq α_{j^{'}}

. Assuming this claim we see via triangle inequality and Claim 9.2.1 that,

d (i^{'}, j) \leq d (i, j) + d (i, j^{'}) + d (i^{'}, j^{'}) \leq α_{j} + 2 α_{j^{'}} \leq 3 α_{j} .

Since

i^{'} \in O^{'}

the nearest client to

j

is within distance

\leq 3 α_{j}

We now prove the claim that

α_{j} \geq α_{j^{'}}

. Let

t = α_{j}

be the time when

j

connects to facility

i

as its witness. Consider two cases. In the first case

d (i, j^{'}) \leq t

which means that

j^{'}

had already reached

i

at or before

t

; in this case

α_{j^{'}} \leq t

since

i

was open at

t

. In the second case

d (i, j^{'}) > t

; this means that

j^{'}

reached

i

strictly after

t

. Since

i

was already open at

t, j^{'}

would not pay to open

i

which implies that

β (i, j^{'}) = 0

but then

(i, j^{'})

would not be a special edge and hence this case cannot arise. This finishes the proof of the claim.

With the preceding two lemmas in place we can bound the total cost of opening facilities in

O^{'}

and connecting clients to them. We will provide a refined statement that turns out to be useful in some applications.

Theorem 9.6.

\sum_{j \in D} d (O^{'}, j) + 3 \sum_{i \in O^{'}} f_{i} \leq 3 \sum_{j \in D} α_{j} \leq 3 {OPT}_{L P} .

In particular the algorithm yields a $3$ -approximation.

Proof. Consider directly connected clients

D_{1}

. We have

D_{1} = ⊎_{i \in O^{'}} Z_{i}

where

Z_{i}

are directly connected to

i

. Via Lemma 9.8 and Claim 9.2.1

\begin{aligned} \sum_{j \in D_{1}} α_{j} & = \sum_{i \in O^{'}} \sum_{j \in Z_{i}} α_{j} \\ = \sum_{i \in O^{'}} \sum_{j \in Z_{i}} (d (i, j) + β_{i, j}) \\ = \sum_{j \in O^{'}} (f_{i} + \sum_{j \in Z_{i}} d (i, j)) \\ \geq \sum_{j \in O^{'}} f_{i} + \sum_{j \in D_{1}} d (O^{'}, j) . \end{aligned}

For indirectly connected clients, via Lemma 9.9, we have

3 \sum_{j \in D_{2}} α_{j} \geq

\sum_{j \in D_{2}} d (O^{'}, j)

. Thus

\begin{aligned} 3 \sum_{j \in D} α_{j} & = 3 \sum_{j \in D_{1}} α_{j} + 3 \sum_{j \in D_{2}} α_{j} \\ \geq 3 \sum_{i \in O^{'}} f_{i} + 3 \sum_{j \in D_{1}} d (O^{'}, j) + \sum_{j \in D_{2}} d (O^{'}, j) \\ \geq 3 \sum_{i \in O^{'}} f_{i} + \sum_{j \in D} d (O^{'}, j) . \end{aligned}

The algorithm is easy and efficient to implement. One of the main advantages of the stronger property that we saw in the theorem is that it leads to a nice algorithm for the

k

-median problem; we refer the reader to Chapter 25 in ^[13:1] for a detailed description. In addition the flexibility of the primal-dual algorithm has led to algorithms for several other variants; see ^[14] for one such example.

9.2.3 Local Search

Local search has been shown to be a very effective heuristic algorithm for facility location and clustering probelms and there is extensive literature on this topic. The first paper that proved constant factor approximation bounds for UCFL is by Korupolu, Plaxtion and Rajaraman ^[15] and it provided a useful template for many future papers. We refer the reader to Chapter 9 in ^[16] for local search analysis for UCFL.

9.3 $k$ -Median

k

-Median has been extensively studied in approximation algorithms due to its simplicity and connection to UCFL. The first constant factor approximation was obtained in ^[17] via LP rounding. We will consider a slight generalization of

k

-Median where the medians are to be selected from the facility set

F

. We describe the LP which is closely related to that for UCFL which we have already seen. The variables are the same:

y_{i}

indicates whether a center is opened at location

i \in F

and

x_{i, j}

indicates whether client

j

is connected to facility

i

. The objective and constraints change since the problem requires one to open at most

k

facilities but there is no cost to opening them.

\begin{aligned} min \sum_{j \in D} \sum_{i \in F} d (i, j) x_{i, j} \\ \sum_{i} x_{i, j} & = 1 \forall j \in D \\ x_{i, j} & \leq y_{i} i \in F, j \in D \\ \sum_{i} y_{i} & \leq k \\ x, y & \geq 0 \end{aligned}

Rounding of the LP for

k

-Median is not as simple as it is for UCFL. This is mainly because one needs a global argument. To understanding this consider the following example. Suppose we have

k + 1

points where the distance between each pair is

1

(uniform metric space). Then the optimum integer solution has cost

1

since we can place

k

medians at

k

points and the remaining point has to travel a distance of

1

. Now consider a fractional solution where

y_{i} = (1 - 1 / (k + 1))

for each point. The cost of this fractional solution is also

1

, however, each client now pays a fractional cost of

1 / (k + 1)

. Any rounding algorithm will make at least one of the clients pay a cost of

1

which is much larger than its fractional cost; thus the analysis cannot be based on preserving the cost of each client within a constant factor of its fractional cost; thus the analysis cannot be based on preserving the cost of each client within a constant factor of its fractional cost.

There are several LP rounding algorithms known. An advantage of LP based approach is that it leads to a constant factor approximation for the Matroid Median problem which is a nice and powerful generalization of the

k

-Median problem; here there is a matroid defined over the facilities and the constraint is that the set of facilities chosen must be independent in the given matroid. One can write a natural LP relaxation for this problem and prove that the LP has a constant integrality gap by appealing to matroid intersection! It showcases the power of classical result in combinatorial optimization. We refer the reader to ^[18],^[19].

9.3.1 Local Search

Local search for center based clustering problems is perhaps one of the most natural heuristics. In particular we will consider the

p

-swap heuristic. Given the current set of

k

centers

S

, the

p

-swap heuristic will consider swapping out up to

p

centers from

S

with

p

new centers. It is easy to see that this local search algorithm can be implemented in

n^{O (p)}

time for each iteration. When

p = 1

we simply refer to the algorithm as (basic) local search. We will ignore the convergence time. As we saw for Max Cut, one can use standard tricks to make the algorithm run in polynomial time with a loss of a

(1 + o (1))

-factor in the approximation bound guaranteed by the worst-case local optimum. Thus the main focus will be on the quality of the local optimum. The following is a very nice result.

Theorem 9.7 (Arya et al. ^[20]). For any fixed $p \geq 1$ the $p$ -swap local search heuristic has a tight worst-case approximation ratio of $(3 + 2 / p)$ for $k$ -Median. In particular the basic local search algorithm yields a $5$ -approximation.

Here we give a proof/analysis of the preceding theorem for

p = 1

, following the simplified analysis presented in ^[21]. See also ^[16:1] and the notes from CMU. Given any set of centers

S

we define

cost (S) = \sum_{j \in D} d (j, S)

to be the sum of the distances of the clients to the centers. Let

S

be a local optimum and let

S^{*}

be some fixed optimum solution to the given

k

-Median instance. To compare

cost (S)

with

cost (S^{*})

the key idea is to set up a clever set of potential swaps between the centers in

S

and centers in

S^{*}

. That is, we consider a swap pair

(r, f)

where

r \in S

and

f \in S^{*}

. Since

S

is a local optimum it must be the case that

cost (S - r + f) \leq cost (S)

. The analysis upper bounds the potential increase in the cost in some interesting fashion and sums up the resulting series of inequalities. This may seem magical, and indeed it is not obvious why the analysis proceeds in this fashion. The short answer is that the analysis ideas required a series of developments with the somewhat easier case of UCFL coming first.

We set up some notation. Let

ϕ : D \to S

be the mapping of clients to facilities in

S

based on nearest distance. Similarly let

ϕ^{*} : D \to S^{*}

the mapping to facilities in the optimum solution

S^{*}

. Thus

j

connects to facility

ϕ (j)

in the local optimum and to

ϕ^{*} (j)

in the optimum solution. We also let

N (i)

denote the set of all clients assigned to a facility

i \in S

and let

N^{*} (i)

denote the set of all clients assigned to a facility

i \in S^{*}

. Let

A_{j} = d (j, S)

and

O_{j} = d (j, S^{*})

be the cost paid by

j

in local optimum and optimal solutions respectively. To reinforce the notation we express

cost (S)

as follows.

cost (S) = \sum_{j \in D} A_{j} = \sum_{i \in S} \sum_{j \in N (i)} d (j, i) .

Similarly

cost (S^{*}) = \sum_{j \in D} O_{j} = \sum_{i \in S^{*}} \sum_{j \in N^{*} (i)} d (j, i) .

Setting up the swap pairs: We create a set of pairs

P

as follows. We will assume without loss of generality that

| S | = | S^{*} | = k

. For technical convenience we also assume

S \cap S^{*} = \emptyset

; we can always create dummy centers that are co-located to make this assumption. Consider the mapping

ρ : S^{*} \to S

where each

i \in S^{*}

is mapped to its closest center in

S

; hence

ρ (i)

for

i \in S^{*}

is the closest center in

S

to it. Let

R_{1}

be the set of centers in

S

that have exactly one center in

S^{*}

mapped to them. Let

R_{0}

be the set of centers in

S

with no centers of

S^{*}

mapped to them. This means that for each

i \in S ∖ (R_{0} \cup R_{1})

there are two or more centers mapped to them. Let

S_{1}^{*} \subseteq S^{*}

be the centers mapped to

R_{1}

. See Figure 9.2. By a simple averaging argument we have the following claim.

Claim 9.3.1.

2 | R_{0} | \geq | S^{*} ∖ S_{1}^{*} |

We create a set of pairs

P

as follows. There will be exactly

k

pairs. For each

f^{*} \in S_{1}^{*}

we add the pair

(r, f^{*})

where

ρ (f^{*}) = r

. For each

f^{*} \in S^{*} ∖ S_{1}^{*}

we add a pair

(r, f^{*})

where

r

is any arbitrary center from

R_{0}

- however we make sure that a center

r \in R_{0}

is in at most two pairs in

P

; this is possible because of Claim 9.3.1.

The pairs satisfy the following property.

Claim 9.3.2. If $(r, f^{*}) \in P$ then for any facility ${\hat{f}}^{*} \neq f^{*}, ρ ({\hat{f}}^{*}) \neq r$ .

The intuition for the pairs is as follows. If

ρ (f^{*}) = {r}

then we are essentially forced to consider the pair

(r, f^{*})

since

r

could be the only center near

f^{*}

with all other centers from

S

very far away. When considering the swap

(r, f^{*})

we can move the clients connecting to

r

f^{*}

. On the other hand if

| ρ^{- 1} (r) | > 1

then

r

is close to several centers in

S^{*}

and may be serving many clients. Thus we do not consider such an

r

in the swap pairs.

Figure 9.2: Mapping between

S^{*}

and

S

with each

f^{*} \in S^{*}

mapped to its closest center in

S

The main technical claim about the swap pairs is the following.

Lemma 9.10. Let $(r, f^{*})$ be a pair in $P$ . Then

0 \leq cost (S + f^{*} - r) - cost (S) \leq \sum_{j \in N^{*} (f^{*})} (O_{j} - A_{j}) + \sum_{j \in N (r)} 2 O_{j} .

We defer the proof of the lemma for now and use it to show that

cost (S) \leq

5 cost (S^{*})

. We sum over all pairs

(r, f^{*}) \in P

and note that each

f^{*} \in S

occurs in exactly one pair and each

r \in S

occurs in at most two pairs. Note that

O_{j} - A_{j}

can be a negative number while

O_{j}

is non-negative number.

\begin{aligned} 0 & \leq \sum_{(r, f^{*}) \in P} (\sum_{j \in N^{*} (f^{*})} (O_{j} - A_{j}) + \sum_{j \in N (r)} 2 O_{j}) \\ \leq \sum_{f^{*} \in S^{*}} \sum_{j \in N^{*} (f^{*})} (O_{j} - A_{j}) + 2 \sum_{r \in S} \sum_{j \in N (r)} 2 O_{j} \\ \leq cost (S^{*}) - cost (S) + 4 cost (S^{*}) . \end{aligned}

This implies the desired inequality.

Figure 9.3: Two cases in proof of Lemma 9.10. Consider the swap pair

(r, f^{*})

In the figure on the left the client

j \in N^{*} (f^{*})

is assigned to

f^{*}

. In the figure on the right the client

j \in N (r) ∖ N^{*} (f^{*})

is is assigned to

r^{'} = ρ ({\hat{f}}^{*})

Proof of Lemma 9.10. Since

S

is a local optimum swapping

r

with

f^{*}

cannot improve the cost and hence we obtain

cost (S + f^{*} - r) - cost (S) \geq 0

. We focus on the more interesting inequality. To bound the increase in cost by removing

r

and adding

f^{*}

S

, we consider a specific assignment of the clients. Any client

j \in N^{*} (f^{*})

is assigned to

f^{*}

(even if it is suboptimal to do so). See Figure 9.3. For such a client

j

the change in cost is

O_{j} - A_{j}

. Now consider any client

j \in N (r) ∖ N^{*} (f^{*})

; since

r

is no longer available we need to assign it to another facility. Which one? Let

ϕ^{*} (j) = {\hat{f}}^{*}

be the facility that

j

is assigned to in the optimum solution. Note that

{\hat{f}}^{*} \neq f^{*}

. We assign

j

r^{'} = ρ ({\hat{f}}^{*})

; from Claim 9.3.2,

ρ ({\hat{f}}^{*}) \neq r

and hence

r^{'} \in S - r + f^{*}

. See Figure 9.3. The change in the cost for such a client

j

d (j, r^{'}) - d (j, r)

. We bound it as follows

\begin{aligned} d (j, r^{'}) - d (j, r) & \leq d (j, {\hat{f}}^{*}) + d ({\hat{f}}^{*}, r^{'}) - d (j, r) (via triangle inequality) \\ \leq d (j, {\hat{f}}^{*}) + d ({\hat{f}}^{*}, r) - d (j, r) (since r^{'} is closest to {\hat{f}}^{*} in S) \\ \leq d (j, {\hat{f}}^{*}) + d (j, {\hat{f}}^{*}) (via triangle inequality) \\ = 2 O_{j} . \end{aligned}

Every other client is assigned to its existing center in

S

. Thus the total increase in the cost is obtained as

\sum_{j \in N^{*} (f^{*})} (O_{j} - A_{j}) + \sum_{j \in N (r) ∖ N^{*} (f^{*})} 2 O_{j} \leq \sum_{j \in N^{*} (f^{*})} (O_{j} - A_{j}) + \sum_{j \in N (r)} 2 O_{j} .

See ^[20:1] for a description of the tight example. The example in the conference version of the paper is buggy.

9.4 $k$ -Means

The

k

-Means problem is very popular in practice for a variety of reasons. In terms of center-based clustering the goal is to choose

k

centers

C = {c_{1}, c_{2}, \dots, c_{k}}

to minimize

\sum_{p} d (p, C)^{2}

. In the discrete setting one can obtain constant factor approximation algorithms via several techniques that follow the approach of

k

-Median. We note that the squared distance does not satisfy triangle inequality, however, it satisfies a relaxed triangle inequality and this is sufficient to generalize certain LP rounding and local search techniques.

In practice the continuous version is popular for clustering applications. The input points

P

are in the Euclidean space

R^{d}

where

d

is typically large. Let

X

be the set of input points where each

x \in X

is now a

d

-dimensional vector. The centers are now allowed to be in ambient space. This is called the Euclidean

k

-Means. Here the squared distance actually helps in a certain sense. For instance if

k = 1

then we can see that the optimum center is simply obtained by

\frac{1}{| P |} \sum_{x \in X} x

- in other words we take the "average". One can see this by considering the problem of finding the center as an optimization problem:

min_{y \in R^{d}} \sum_{x \in X} ∥ x - y ∥_{2}^{2} = min_{y \in R^{d}} \sum_{x \in X} {(x_{i} - y_{i})}^{2}

. It can be seen that we can optimize in each dimension separately and that the optimum in dimension

i

, can be seen to be

y_{i}^{*} = \frac{1}{| X |} x_{i}

. Surprisingly, hardness results for Euclidean

k

-Means were only established in the last decade. NP-hardness even when

d = 2

is established in ^[22], and APX-hardness for high dimensions was shown in ^[23].

9.5 Lloyd's algorithm, $D^{2}$ -sampling and $k$ -Means $+ +$

Llyod's algorithm is a very well-known and widely used heuristic for the Euclidean

k

-Means problem. It can viewed as a local search algorithm with an alternating optimization flavor. The algorithm starts with a set of centers

c_{1}, c_{2}, \dots, c_{k}

which are typically chosen randomly in some fashion. The centers define clusters

C_{1}, C_{2}, \dots, C_{k}

based on assigning each point to its nearest center. That is

C_{i}

is the set of all points in

X

that are closest to

c_{i}

(ties broken in some fashion). Once we have the clusters, via the observation above, one can find the best center for that cluster by taking the mean of the points in the cluster. That is, for each

C_{i}

we find a new center

c_{i}^{'} = \frac{1}{| C_{i} |} \sum_{x \in C_{i}} x

(if

C_{i}

is empty we simply set

c_{i}^{'} = c_{i}

). Thus we have a new set of centers and we repeat the process until convergence or some time limit. It is clear that the cost can only improve by recomputing the centers since we know the optimum center for

k = 1

is obtained by using the average of the points.

\underline{Lloyds- k -Means (X, k)}

1. Seeding: Pick

k

centers

c_{1}, c_{2}, \dots, c_{k}

2. repeat

A. Find clusters

C_{1}, C_{2}, \dots, C_{k}

where

C_{i} = {x \in X ∣ c_{i}

is closest center to

x}

\cos t = \sum_{i = 1}^{k} \sum_{x \in C_{i}} d {(x, c_{i})}^{2}

C. For

i = 1

k

c_{i} = \frac{1}{| C_{i} |} \sum_{x \in C_{i}} x

3. Until cost improvement is too small

4. Output clusters

C_{1}, C_{2}, \dots, C_{k}

There are two issues with the algorithm. The first issue is that the algorithm can, in the worst-case, run for an exponential number of iterations. This issue is common for many local search algorithms and as we discussed, it can be overcome by stopping when cost improvement is too small in a relative sense. The second issue is the more significant one. The algorithm can get stuck in a local optimum which can be arbitrarily bad when compared to the optimum solution. See figure below for a simple example.

Figure 9.4: Example demonstrating that a local optimum for Lloyd's algorithm can be arbitrarily bad compared to the optimum clustering. The green clusters are the optimum ones and the red ones are the local optimum.

$D^{2}$ -sampling and $k$ -Means++: To overcome the bad local optima it is common to run the algorithm with random starting centers. Arthur and Vassilvitskii ^[24] suggested a specific random sampling scheme to initialize the centers that is closely related to independent work in ^[25]. This is called

D^{2}

sampling.

\underline{D^{2} -sampling- k -Means + + (X, k)}

S = {c_{1}}

where

c_{1}

is chosen uniformly from

X

2. for

i = 2

k

A. Choose

c_{i}

randomly where

P [c_{i} = x] ≃ d (x, S)^{2}

S \leftarrow S \cup {c_{i}}

3. Output

S

k

-Means ++ is Lloyd's algorithm intialized with

k

centers obtained from

D^{2}

sampling.

Theorem 9.8 (^[24:1]). Let $S$ be the output of Lloyd's algorithm initialized with $D^{2}$ sampling. Then $E [cost (S)] \leq 8 (\ln k + 2) OPT$ . Moreover there are examples showing that it is no better than $2 \ln k$ competitive.

The analysis establishes that the seeding already creates a good approximation, so in a sense the local search is only refining the initial approximation. ^[26], ^[27] show that if one uses

O (k)

centers, initialized according to

D^{2}

sampling, then the local optimum will yield a constant factor approximation with constant probability; note that this is a bicriteria approximation where the number of centers is a constant factor more than

k

and the cost is being compared with respect to the optimum cost with

k

centers. The authors also show that there is a subset of

k

centers from the output of the algorithm that yields a constant factor approximation. One can then run a discrete optimization algorithm using the centers. Another interesting result based on

D^{2}

-sampling ideas yields a PTAS but the running time is of the form

O (n d 2^{\tilde{O} (k^{2} / ϵ)})

^[28]. See ^[29] for a scalable version of

k

-Means

+ +

Tomás Feder and Daniel Greene. “Optimal algorithms for approximate clustering”. In: Proceedings of the twentieth annual ACM symposium on Theory of computing. 1988, pp. 434–444. ↩︎
Dorit S Hochbaum and David B Shmoys. “A unified approach to approximation algorithms for bottleneck problems”. In: Journal of the ACM (JACM) 33.3 (1986), pp. 533–550. ↩︎
Viswanath Nagarajan, Baruch Schieber, and Hadas Shachnai. “The Euclidean k-supplier problem”. In: Mathematics of Operations Research 45.1 (2020), pp. 1–14. ↩︎
Deeparnab Chakrabarty and Maryam Negahbani. “Generalized center problems with outliers”. In: ACM Transactions on Algorithms (TALG) 15.3 (2019), pp. 1–14. ↩︎
Ján Plesnık. “A heuristic for the p-center problems in graphs”. In: Discrete Applied Mathematics 17.3 (1987), pp. 263–268. ↩︎
Tanvi Bajpai, Deeparnab Chakrabarty, Chandra Chekuri, and Maryam Negahbani. “Revisiting Priority $k$ -Center: Fairness and Outliers”. In: arXiv preprint arXiv:2103.03337 (2021). ↩︎
Sariel Har-Peled and Manor Mendel. “Fast construction of nets in lowdimensional metrics and their applications”. In: SIAM Journal on Computing 35.5 (2006), pp. 1148–1184. ↩︎
Sariel Har-Peled and Benjamin Raichel. “Net and prune: A linear time algorithm for euclidean distance problems”. In: Journal of the ACM (JACM) 62.6 (2015), pp. 1–35. ↩︎
Shi Li. “A 1.488 approximation algorithm for the uncapacitated facility location problem”. In: Information and Computation 222 (2013), pp. 45–58. ↩︎
Sudipto Guha and Samir Khuller. “Greedy strikes back: Improved facility location algorithms”. In: Journal of algorithms 31.1 (1999), pp. 228–248. ↩︎
Kamal Jain and Vĳay V Vazirani. “Approximation algorithms for metric facility location and k-median problems using the primal-dual schema and Lagrangian relaxation”. In: Journal of the ACM (JACM) 48.2 (2001), pp. 274–296. ↩︎
Reuven Bar-Yehuda and Shimon Even. “A linear-time approximation algorithm for the weighted vertex cover problem”. In: Journal of Algorithms 2.2 (1981), pp. 198–203. ↩︎
Vĳay V Vazirani. Approximation algorithms. Springer Science & Business Media, 2013. ↩︎ ↩︎
Yair Bartal, Moses Charikar, and Danny Raz. “Approximating Min-Sum k-Clustering in Metric Spaces”. In: Proceedings of the Thirty-Third Annual ACM Symposium on Theory of Computing. STOC ’01. Hersonissos, Greece: Association for Computing Machinery, 2001, pp. 11–20. isbn: 1581133499. doi: 10.1145/380752.380754. url: https://doi.org/10.1145/380752.380754. ↩︎
Madhukar R Korupolu, C Greg Plaxton, and Rajmohan Rajaraman. “Analysis of a local search heuristic for facility location problems”. In: Journal of algorithms 37.1 (2000), pp. 146–188. ↩︎
David P Williamson and David B Shmoys. The design of approximation algorithms. Cambridge university press, 2011. ↩︎ ↩︎
Moses Charikar, Sudipto Guha, Éva Tardos, and David B Shmoys. “A constant-factor approximation algorithm for the k-median problem”. In: Journal of Computer and System Sciences 65.1 (2002), pp. 129–149. ↩︎
Ravishankar Krishnaswamy, Amit Kumar, Viswanath Nagarajan, Yogish Sabharwal, and Barna Saha. “The matroid median problem”. In: Proceedings of the twenty-second annual ACM-SIAM symposium on Discrete Algorithms. SIAM. 2011, pp. 1117–1130. ↩︎
Chaitanya Swamy. “Improved approximation algorithms for matroid and knapsack median problems and applications”. In: ACM Transactions on Algorithms (TALG) 12.4 (2016), pp. 1–22. ↩︎
Vĳay Arya, Naveen Garg, Rohit Khandekar, Adam Meyerson, Kamesh Munagala, and Vinayaka Pandit. “Local search heuristics for k-median and facility location problems”. In: SIAM Journal on computing 33.3 (2004), pp. 544–562. ↩︎ ↩︎
Anupam Gupta and Kanat Tangwongsan. Simpler Analyses of Local Search Algorithms for Facility Location. 2008. arXiv: 0809.2554 [cs.DS]. ↩︎
Meena Mahajan, Prajakta Nimbhorkar, and Kasturi Varadarajan. “The planar k-means problem is NP-hard”. In: Theoretical Computer Science 442 (2012), pp. 13–21. ↩︎
Pranjal Awasthi, Moses Charikar, Ravishankar Krishnaswamy, and Ali Kemal Sinop. “The Hardness of Approximation of Euclidean k-Means”. In: 31st International Symposium on Computational Geometry (SoCG 2015). Ed. by Lars Arge and János Pach. Vol. 34. Leibniz International Proceedings in Informatics (LIPIcs). Dagstuhl, Germany: Schloss Dagstuhl–LeibnizZentrum fuer Informatik, 2015, pp. 754–767. isbn: 978-3-939897-83-5. doi: 10.4230/LIPIcs.SOCG.2015.754. url: http://drops.dagstuhl.de/opus/ volltexte/2015/5117. ↩︎
Sergei Vassilvitskii and David Arthur. “k-means++: The advantages of careful seeding”. In: Proceedings of the eighteenth annual ACM-SIAM symposium on Discrete algorithms. 2006, pp. 1027–1035. ↩︎ ↩︎
Rafail Ostrovsky, Yuval Rabani, Leonard J Schulman, and Chaitanya Swamy. “The effectiveness of Lloyd-type methods for the k-means problem”. In: Journal of the ACM (JACM) 59.6 (2013), pp. 1–22. ↩︎
Ankit Aggarwal, Amit Deshpande, and Ravi Kannan. “Adaptive sampling for k-means clustering”. In: Approximation, Randomization, and Combinatorial Optimization. Algorithms and Techniques. Springer, 2009, pp. 15–28. ↩︎
Nir Ailon, Ragesh Jaiswal, and Claire Monteleoni. “Streaming k-means approximation.” In: NIPS. Vol. 22. 2009, pp. 10–18. ↩︎
Ragesh Jaiswal, Amit Kumar, and Sandeep Sen. “A simple D 2-samplingbased PTAS for k-means and other clustering problems”. In: Algorithmica 70.1 (2014), pp. 22–46. ↩︎
Bahman Bahmani, Benjamin Moseley, Andrea Vattani, Ravi Kumar, and Sergei Vassilvitskii. “Scalable K-Means++”. In: Proc. VLDB Endow. 5.7 (Mar. 2012), pp. 622–633. issn: 2150-8097. doi: 10.14778/2180912.2180915. url: https://doi.org/10.14778/2180912.2180915. ↩︎

Approximation Algorithms: Clustering and Facility Location

9.1 k k kkk-Center

9.1.1 Gonzalez's algorithm and nets in metric spaces

9.1.2 Hochbaum-Shmoys bottleneck approach

9.1.3 Related Problems and Discussion

9.2 Uncapacitated Facility Location

9.2.1 LP Rounding

9.2.2 Primal-Dual

9.2.3 Local Search

9.3 k k kkk-Median

9.3.1 Local Search

9.4 k k kkk-Means

9.5 Lloyd's algorithm, D 2 D 2 D^(2)D^{2}D2-sampling and k k kkk-Means + + + + ++++++

Recommended for you

Report Article

9.1 $k$ -Center

9.3 $k$ -Median

9.4 $k$ -Means

9.5 Lloyd's algorithm, $D^{2}$ -sampling and $k$ -Means $+ +$