Faster and Simpler Greedy Algorithm for $k$-Median and $k$-Means

We see our main contribution as a deeper understanding of the greedy algorithm from Mettu and Plaxton. This results in some novel tradeoff between approximation and running time for $(k,z)$-clustering, that we compare below with other existing techniques.

We note first that our $n^{1+1/c}$ algorithms (in graphs and Euclidean space) results do not beat the state-of-the-art, in terms of approximation and running time: indeed, in sparse graph, an algorithm from [42] (which uses the primal-dual techniques as subroutine) runs in near-linear time and computes a constant-factor approximation to $k$-median. An equivalent result is not shown for $k$-means, or more generally $(k,z)$-clustering. While it is likely that the analysis would follow, using recent results on the primal-dual method for $k$-means, our result gives the first formal proof for a fast constant-factor approximation to $k$-means in graphs.

As any set of points in Euclidean space can be turned into a sparse graph in time $O\left(n^{1+1/c+o(1)}d\right)$ while preserving the pairwise distance up to a factor $O(c)$ (using spanners [31]), a constant-factor approximation for Euclidean inputs running in time $O\left(n^{1+1/c+o(1)}d\right)$ follows from Thorup’s result. Again, this formally holds for $k$-median, but not $k$-means.

However, in the own words of Thorup, this algorithm “is rather complicated and hence unlikely to be of direct practical relevance.” Our contribution here is to present a different algorithm, greedy and much simpler than the one of [42]. In addition, our algorithms works for the incremental version of the problems as well.

In the realm of linear time, we are not aware of any sub-polynomial approximation algorithm for Euclidean $(k,z)$-clustering: only the algorithm from [11] runs in near-linear time but yields a $O(k^4)$-approximation. The paper [26] claims that an algorithm from [22] would produce a $O(polylogn)$ approximation in near-linear time: however, their proof is not correct. More precisely, [22] present a near linear-time algorithm with no approximation guarantee: to achieve a $polylog(n)$-approximation ratio, they need a rejection-sampling step that adds a $n^{1+o(1)}$ running-time – as opposed to our near linear $npolylog(n)$. The other algorithm indeed runs in near-linear time, but does not have any approximation guarantee: indeed, their multi-tree embedding is shown to preserve squared distances, in expectation, up to a polylog factor. However, this embedding does not yield a metric space (in particular, it does not respect triangle inequality²²2The multi-tree embedding consists of the following: build several quadtrees, and define the distance between two points to be the minimum distance accross all quadtrees. This minimum may not respect triangle inequality. As an illustration, consider three points on the line, $a$ at $x=0$, $b$ at $x=2$ and $c$ at $x=4$. Consider two quadtrees, one with a split at $x=1$ (making all points with very far from all points with $x>1$, say distance $10$, while points with or $x>1$ are at their exact distance), the other one with a split at $x=3$. In both trees, $a$ and $c$ are at distance $10$, hence they are at distance 10 in the multi-tree. However, in the first one $b$ and $c$ are at distance $2$, and in the second $a$ and $b$ are at distance $2$: hence, with triangle inequality $a$ and $c$ should be at distance $4$.): we are not aware of any approximation algorithm in this setting – in particular, the $k$-means++ algorithm employed by [22] does not.

Our near linear-time approximation algorithm is therefore new, and, importantly, corrects the claim made in [26].

Jiang, Jin, Lou and Lu [32] took inspiration from an initial pre-print of our paper to show how to implement fast local search, in spaces that admit sparse spanners: they get a $O(c)$-approximation for $(k,z)$-clustering in time $n^{1+1/c+o(1)}$. Hence, they improve our tradeoff between run-time and approximation ratio, and achieve the same result as combining spanners and Thorup’s algorithm. Compared to the latter, their local search procedure has the advantage being much simpler, and their proof directly generalizes to $k$-means. Our initial pre-print only mentioned results in Euclidean space, while [32] applies to any metric with sparse spanners. We concurrently generalized our techniques to graphs – hence, to any metric that admits efficiently computable sparse spanners.

We note that this paper explicitly cites our work as a motivation: we present a $poly(c)$-approximation in time $n^{1+1/c+o(1)}$, and their goal was to bring the approximation-ratio to $O(c)$, closer to the “LSH tradeoff”. In addition, we note that their algorithm fails to get any near-linear time run-time, as they deeply rely on LSH.

