<p>We improve Chandra's bound, as he conjectured was possible, giving an approximation algorithm that opens <span class="math-container">$f(k,\epsilon)=O(k\log (1/\epsilon))$</span> facilities to obtain assignment cost at most <span class="math-container">$1+\epsilon$</span> times the optimal with <span class="math-container">$k$</span> facilities.</p>
<p><strong>Theorem 1.</strong> <em>There is a polynomial-time algorithm that, given any metric <span class="math-container">$k$</span>-medians instance, returns a solution that uses <span class="math-container">$O(k\log(1/\epsilon))$</span> facilities with assignment cost <span class="math-container">$1+O(\epsilon)$</span> times the optimum with <span class="math-container">$k$</span> facilities.</em></p>
<p><em>Proof.</em> We use the same idea that Chandra does, but instead of using Lin and Vitter's filtering technique, we adapt the improved rounding scheme from <a href="https://scholar.google.com/scholar?cluster=502252032625720299&amp;hl=en&amp;as_sdt=0,22" rel="nofollow noreferrer">K-medians, facility location, and the Chernoff-Wald bound</a>.</p>
<p>Fix any optimal solution <span class="math-container">$(x, y)$</span> to the standard <span class="math-container">$k$</span>-median linear-program relaxation:
<span class="math-container">$$
\begin{align*}
&amp; &amp; \min \textstyle \sum_{j\in C}\sum_{f\in L} d(j, f) x_{jf} \\
&amp; (\forall j\in C) &amp; \textstyle \sum_{f\in L} x_{jf}  = 1, \\
&amp;  &amp; \textstyle \sum_{f\in L} y_f = k, \\
&amp; (\forall j\in C, f\in L) &amp; 0 \le x_{jf} \le y_f.
\end{align*}
$$</span>
Let <span class="math-container">$\alpha_j = \sum_{f\in L} d(j, f) x_{jf}$</span> be the assignment cost incurred by <span class="math-container">$(x, y)$</span> for client <span class="math-container">$j$</span>.
Let OPT <span class="math-container">$= \sum_j \alpha_j$</span> be the total assignment cost.
The LP is a relaxation of the actual instance, so the true optimum is at least OPT.</p>
<p>The rounding scheme has two phases.
In the first phase it generates a set <span class="math-container">$S\subseteq L$</span> of facilities and a partial assignment of clients to facilities in <span class="math-container">$S$</span>, by sampling <span class="math-container">$h=\lceil k\ln(1/\epsilon)\rceil$</span> times as follows:</p>
<ol>
<li>for <span class="math-container">$t=1,2,\ldots, h$</span>:</li>
<li><span class="math-container">$~~~~$</span> sample a random facility <span class="math-container">$f_t\in L$</span> from the distribution defined by <span class="math-container">$y/k$</span></li>
<li><span class="math-container">$~~~~$</span> for each unassigned client <span class="math-container">$j\in C$</span>, with probability <span class="math-container">$x_{jf_t}/y_{f_t}$</span>, assign <span class="math-container">$j$</span> to facility <span class="math-container">$f_t$</span></li>
<li>let <span class="math-container">$A$</span> be the set of assigned clients</li>
</ol>
<p>In each iteration, the chance of assigning a given unassigned client <span class="math-container">$j$</span> to a given facility <span class="math-container">$f$</span> is <span class="math-container">$\frac{y_f}{k} \frac{x_{jf}}{y_f} = \frac{x_{jf}}{k}$</span>.</p>
<p>In each iteration, the probability of assigning a given unassigned client <span class="math-container">$j$</span> to some facility is <span class="math-container">$\sum_{f} \frac{x_{jf}}{k} = \frac 1 k$</span>.  So the chance that any given client is not in <span class="math-container">$A$</span> (not assigned during the first phase) is <span class="math-container">$p = (1-1/k)^h \le \exp(-h/k) \le \epsilon$</span>.</p>
<p>For any given client <span class="math-container">$j$</span>, conditioned on <span class="math-container">$j$</span> being in <span class="math-container">$A$</span> (assigned to some facility during the first phase), the expected cost of its assignment is <span class="math-container">$\sum_f x_{jf} d(j, f) = \alpha_j$</span>.
Thus the expected cost of the first phase is <span class="math-container">$\sum_j (1-p)\alpha_j = (1-p)$</span>OPT.</p>
<p>In the second phase, (following Chandra's approach) use any existing <span class="math-container">$O(1)$</span>-approximation algorithm to find an <span class="math-container">$O(1)$</span>-approximate solution (w.r.t. the induced LP, and using here that the problem is metric) to the <span class="math-container">$k$</span>-medians instance obtained from the given instance by removing all clients that were assigned in the first phase.  This step opens at most <span class="math-container">$k$</span> facilities and (because the induced LP's solution costs at most <span class="math-container">$\sum_{j\not \in A} \alpha_j$</span>)
incurs assignment cost <span class="math-container">$O(\sum_{j\not\in A} \alpha_j)$</span>.
So by linearity of expectation the expected cost of this phase is
<span class="math-container">$O(\sum_j \Pr[j\not\in A] \alpha_j) = O(\sum_j p \alpha_j) = O(p \,\text{OPT})$</span>.</p>
<p>Thus, in total at most <span class="math-container">$k + 1 + k\ln(1/\epsilon)$</span> facilities are opened, and the expected assignment cost is at most
<span class="math-container">$(1-p)\text{OPT} + O(p\,\text{OPT}) = (1+O(\epsilon))$</span>OPT.
<span class="math-container">$~~~\Box$</span></p>