Another preprint appeared on arxiv [25] after the initial release of the present article. They claim a running time $\tilde{O}(nd)$: however, in their proof, they have a total running of $\mathrm{\Omega }(k^4)$; hence, the claim of near-linear running time is only valid under the assumption $n\ge k^4$ (see their proof of Theorem 3.23).

The $k$-median and $k$-means problems are NP-hard even when the input is in the Euclidean plane ${ℝ}^2$ [38, 36]. However, in low-dimensional spaces, it is possible to compute a $(1+\epsilon )$-approximation, for any $\epsilon >0$, in time $f(\epsilon ,d)\tilde{O}(n)$ [17]. If the target running-time is polynomial in the dimension $d$, the problems becomes NP-hard to approximate: within a factor of 1.015 for $k$-median and 1.06 for $k$-means [20]. The best approximation ratios are $1+\sqrt{2}$ for $k$-median and 5.912 for $k$-means, based on a primal-dual algorithm running in large polynomial time [16]. For faster and practical algorithms, [22] improves the running time of $k$-means++ to almost linear, while roughly preserving the approximation guarantee, and [35] that improves the approximation guarantee to $O(1)$ albeit with a running time of $O(nkd)$. For linear-time algorithm, the one from [11] stays the best we know of. Embedding the input metric into a quadtree yields an expected distortion on distances of $poly(d)\log \mathrm{\Delta }$: hence, combined with dimension-reduction to turn $d$ into $\log k$, with reduction of the diameter to $poly(n)$, one can easily compute a $polylog(k)\log n$-approximation to $k$-median. However, this does not work for $k$-means as the quadtree does not preserve squared distances (see e.g. [17]).

Several sketching techniques are applicable to clustering in Euclidean space: it is possible to reduce the dimension to $O({\epsilon }^{-2}\log k)$ in near-linear time $\tilde{O}(nd)$, while preserving the cost of any clustering up to a multiplicative $(1\pm \epsilon )$ factor. It is also possible to build coresets in time $\tilde{O}(nd{k}^{o(1)}+n\log \log \mathrm{\Delta })$, which reduces the number of distinct points to $O(k{\epsilon }^{-2-z})$ (see [26] for the specific running time,³³3As we noted previously, this paper uses polylogarithmic approximation in near-linear time, which was actually not known; we state their run-time when using an almost-linear time algorithm instead. which uses the coreset algorithms from [23, 21]). Combining those techniques with e.g. [35], it is possible to compute an $O(1)$-approximation to $(k,z)$-clustering in time $\tilde{O}(nd{k}^{o(1)}+k^2)$.

Beyond Euclidean space, $k$-median is NP-hard to approximate within a factor of $1+2/e$ and $k$-means within $1+8/e$ [29].

For the incremental version of $k$-median, the best known approximation ratio is 7.656 for general metric spaces [12] and 7.076 for Euclidean spaces [41]. The approximation ratio cannot be better than 2.01 [12].

Last, Mettu and Plaxton [39] introduced an algorithm related to recursive greedy, for the Facility Location problem. In this algorithm as well they weight balls with a value: [6] showed how to use a value proportional to the number of points in the ball, to estimate the optimal cost of Facility Location problem in a streaming setting.

Instead of ${\mathrm{Value}}_{\mathrm{MP}}(B(x,r))={\sum }_{y\in B(x,r)}{(r-\text{dist}(x,y))}^z$, we will use a function $Value$ approximating $r^z\cdot |B(x,r)|$ as follows: first, count the number of points inside the ball $B(x,r)$ up to a factor $3$, and then multiply by $r^z$. Formally, we show that it is enough to use a function $Value$ that satisfies

This is our key new insight, that allows for fast implementation: it is indeed much easier to approximately count the number of points in a ball than to evaluate ${\mathrm{Value}}_{\mathrm{MP}}$.

In addition to this key simplification, we allow for some approximations in the computation of the different ball considered by the algorithm.

Allowing approximation of balls to redefine children. We say that $N(x,r)$ is a $c$-approximate ball of $x$ at radius $r$ if it satisfies $B(x,r)\subseteq N(x,r)\subseteq B(x,c\cdot r)$. When $c$ is fixed, we simply say that $N(x,r)$ is an approximate ball. The algorithm uses this notion instead of the “child” used in Algorithm 1: it considers all balls (of radius $r/2c$) that are centered at points that are in an approximate ball of $x$ (of radius $r$).

Forbidding balls. To select the starting ball at the beginning of an iteration, we move away from the notion of $isolated$. Instead, our algorithm maintain a set of available balls. Initially, those are all balls; when center $c_j$ is placed at the end of $j$-th iteration, our algorithm forbids (i.e., remove them from the set of available balls) all balls that are too close to $c_j$. More precisely, the algorithm computes, for all $r$ powers of $2c$ such that $1\le r\le \mathrm{\Delta }$, an approximate ball $N(c_j,100c^4\cdot r)$. It removes from the set of available balls all balls of the form $B(p,r)$ with $p\in N(c_j,100c^4\cdot r)$. There are $O(n\log \mathrm{\Delta })$ such balls.

Reducing the number of balls. Perhaps not surprisingly, the radius of balls can be rounded, in order to limit the number of distinct balls to consider. For this, we assume for simplicity that the diameter $\mathrm{\Delta }$ is a power of $2c$. The algorithm will consider only balls of the form $B(x,r)$, where $x\in P$ is an input point and $r$ is a power of $2c$, such that⁵⁵5Although points are at distance at least $1$, it is important for our algorithm to consider balls with smaller radius $1/{(2c)}^7$ in order to be sure that, around any point, there is a ball available unless there is a center at the point. $1/{(2c)}^7\le r\le \mathrm{\Delta }$.

We give the pseudocode of our modified algorithm in Algorithm 2.

Let $B(y,r_y)$ be a descendant of $B(x,r_x)$, and let $x_0=x,\mathrm{\dots },x_{\mathrm{\ell }}=y$ be the centers of the sequence of balls $(B_0,\mathrm{\dots },B_{\mathrm{\ell }})$ from the definition of descendant. For every $i$, we have $\text{dist}(x_i,x_{i+1})\le \frac{10c^2\cdot r_x}{{(2c)}^i}$. By triangle inequality, this implies

$◀$

A consequence of this fact is that any ball appears at most once in an approximate ball of the Center Selection loop:

We split the proof into two parts: first, if a ball appears in an approximate ball, then it is forbidden at the next Forbidding step. Second, if a ball is available, then all its potential descendants are available. Combined, those two results conclude our lemma.

Let $y$ be the center of the ball $B^j$, and let its radius be $r_y=2c\cdot r_x$. Fact 7 ensures that both $x$ and $c_j$ are at distance at most $20c^2\cdot 2c\cdot r_x$ from $y$: therefore, , hence the ball $B(x,r_x)$ is in the set $N(c_j,100c^4{r}_x)$ and is forbidden on line 12 when $c_j$ is selected as a center. $◀$

Let $B(x,r_x)$ be any ball, and let $B(y,r_y)$ be a potential descendant of $B(x,r_x)$. Suppose that $B(y,r_y)$ is not available at the beginning of an iteration of the loop in line 4. Then there exists a center $c_i$ selected by the algorithm such that $y\in N(c_i,100c^4\cdot r_y)$ and therefore $\text{dist}(y,c_i)\le 100c^5\cdot r_y$. Because $B(y,r_y)$ is a descendant of $B(x,r_x)$, we know by Fact 7 that $\text{dist}(x,y)\le 20c^2\cdot r_x$. Using the triangle inequality, we get $\text{dist}(x,c_i)\le \text{dist}(x,y)+\text{dist}(y,c_i)\le 20c^2\cdot r_x+100c^5\cdot r_y$. We also know that $r_y\le r_x/2c$ and therefore $\text{dist}(x,c_i)\le (50c^4+20c^2)\cdot r_x\le 100c^4\cdot r_x$. Hence $x\in B(c_i,100c^4\cdot r_x)\subset N(c_i,100c^4\cdot r_x)$, and $B(x,r_x)$ is also forbidden at the $i$-th iteration. $◀$ Combining those two facts concludes the proof. $◀$