CC BY-SA 4.0
2021-06-15T16:46:48.733Z
Nice. I was going to try the randomized algorithm/analysis but you were too fast :). One should also check and compare with Hsu-Telgarsky analysis who basically do greedy. I will look later today.
Chandra Chekuri
CC BY-SA 4.0
2021-06-15T17:07:34.870Z
The paper of Makarychev, Makarychev, Ward essentially has this implicitly. Their main goal is k-means in high-dimensional Euclidean space and in fact they have to deal with the issue of centers being feasible in the ambient space. https://arxiv.org/abs/1507.04227.
Chandra Chekuri
CC BY-SA 4.0
2021-06-15T17:52:55.000Z
@ChandraChekuri Just want to add something more. The centers being feasible in the ambient space do not cause a problem in the Euclidean space, since an $\alpha$-approximation for the discrete k-median gives an $\alpha(1+\epsilon)$-approximation in the continuous Euclidean space. The reduction follows from the weak-coreset construction of [Feldman et al.](http://people.csail.mit.edu/dannyf/kmeancoreset.pdf) :-)
Inuyasha Yagami
CC BY-SA 4.0
2021-06-15T21:06:00.457Z
This may generalize to the problem with outliers since it basically relies on having a constant factor approximation with respect to LP for an underlying problem. Perhaps there is a simple framework here that may work for a class of problems.
Chandra Chekuri
CC BY-SA 4.0
2021-06-15T21:27:45.850Z
The Hsu-Telgarsky paper also works I believe. They start with a constant factor solution and apply Greedy to reduce the cost and they essentially only use sub/supermodularity of the cost which is the essence. Anyway, now I believe that both Makarychev, Makarychev, Ward paper and Hsu-Telgarsky paper have the desired result with the latter being more general in some sense.
Chandra Chekuri
CC BY-SA 4.0
2021-06-16T12:03:06.277Z
Thanks Neal. I have a doubt in computation of expected cost in phase 1. Is the expected cost of client $j$ is derived in this way: $= \frac{\alpha_j}{k} + \left( 1-\frac{1}{k} \right) \cdot \frac{\alpha_j}{k} + \dotsc + \left( 1-\frac{1}{k} \right)^{h-1} \cdot \frac{\alpha_j}{k} =\left( 1-\left( 1-\frac{1}{k} \right) ^{h} \right) \cdot \alpha_j = (1-p) \cdot \alpha_j$. or there is some shorter way that I am missing? Thanks!
Inuyasha Yagami
CC BY-SA 4.0
2021-06-16T13:15:12.887Z
Given $j$, you can also just observe that, for any given $f'$ and iteration $t$, the probability that $j$ is assigned to $f'$ in iteration $t$, given that $j$ is assigned to some facility in iteration $t$, is $\frac{(x_{jf'}/k)}{1/k} = x_{jf'}$.  So, given that $j$ is assigned in some iteration, the probability that it is assigned to $f'$ is $x_{jf'}$.  Then the claim follows by linearity of expectation.
Neal Young
CC BY-SA 4.0
2021-06-16T14:00:38.977Z
Thanks! Just one last thing. I think there is typo in the main theorem and the first paragraph. It should be $O( k \log(1/\epsilon))$.
Inuyasha Yagami
CC BY-SA 4.0
2021-09-30T22:25:29.147Z
Thanks for the bounty!
Neal Young
CC BY-SA 4.0
2021-06-15T15:12:42.993Z
2021-06-16T19:48:34.243Z
2021-06-16T19:48:34.243Z
Neal Young
<p>Historically Lin and Vitter showed that one can obtain a <span class="math-container">$(2+\epsilon, k(1+1/\epsilon))$</span> bi-criteria approximation via a simple filtering trick with respect to the LP. This was before a constant factor approximation was shown via the LP. After that there seemed to be less interest in a bi-criteria PTAS even though it was explored extensively in geometric settings. Here is a sketch for <span class="math-container">$k$</span>-median that seems to give a <span class="math-container">$(1+\epsilon)$</span>-approximation using <span class="math-container">$O(k\log(1/\epsilon)/\epsilon)$</span> centers.</p>
<p>Consider solving the standard LP for <span class="math-container">$k$</span>-median. Let <span class="math-container">$y_a$</span> be the amount of center opened at location <span class="math-container">$a$</span>. For each point <span class="math-container">$v$</span> we can define <span class="math-container">$\alpha_v$</span> to be the distance paid by <span class="math-container">$v$</span> in the LP. Note that the total LP cost is <span class="math-container">$\sum_v \alpha_v$</span> which is a lower bound on the optimum integer cost. For each <span class="math-container">$v$</span> define <span class="math-container">$B_v$</span> to be a ball of radius <span class="math-container">$(1+\epsilon)\alpha_v$</span>. Via Lin-Vitter style filtering idea, we can obtain a new LP solution that uses <span class="math-container">$O(k/\epsilon)$</span> centers with the property that for each <span class="math-container">$v$</span>, there is a total fractional center of at least <span class="math-container">$1$</span> in <span class="math-container">$B_v$</span>. We now set up a weighted Set Cover instance. We want to find centers such that every ball <span class="math-container">$B_v$</span> is hit by a center. We have a fractional solution with <span class="math-container">$O(k/\epsilon)$</span> centers. We can show that Greedy with <span class="math-container">$O(k \log (1/\delta)/\epsilon)$</span> centers will cover at least <span class="math-container">$(1-\delta)$</span> fraction of the balls. We can do this in the weighted sense as well where weight of <span class="math-container">$v$</span> is <span class="math-container">$\alpha_v$</span>. So we run Greedy with <span class="math-container">$\delta = \epsilon$</span> and obtain a collection of <span class="math-container">$O(k \log (1/\epsilon)/\epsilon)$</span> centers. Let <span class="math-container">$H$</span> be the set of points <span class="math-container">$v$</span> such that <span class="math-container">$B_v$</span> is not hit by the chosen centers. We have <span class="math-container">$\sum_{v\in H} \alpha_v \le  \epsilon \sum_v \alpha_v$</span>. For these uncovered points we can use a constant factor (or even bicriteria) approximation for <span class="math-container">$k$</span>-median via the LP. Thus the total cost will be <span class="math-container">$(1+O(\epsilon))\sum_v \alpha_v$</span> and the total number of centers used will be <span class="math-container">$O(k \log(1/\epsilon)/\epsilon)$</span>.</p>
<p>I wonder if one can use only <span class="math-container">$O(k\log(1/\epsilon))$</span> centers. I think it may be possible by being more careful.</p>

CC BY-SA 4.0
2021-06-15T13:38:26.353Z
Chandra, very nice.  BTW which result by Lin and Vitter do you have in mind? (I know about their paper "$\epsilon$-Approximations with Packing Constraint Violation", but that is for non-metric and gives a $(1+\epsilon, k(1+1/\epsilon)(1+\ln n))$-bicriteria approximation [Theorem 1]), which was later improved to $(1+\epsilon, O(k\log(n/\epsilon)))$ and then to $(1, O(k \log n))$ [here](https://dx.doi.org/10.1007/s00453-007-9005-x), Theorem 5.)  I'll think about whether we can reduce the number of centers for metric, using those and your ideas.
Neal Young
CC BY-SA 4.0
2021-06-15T13:00:37.117Z
2021-06-15T13:00:37.117Z
Chandra Chekuri
<p>The <span class="math-container">$k$</span>-median problem is defined as follows: Given a set <span class="math-container">$C$</span> of clients and a set <span class="math-container">$L$</span> of facility locations defined over a distance metric <span class="math-container">$d$</span>. Open a set <span class="math-container">$F$</span> of <span class="math-container">$k$</span> facility in <span class="math-container">$L$</span> such that the total cost of assigning the clients to their closest facilities is minimized, i.e.,</p>
<p><span class="math-container">$$minimize \quad \Phi(F,C) \equiv \sum_{j = 1}^{|C|} \, \min_{f \in F} \big\{ d(j,f) \big\}$$</span></p>
<p>Suppose the optimal cost of the instance is denoted by <span class="math-container">$OPT$</span>.</p>
<p>Is there any bi-criteria approximation algorithm for the problem that opens <span class="math-container">$f(k,\epsilon)$</span> facilities and gives <span class="math-container">$(1+\epsilon)$</span>-approximation with respect to <span class="math-container">$OPT$</span> for any computable function <span class="math-container">$f$</span>? The running time of the algorithm is allowed to be <span class="math-container">$g(k,\epsilon) \cdot n^{O(1)}$</span> for any computable function <span class="math-container">$g$</span>.</p>
<p>Note: I know that in the Euclidean space the problem does have such a type of algorithm. But I am particularly interested in general metric spaces.</p>

CC BY-SA 4.0
2021-06-15T02:15:40.400Z
Are you asking for a polynomial-time algorithm?
Neal Young
CC BY-SA 4.0
2021-06-15T03:14:43.197Z
@NealYoung That was an important point. I have edited the question. Thank You!
Inuyasha Yagami
CC BY-SA 4.0
2021-06-15T12:37:14.620Z
I though such a result should be known but I couldn't find it. A paper by Hsu and Telgarsky may be helpful to look at. https://arxiv.org/abs/1607.06203
Chandra Chekuri
<approximation-algorithms><clustering>
Is there any Bi-criteria PTAS for Metric $k$-Median?
Inuyasha Yagami
# Is there any Bi-criteria PTAS for Metric $k$-Median?

(Asked by: username_0 on 2021-06-14)

The $k$-median problem is defined as follows: Given a set $C$ of clients and a set $L$ of facility locations defined over a distance metric $d$. Open a set $F$ of $k$ facility in $L$ such that the total cost of assigning the clients to their closest facilities is minimized, i.e.,

$$minimize \quad \Phi(F,C) \equiv \sum\_{j = 1}^{|C|} \, \min\_{f \in F} \big\{ d(j,f) \big\}$$

Suppose the optimal cost of the instance is denoted by $OPT$.

Is there any bi-criteria approximation algorithm for the problem that opens $f(k,\epsilon)$ facilities and gives $(1+\epsilon)$-approximation with respect to $OPT$ for any computable function $f$? The running time of the algorithm is allowed to be $g(k,\epsilon) \cdot n^{O(1)}$ for any computable function $g$.

Note: I know that in the Euclidean space the problem does have such a type of algorithm. But I am particularly interested in general metric spaces.

**Tags:** approximation-algorithms, clustering

> **Comments:**
> **username_1**: Are you asking for a polynomial-time algorithm?  
> **username_0**: @NealYoung That was an important point. I have edited the question. Thank You!  
> **username_2**: I though such a result should be known but I couldn't find it. A paper by Hsu and Telgarsky may be helpful to look at. https://arxiv.org/abs/1607.06203  

---

## 2 answers

### Answer by username_1 on 2021-06-15. Score: 4 | Accepted answer

We improve Chandra's bound, as he conjectured was possible, giving an approximation algorithm that opens $f(k,\epsilon)=O(k\log (1/\epsilon))$ facilities to obtain assignment cost at most $1+\epsilon$ times the optimal with $k$ facilities.

**Theorem 1.** *There is a polynomial-time algorithm that, given any metric $k$-medians instance, returns a solution that uses $O(k\log(1/\epsilon))$ facilities with assignment cost $1+O(\epsilon)$ times the optimum with $k$ facilities.*

*Proof.* We use the same idea that Chandra does, but instead of using Lin and Vitter's filtering technique, we adapt the improved rounding scheme from [K-medians, facility location, and the Chernoff-Wald bound](https://scholar.google.com/scholar?cluster=502252032625720299&hl=en&as_sdt=0,22).

Fix any optimal solution $(x, y)$ to the standard $k$-median linear-program relaxation:
$$
\begin{align\*}
& & \min \textstyle \sum\_{j\in C}\sum\_{f\in L} d(j, f) x\_{jf} \\
& (\forall j\in C) & \textstyle \sum\_{f\in L} x\_{jf} = 1, \\
& & \textstyle \sum\_{f\in L} y\_f = k, \\
& (\forall j\in C, f\in L) & 0 \le x\_{jf} \le y\_f.
\end{align\*}
$$
Let $\alpha\_j = \sum\_{f\in L} d(j, f) x\_{jf}$ be the assignment cost incurred by $(x, y)$ for client $j$.
Let OPT $= \sum\_j \alpha\_j$ be the total assignment cost.
The LP is a relaxation of the actual instance, so the true optimum is at least OPT.

The rounding scheme has two phases.
In the first phase it generates a set $S\subseteq L$ of facilities and a partial assignment of clients to facilities in $S$, by sampling $h=\lceil k\ln(1/\epsilon)\rceil$ times as follows:

1. for $t=1,2,\ldots, h$:
2. $~~~~$ sample a random facility $f\_t\in L$ from the distribution defined by $y/k$
3. $~~~~$ for each unassigned client $j\in C$, with probability $x\_{jf\_t}/y\_{f\_t}$, assign $j$ to facility $f\_t$
4. let $A$ be the set of assigned clients

In each iteration, the chance of assigning a given unassigned client $j$ to a given facility $f$ is $\frac{y\_f}{k} \frac{x\_{jf}}{y\_f} = \frac{x\_{jf}}{k}$.

In each iteration, the probability of assigning a given unassigned client $j$ to some facility is $\sum\_{f} \frac{x\_{jf}}{k} = \frac 1 k$. So the chance that any given client is not in $A$ (not assigned during the first phase) is $p = (1-1/k)^h \le \exp(-h/k) \le \epsilon$.

For any given client $j$, conditioned on $j$ being in $A$ (assigned to some facility during the first phase), the expected cost of its assignment is $\sum\_f x\_{jf} d(j, f) = \alpha\_j$.
Thus the expected cost of the first phase is $\sum\_j (1-p)\alpha\_j = (1-p)$OPT.

In the second phase, (following Chandra's approach) use any existing $O(1)$-approximation algorithm to find an $O(1)$-approximate solution (w.r.t. the induced LP, and using here that the problem is metric) to the $k$-medians instance obtained from the given instance by removing all clients that were assigned in the first phase. This step opens at most $k$ facilities and (because the induced LP's solution costs at most $\sum\_{j\not \in A} \alpha\_j$)
incurs assignment cost $O(\sum\_{j\not\in A} \alpha\_j)$.
So by linearity of expectation the expected cost of this phase is
$O(\sum\_j \Pr[j\not\in A] \alpha\_j) = O(\sum\_j p \alpha\_j) = O(p \,\text{OPT})$.

Thus, in total at most $k + 1 + k\ln(1/\epsilon)$ facilities are opened, and the expected assignment cost is at most
$(1-p)\text{OPT} + O(p\,\text{OPT}) = (1+O(\epsilon))$OPT.
$~~~\Box$

> **Comments:**
> **username_2**: Nice. I was going to try the randomized algorithm/analysis but you were too fast :). One should also check and compare with Hsu-Telgarsky analysis who basically do greedy. I will look later today.  
> **username_2**: The paper of Makarychev, Makarychev, Ward essentially has this implicitly. Their main goal is k-means in high-dimensional Euclidean space and in fact they have to deal with the issue of centers being feasible in the ambient space. https://arxiv.org/abs/1507.04227.  
> **username_0**: @ChandraChekuri Just want to add something more. The centers being feasible in the ambient space do not cause a problem in the Euclidean space, since an $\alpha$-approximation for the discrete k-median gives an $\alpha(1+\epsilon)$-approximation in the continuous Euclidean space. The reduction follows from the weak-coreset construction of [Feldman et al.](http://people.csail.mit.edu/dannyf/kmeancoreset.pdf) :-)  
> **username_2**: This may generalize to the problem with outliers since it basically relies on having a constant factor approximation with respect to LP for an underlying problem. Perhaps there is a simple framework here that may work for a class of problems.  
> **username_2**: The Hsu-Telgarsky paper also works I believe. They start with a constant factor solution and apply Greedy to reduce the cost and they essentially only use sub/supermodularity of the cost which is the essence. Anyway, now I believe that both Makarychev, Makarychev, Ward paper and Hsu-Telgarsky paper have the desired result with the latter being more general in some sense.  
> **username_0**: Thanks Neal. I have a doubt in computation of expected cost in phase 1. Is the expected cost of client $j$ is derived in this way: $= \frac{\alpha_j}{k} + \left( 1-\frac{1}{k} \right) \cdot \frac{\alpha_j}{k} + \dotsc + \left( 1-\frac{1}{k} \right)^{h-1} \cdot \frac{\alpha_j}{k} =\left( 1-\left( 1-\frac{1}{k} \right) ^{h} \right) \cdot \alpha_j = (1-p) \cdot \alpha_j$. or there is some shorter way that I am missing? Thanks!  
> **username_1**: Given $j$, you can also just observe that, for any given $f'$ and iteration $t$, the probability that $j$ is assigned to $f'$ in iteration $t$, given that $j$ is assigned to some facility in iteration $t$, is $\frac{(x_{jf'}/k)}{1/k} = x_{jf'}$.  So, given that $j$ is assigned in some iteration, the probability that it is assigned to $f'$ is $x_{jf'}$.  Then the claim follows by linearity of expectation.  
> **username_0**: Thanks! Just one last thing. I think there is typo in the main theorem and the first paragraph. It should be $O( k \log(1/\epsilon))$.  
> **username_1**: Thanks for the bounty!  

---

### Answer by username_2 on 2021-06-15. Score: 1

Historically Lin and Vitter showed that one can obtain a $(2+\epsilon, k(1+1/\epsilon))$ bi-criteria approximation via a simple filtering trick with respect to the LP. This was before a constant factor approximation was shown via the LP. After that there seemed to be less interest in a bi-criteria PTAS even though it was explored extensively in geometric settings. Here is a sketch for $k$-median that seems to give a $(1+\epsilon)$-approximation using $O(k\log(1/\epsilon)/\epsilon)$ centers.

Consider solving the standard LP for $k$-median. Let $y\_a$ be the amount of center opened at location $a$. For each point $v$ we can define $\alpha\_v$ to be the distance paid by $v$ in the LP. Note that the total LP cost is $\sum\_v \alpha\_v$ which is a lower bound on the optimum integer cost. For each $v$ define $B\_v$ to be a ball of radius $(1+\epsilon)\alpha\_v$. Via Lin-Vitter style filtering idea, we can obtain a new LP solution that uses $O(k/\epsilon)$ centers with the property that for each $v$, there is a total fractional center of at least $1$ in $B\_v$. We now set up a weighted Set Cover instance. We want to find centers such that every ball $B\_v$ is hit by a center. We have a fractional solution with $O(k/\epsilon)$ centers. We can show that Greedy with $O(k \log (1/\delta)/\epsilon)$ centers will cover at least $(1-\delta)$ fraction of the balls. We can do this in the weighted sense as well where weight of $v$ is $\alpha\_v$. So we run Greedy with $\delta = \epsilon$ and obtain a collection of $O(k \log (1/\epsilon)/\epsilon)$ centers. Let $H$ be the set of points $v$ such that $B\_v$ is not hit by the chosen centers. We have $\sum\_{v\in H} \alpha\_v \le \epsilon \sum\_v \alpha\_v$. For these uncovered points we can use a constant factor (or even bicriteria) approximation for $k$-median via the LP. Thus the total cost will be $(1+O(\epsilon))\sum\_v \alpha\_v$ and the total number of centers used will be $O(k \log(1/\epsilon)/\epsilon)$.

I wonder if one can use only $O(k\log(1/\epsilon))$ centers. I think it may be possible by being more careful.

> **Comments:**
> **username_1**: Chandra, very nice.  BTW which result by Lin and Vitter do you have in mind? (I know about their paper "$\epsilon$-Approximations with Packing Constraint Violation", but that is for non-metric and gives a $(1+\epsilon, k(1+1/\epsilon)(1+\ln n))$-bicriteria approximation [Theorem 1]), which was later improved to $(1+\epsilon, O(k\log(n/\epsilon)))$ and then to $(1, O(k \log n))$ [here](https://dx.doi.org/10.1007/s00453-007-9005-x), Theorem 5.)  I'll think about whether we can reduce the number of centers for metric, using those and your ideas.  

---