First, by definition, the running time of computing all values in line 2 of the algorithm takes time $T_{Value}$.

For analyzing the Forbidding loop, we note the following: instead of computing $N(c_i,100c^{4}r)$ in line 14, it is merely enough to compute the set of balls that have not been forbidden yet in that set. To do so with the datastructure from the theorem statement, one can merely do the following: for each valid $r$, initialize the datastructure with all points in the ground set. Then, after each computation of an approximate ball $N(c_i,100c^{4}r)$ of the current set, remove all $x\in N(c_i,100c^{4}r)$ from the ground set of the corresponding $r$ in line 15. Hence, the preprocessing time is $T_{Pre}\cdot O(\log \mathrm{\Delta })$, and the total running time is $O((T_N+T_{rm})|P|\log \mathrm{\Delta })$. Indeed, there are $O(\log \mathrm{\Delta })$ different radius $r$, and at each level any point appears at most once before being removed. The guarantees of the datastructure therefore ensure that, at any given level, the total running time is bounded by $O((T_N+T_{rm})|P|)$.

Last, it remains to analyze the Center Selection loop. For this, we can use the datastructure from the lemma, without removing any point: Lemma 8 ensures that each ball $B(x,r)$ appears only once. Hence, the total running time is again $O((T_{Pre}+T_N|P|)\log \mathrm{\Delta }))$. $◀$

This is a direct application of Lemma 12: given $r$ and $\delta =1/n^2$, compute $f_i(p)$ for all $i$ and $p$, in time $O\left(d{n}^{1+3/c^2+o(1)}\right)$. First, to remove a point $p$ from $P$, simply remove it from all the tables $T_i[f_i(p)]$ for $i=1,\mathrm{\dots },\mathrm{\ell }$: this takes time $O(\mathrm{\ell })=O\left(n^{3/c^2}\right)$.

To answer a query given a point $p$, compute $T_i[f_i(p)]$ for all $i$, in time $O(|T_i[f_i(p)]|)$ and define $N(p,r):={\cup }_{i=1}^{\mathrm{\ell }}{T}_i[f_i(p)]$. The running time to compute the union is at most $\mathrm{\ell }\cdot O(|N(p,r)|)=O(n^{3/c^2}|N(p,r)|)$. The first two bullets of Lemma 12 ensure the desired accuracy guarantee. $◀$

Hence, in the vocabulary of Theorem 6, $T_{rm}=T_N=O\left(n^{3/c^2}\right)$. It only remains to compute the values: this can be easily done combining Lemma 12 with the sketching techniques of Lemma 3, as follows:

We show how to compute, for all $p\in P$, an approximate count of the number of points in $B(p,r)$, namely a value $Count(p,r)$ such that $|B(p,r)\cap P|/3\le Count(p,r)\le 3r\cdot |B(p,c\cdot r)\cap P|$. Multiplying $Count$ by $r^z$ gives the lemma.

To build the estimates $Count(p,r)$, the first step of the algorithm is to compute $f_i(p)$, for all $i\in \{1,\mathrm{\dots },\mathrm{\ell }\}$ and all $p\in P$, using Lemma 12 with $r$ and $\delta =1/n^2$. This takes time $O\left(d{n}^{1+3/c^2+o(1)}\right)$. Due to Lemma 12, we have the guarantee that, with probability $1-1/n^2$,

Therefore, it is merely enough to estimate $\left|{\cup }_{i=1}^{\mathrm{\ell }}{T}_i[f_i(p)]\right|$ using Lemma 3 (with $S_{i,u}=T_i[u]$, $t=n$, and for all $p$, $Q_p=\{f_i(p),i=1,\mathrm{\dots },\mathrm{\ell }\}$). As each point $p$ is in at most $\mathrm{\ell }$ sets $S_{i,u}$, and each $Q_p$ has size at most $\mathrm{\ell }$, the running time of the algorithm from Lemma 3 is $O(\mathrm{\ell }n)$

Hence, the overall running time is $O\left(d{n}^{1+3/c^2+o(1)}\right)+O(n\mathrm{\ell }\log (n))=O\left(d{n}^{1+3/c^2+o(1)}\right)$. $◀$

Thus, Theorem 6 implies: