Multi-Swap k-Means++
Authors: Lorenzo Beretta
Vincent Cohen-Addad
Silvio Lattanzi
Nikos Parotsidis
Venue: NeurIPS 2023 poster (2023)
Keywords: Clustering
k-means
approximation algorithms

Abstract: The $k$-means++ algorithm of Arthur and Vassilvitskii (SODA 2007) is often the practitioners' choice algorithm for optimizing the popular $k$-means clustering objective and is known to give an $O(\log k)$-approximation in expectation. To obtain higher quality solutions, Lattanzi and Sohler (ICML 2019) proposed augmenting $k$-means++ with $O(k \log \log k)$ local-search steps obtained through the $k$-means++ sampling distribution to yield a $c$-approximation to the $k$-means clustering problem, where $c$ is a large absolute constant. Here we generalize and extend their local-search algorithm by considering larger and more sophisticated local-search neighborhoods hence allowing to  swap multiple centers at the same time. Our algorithm achieves a $9 + \varepsilon$ approximation ratio, which is the best possible for local search. Importantly we show that our algorithm is practical, namely easy to implement and fast enough to run on a variety of classic datasets, and outputs solutions of better cost.

Multi-Swap k-Means++
Lorenzo Beretta' University of Copenhagen
Vincent Cohen-Addad
Silvio Lattanzi Google Research
Google Research cohenaddad@rooglee com
loremzo@berretta@@aallcoom
silviol@googe com
Nikos Parotsidis Google Research ni nikosp@gogge.. com
Abstract
The k-means++ algorithm of Arthur and Vassilvitskii (SODA 2007) is often the practitioners' choice algorithm for optimizing the popular k-means clustering ob- jective and is known to give an Ollog k) aapproximatio in expectation. To obtain higher quality solutions, Lattanzi and Sohler (ICML 2019) proposed augmenting k-means++ with O log log k) local search steps obtained through the k-means++ sampling distribution to yield c-approximation the k-means clustering problem, where is a large absolute constant. Here we generalize and extend their local search algorithm by considering larger and more sophisticated local search neigh- borhoods hence allowing to swap multiple centers at the same time. Our algorithm achieves a 9 + approximation ratio, which is the best possible for local search. Importantly show that our approach yields substantial practical improvements, we show significant quality improvements over the approach of Lattanzi and Sohler (ICML 2019) several datasets.
Introduction
Clustering is central problem in unsupervised learning. In clustering one is interested in grouping together "similar" object and separate "dissimilar" one. Thanks to its popularity many notions of clustering have been proposed overtime. In this paper, focus metric clustering and of the most studied problem in the area: the Euclidean h-means problem.
In the Euclidean k-mmaas problem one is given in input set of points P in IRd The goal of the problem is find a set of centers so that the sum of the square distances the centers is minimized. More formally, we are interested in finding a set C of l points in B such that Epep mintee cl??, where with ||p cll we denote the Euclidean distance between and e.
The k-means problem has a long history, in statistics and operations research. For Euclidean k- means with running time polynomial in both n, k and d, 5.912 1-pp-ppproiimttiionn was recently shown in Cohen- et al. 2022a]. improving upon Kanungo et al. [2004], Ahmadian et al. (201)) Grandoni al. [2022 by leveraging the properties the Euclidean metric. In terms of lower bounds. the first to show that the high- digensinnnsi k-means problems were APX-hard were Guruswami and
Indyk] (2003). and later (Aaasthi et al. [2015) showed that the APX-hardness holds even if the centers can be placed arbitrarily in _ The inapproximability bound was later slightly improved by Lee et al. 2017 until the recent best known bounds (Cohn-AAdaa and Karthik C. S. [2019) [ohen-Addad et al. 2022d 2021] that showed that it is NP-hard to achieve a better than I (06-approximation and hard to approximate better than 1.36 assuming a stronger conjecture. From a more practical point of view, Arthur and Vassilvitskii 2009) showed that the widely-used popular heuristic of Lloyd Lloyd
Authors are ordered alphabetical order.
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

[957] can lead to solutions with arbitrarily bad approximation guarantees, but can be improved by simple seeding strategy, called k-means++, so as to guarantee that the output is within an (log k) factor of the optimum Arthur and Vassilvitskii [2007].
Thanks its simplicity k-means++ is widely adopted in practice. In effort to improve its performances Lattanzi and Sohler (2019), Choo et al. (2020) combine k-means++ and local search to efficiently obtain a constant approximation algorithm with good practical performance. These two studies show that one can use the k-means++ distribution in combination with a local search algorithm to get the best of both worlds: a practical algorithm with constant approximation guarantees.
However, the constant obtained in Lattanzi and Sohler [2019), [Cho et al. [2020] is very large (several
thousands in theory) and the question as whether one could obtain practical algorithm that would
efficiently match the 9 + e-approximation obtained by the n ((d/e) algorithm of Kanungo et al. [2004]
has remained open. Bridging the gap between the theoretical approach of Kanungo et al. [2004] and k-means++ has thus been a long standing goal.
Our Contributions. We make significant progress on the above line of work.
We adapt techniques from the analysis of Kamungo et al. [2004] to obtain a tighter analysis of the algorithm in Lattanzi and Sohler [2019]. In particular in Corollary we show that their algorithm achieves an approximation of ratio of Rs 26.64.
We extend this approach multi-swaps, where we allow swapping more than one center at each iteration of local search, improving significantly the approximation to R 10.48 in time (nd (A))(
Leveraging ideas from Codan et al. 2021]. we design a better local search swap that
improves the approximation further to 9 + E (see Theorem 12). This new algorithm matches the 9 aapproximmaio achieved by the local search algorithm in Kaanungo et al. [2004] but it is significantly more efficient. Notice that 9 is the best approximation achievable through local search algorithms, as proved in Kanungo al. (2004).
We provide experiments where we compare against k-means++ and Lattanzi and Sohler [2019]. We study a variant of our algorithm that performs very competitively with our theoretically sound algorithm. The variant is very efficient and still outperforms previous work in terms of solution quality, even after the standard postprocessing using Lloyd.
Additional Related Work. We start by reviewing the approach of Kanungo al. 2004] and a possible adaptation to our setting. The bound of 9 + on the approximation guarantee shown by Kanungo et al. [2004] is for the following algorithm: Given a set S of k centers, if there is a set S of at most 2/e points in Rd together with set S of | points in S such that S \ S S+ achieves better k-means cost than S, then set S == \ S U and repeat until convergence. The main drawback of the algorithm is that asks whether there exists set S of points in Rdd that could be swapped with elements of 5 to improve the cost. Identifying such a set, even of constant size, is already non-trivial. The best way of doing so is through the following path: First compute a coreset using the state-of-the-ee the-art coreset construction of Cohen-AAdad et al. [20226) and apply the dimensionality reduction of |eechetti et al. [2019] Makarychev et al. [2019), hence obtaining a set of Ö(k/E") points in dimension O(log k k/e²). Then, compute grids using the discretization framework of Matousek [0000 to identify a set of E od4) ~ kO(e=> log(1/e)) grid points that contains oottimm centers. Now, run the local search algorithm where the sets S are chosen from the grid points by brute-force enumeration over all possible subsets of grid points of size at most, say s. The running time of the whole algorithm with swaps of magnitude s, i.c.: |S+| S s, hence becomes kO(s *log(2/e)) for an approximation of (1 + (9 + 2/s), meaning a dependency in k of k0 3 log(1/=)) to achieve a 9 + c-approximation. Our results improves upon this approach in two ways: (1) it improves over the above theoretical bound and (2) does so through an efficient and implementable, i.e.: practical, algorithm.
Recently, Grunau et al. [003) looked at how much applying a greedy rule on top of the k-means++ heuristic improves its performance. The heuristic is that at each step, the algorithm samples & centers and only keeps the one that gives the best improvement in cost. Interestingly the authors prove that
from theoretical standpoint this heuristic does not improve the quality of the output. Local search algorithms for k-median and k-means have also been studied by Gupta and Tangwongsan (2088 who
drastically simplified the analysis of Arya et al. 2004]. Cohen- and Schwiegelshohn 2017

---

demonstrated the power of local search for stable instances. Friggstad et al. 2019), Codan-Addad al. [019] showed that local search yields PTAS for Euclidean inputs of bounded dimension (and
doubling metrics) and minor-free metrics. Codhen-Addad 20188) showed how to speed up the local
search algorithm using kd-trees (i.e.: for low dimensional inputs).
For fixed k, there are several known approximation schemes, typically using small coresets Becchetti al. 2019], Feldman and Langberg [201T], Kumar et al. 22100. The state-of-th-art approaches
are due to [Bhattcharra et al. [2020]. Jaiswal et al. 2014 The best known coreset construction
remains Cohen-Addan
al. [0022cbb
the constraint on the number of output centers is relaxed, then we talk about bicriteria approxima- tions and k-means has been largely studied Bandyapadhyay and Varadarajan [2016]. Charikar and
Guha 2005 Cohen-Addad and Mathieu (2015). Korupolu et al. [2000). Makarychev et al. 2016
Preliminaries
Notation. We denote with P S RRd the set of input points and let n = [P). Given a point set Q C P we use uu(Q) to denote the mean of points in Q. Given point P and set of centers A we denote with Alp] the closest center in to p (ties are broken arbitrarily). We denote with C the of centers currently found by our algorithm and with O an optimal set of centers. Therefore, given E P. we denote with Clp) and O [p] its closest ALG-center and OPT-center respectively. We denote by cost(Q, A) the cost of points in Q S P w.r.t. the centers in A, namely
cost (e, 4) = Σ qEQ
min 11 cll2. cEA
We use ALG and OPT shorthand for cost(P,C) and cost O*) respectively. When we sample points proportionally to their current cost (namely, sample with probability cost C) costtt, we call this the 2 distribution. When using (-(-) and (2)(.) we mean that & is considered constant. We use ŌJ) to hide polylogarithmic factors in f. The following lemma is folklore.
Lemma 1. Given a point set Q C P and a point p E P we have
cos cost(Q.p) p) = cost(Q, (Q)) + [Q] lp 2(
Let O be an optimal cluster, we define the radius of O as P3 such that o 10; = co,t(O;, 0,), where ((;)). We define the ò-core of the optimal cluster O; as the set of points p E O that lie in a ball of radius (1+ d)p; centered in og. In symbols, core(O;) = P n Bo,,, (1 8)p.). Throughout the paper, is always small constant fixed upfront, hence we omit it.
Lemma 2. Let O be an optimal cluster and sample 9 E O according to the D2 -isttrribtttion restricted to O7. If cossttO > (2 + 35) cost(O;, o:) then Pr[g E core(O?)) 12g(1).
Proof. Define a == co C) /cost(O o;) > + 35. Thanks to Lemma | for each E C we have ||c 0:112 2 (a 1)p? Therefore, for each y E core(O,) and every c E C we have
cost (,, c) = Iyy c|1
as(ap))
Moreover, by a Markov's inequality argument we have 10; core(O s - 110 and thus [core(O;)] 2 Is((O;)). Combining everything we get
cost(core(O)), C) 2 |core(O;) |- min cEC ()))
cost(y, 2a((O)1) 2s(app)
and |0 app cost(O C), hence cost(core(O)) lgg (c(ccos(O
3 Multi-Swap k-Means++
The single-swap local search (SSLS) k-means++ algorithm in Lattanzi and Sohler 2019 works as follows. First, & centers are sampled using k-means++ (namely, they are sampled one by one according to the 2 distribution, updated for every new center). Then, O log log k) steps of local search follow. In each local search step a point g P is Da-sapdeed, then let be the center among
3

---

the current centers C such that cost(P, (C)) u U is minimum. If cost(P. (C) (ct) U (a)) < cost(P, C) then we swap e and 9. or more formally we set C t (C) (c (c])) {q}.
We extend the SSLS so that we allow to swap multiple centers simultaneously and call this algorithm multi-swap local search (MSLS) k-means++. Swapping multiple centers at the same time achieves lower approximation ratio, in exchange for a higher time complexity. In this section, we present and analyse the p-swap local search (LS) algorithm for a generic number of p centers swapped each step. For any constant d > O, we obtain an approximation ratio ALG/OPE = n7² + d where
rf - (2 + 2/p)n (4 + 2/p) = 0.
The Algorithm. First, we initialize our set of centers using k-means++. Then, we run OOddk" local search steps, where a local search step works as follows. We D²~sample a set In = {q1 ... 9p} of points from P (without updating costs). Then, we iterate over all possible sets Out { ... Cp)
of p distinct elements in C U In and select the set Out such that performing the swap (In, Out) maximally improves the cos If this choice of Out improves the cost, then we perform the swap (In, Out), else we do not perform any swap for this step.
Theorem 3. For any 8 > O, the p-swap local search algorithm above runs in (ndk22) time and, with constant probability, finds an (1² + 8appprozzimttinnn of k-means, where y satisfies Equation i
Notice that the SSLS algorithm of Lattanzi and Sohler 2019 exactly the p-swap LS algorithm above for p == 1.
Corollary 4. The single-swap local search in Lattanzi and Sohler 2019 Choo et al. 2020 achieves an approximation ratio- < 26.64.
Corollary 5. For (1) large enough, multi-swap swap local search achieves an approximation ratio < 10.48 in time Ond- poly(k)).
3.1 Analysis of Multi-Swap k-means++
In this section we prove Theorem Our main stepping stone is the following lemma.
Lemma 6. Let ALG denote the cost at some point in the execution of MSLS. As long as 1/2 + local search step improves the cost by a faccoor (21/k) with probability 011/k² ((1/kP-1). 1)
>
Proof of Theorem 3 First, we show that O(k" log log k) local steps suffice to obtain the desired approximation ratio, with constant probability. Notice that local search step can only improve the cost function, so it is sufficient to show that the approximation ratio is achieved at some point in time. We initialize our centers using k-means+++, which gives a Olog ->pappoxximation in expectation. Thus, using Markov's 's inequality the approximation guarantee Ollog k) holds with arbitrary high constant probability. We say that a local-search step is successful if it improves the cost by a factor of at least 1 - ((1kk). Thanks to Lemma 6 we know that unless the algorithm has already achieved the desired approximation ratio then a local -saachh step is successful with probability (211 from (log k) to nj² 6 need O log log k) successful local search steps. Standard concentration bounds on the value of a Negative Binomial random variable show that, with high probability, the number of trial to obtain Ol log log k) successful local-searde steps is Okp log log k). Therefore, after Okp log log k) local-search steps we obtain an approximation ratio of r? + 8.
To go
To prove the running time bound it is sufficient to show that a local search step can be performed in time Ö(ndk"-1). This is possible if we maintain, for each point x E P. dynamic sorted dictionary storing the pairs (cost(x, (a), .) for each C; E C. Then we can combine the exhaustive search over all possible size-p subsets of C U In and the computation of the new cost function using time Ondkk=1 log k)). To do so, we iterate over all possible size-(p - 1) subsets 7 of C U In and update all costs as if these centers were removed, then for each point x P we compute how much its cost increases if we remove its closest center C₂ in (CU Tm) Z and charge that amount to ca- In the end, we consider Out- Z U (c} where is the cheapest-to-removs center found in this way.
The rest of this section is devoted to proving Lemma 6 For convenience, we prove that Lemma holds whenever ALG/OPT > r + O8), which is wlog by rescaling Recall that we now focus on
II In n Out * 0 then we are actually performing the swap (In) \ Out, Out \ In) of size < p.
Also known dynamic predecessor search data structure.

---

given step of the algorithm, and when we say current cost, current centers and current clusters we refer the state of these objects the end of the last local-search step before the current one. Let Oi... Og be optimal clustering of P and let O {o: = n₂ for i = 1... k} be the set of optimal centers of these clusters. We denote with C1 ... Ck the current set of clusters at that stage of the local search and with C = {c1...cu ) ) the set of their respective current centers.
We say that G; captures oj if o is the closest current center to 0j. namely c; C[oz]. We say that is busy if it captures more than p optimal centers, and we say it is lonely if it captures no optimal center. Let O {o= cost(O,, C) > 8 ALGG/} and Ć = C [Cloi] or E E \\ Õ}. For ease of notation, simply assume that O {or ... on and C {or ... c.r). Notice that h' h.
Weighted ideal multi-swaps. Given In C P and Out S Č of the same size we say that the swap (In, Ouf) is an ideal swap if In S O. We now build a set of weighted ideal multi-swaps S. First, suppose wlog that {cr... (c)...c)] } is the set of current centers in C that are neither lonely nor busy. Let L be the set lonely centers C. For each = 1...t, we do the following. Let In be the set of optimal centers in Ô captured by G. Choose a set L; of [In] - centers from L, set (In, Oul) and add it to S. For each busy center C; {cI+1... cy } let be the set of optimal centers Ô captured by C, pick a set L; of A| lonely current centers from L (a counting argument shows that this is always possible). Set + C \ L, For each 0j E A and eg E L; assign weight 1/[[A| - to (oj, (e) and add it to S. Suppose we are left with l centers of .6. of included any de in any swap yet. However, since C O
f. < L L; and define Oul L, U {c:}. Assign weight
O such that Clo & C. Apparently, we have not left with least l > l lonely centers 6 ....... dp C. For each 1... we assign weight to (o,, c4) and add to S.
Observation 7. The process above generates set of weighted ideal multi-swaps swaps such that: (i) Every swap has size at most p: (ii) The combined weights of swaps involving an optimal center o E is (iii) The combined weights of swaps involving current center is at most 1 1/p.
Consider an ideal swap (In, Out). Let 0j = Uo,emn O and Cout = Us, Cy. Define the reassignment cost Reassign gIn, Out) the increase in cost of reassigning points in Cout \ 0 In to centers Out. Namely,
Reassign(In, Out) cost(Cous \ Out) cost(Cout \ C),
We take the increase in cost of the following reassignment as an upper bound to the reassignment cost. For each p E Cout 0j, we consider its closest optimal center O' [p] and reassign p to the current center that is closest to O [p), namely C(O
[p]] In formulas, we have
Reassign(In, Out) S
Σ
cost(p,
cost(p.Cip)
peCout O
cost(p,Cion (p]) cost(p,Clp)
peCout
Indeed, by the way we defined our ideal swaps we have C(O & Out for each p & O and this reassignment is valid. Notice that the right hand side in the equation above does not depend on In.
Lemma 8. EpEp cost(p, (pl) & 20PT + ALG + 2/ALG//PP
Proof. Deferred to the supplementary material.
Lemma 9. The combined weighted reassignment costs of all ideal multi-swaps in S is at most
(2 + 2/p) (OPT + (ALG//PP)
5

---

Proof. Denote by w(In, Out) the weight associated with the swap (In, Out).
Σ
w(In, Out)
Reassign(In, Out) -
(Un,Out)ES
w wInn, Out).
cost(p: costp
(In,Out)ES
pECout
(1+ 1/p)-
Σ Σ
cost(p, cost",/)))ccts(.CCpp)) -
GEE PEC,
(1+ 1/p)
Σ PEP
cost(p,c(0'1) ALG
The second inequality uses (tit) from Observation
7
Applying Lemma 8 completes the proof.
Recall the notions of radius and core of an optimal cluster introduced in Section 2 We say that swap (In, Out) is strongly improving if cost(P, (C In) \ Out) S (1 - 5/k) cost(P,C). Let
In {o1....o o, S Ô and Out = =fI....c. C Ć we say that an ideal swap (In, Out) is good if for every 91 E core(o)) ...4. 9s core(o,) the swap (Q, Out) strongly improving, where
Q - (-)...%s)) We call an ideal swap bad otherwise. We say that an optimal center 0; E O is good if that's the for at least
of the ideal swaps it belongs to, otherwise we say that it is bad. Notice that each optimal center in Ö is assigned to a swap in S, with G the union of cores of good optimal centers in O
either good or bad. Denote
Lemma 10. If an ideal swap (In, Out) is bad, then we have
cost(O C) < (2 6) cost(O]n O²) + Reassign(In, Out) ÖALG/K.
Proof. Let In for 0z). 0z) Q
{9I...9e} 9s} such that q1 core(o1).... 9s core(o.).
Then, by Lemma I cost In: ②) < (2 - ())oost( O*). Moreover, Reassign (nn, Oul) = cost(P (P) 0 C Oul) cost( P Oi because points in P Cout are not affected by the swap. Therefore, cost(P, (C Q) Out) S (2 + 8) cost In: O²) + Reassign [n, Ond) +
cost(P) (P) On Suppose by contradiction that Equation 2 does not hold, then
cost(P,C) cost (P, (CU 2) Out)
cost(P) Of,,C) + cost(O7,..C) - cost (CU Q) Out) 2 8ALG/k.
Hence, Out) is strongly improving and this holds for any choice of Q, contradiction.
Lemma 11. /FLG//PPT>> n² + 5 then cost(G,C) l2g Co((ost()))), Thus, if we D2-sample we have Plg E G] = (1).)
Proof. First, we observe that the combined current cost all optimal clusters in O' \ Ö is at most k jALG/k kk 8ALG. Now, prove that the combined current cost of all Of such that o is bad is < (1 ALG. Suppose, by contradiction, that it is not the case, then we have:
(1 26)A
Σ Bad o.Eō
cost(Or,
Σ
w(In, Oat) cost(O],..))
Bad (Jn,Out)ES
w(Tn, Out) ((2 + 8)cost(O) 0²) Reassign(Tn, Out) 5ALG/KK) <
OIn,Out)
(2 + 8) + (2 22+2/))OPT + (2 +
BALG.
The second and last inequalities make use of Observation 7 The third inequality uses Lemma 10
Setting 7 ALL/OPT we obtain the inequality ,72 (2 + 2/p + 0())) (4+)/==0))) (2++))/==0(0))((4+))))-)))00 Hence, we obtain a contradiction in the previous argument as long as 7) (4 + 2/p + O()) 0. A contradiction there implies that at least an 5-fraction of the current cost due to points in Gaond o O,, We combine this with Lemma and conclude that the total current
cost of G UGeod o, COre(O;) is 2s (s(cost(P,C))

---

Finally, we prove Lemma Whenever 91 E G we have that 91 E core(o;) for some good 01 Then, for some s S p we can complete on with 02 0,. such that In = {or... o,} belongs to a good swap. Concretely, there exists Oul C C such that (In, Out) is a good swap. Since In Ô we have cost(O,,C) > 8OPT/k: for all o; E In, which combined with Lemma 2 gives that for i = 2.... Pl.: E core(o:)] 2 (1/k.)). Hence, we have Plu: E core(o;) for I...s] 2 2.. (1/kp-1)) Whenever we sample @1 ...... 2, from core(o1)... ... core(os), we have that (2, Oul) is strongly improving. Notice, however, that (2, Oul) is s-swap and we may have < p. Nevertheless, whenever we sample q1 ...9. 9, followed by any sequence 9s+1 ... Ap it is enough to choose Out' = Out U [9=+1.... ... Jpp to obtain that (191...pp)) Out') is an improving p-swap.
4 A Faster (9 + -++ Approximation Local Search Algorithm
The MSLS algorithm from Section 3 achieves approximation ratio of 17² + where n²² - (2 2/p)n - (4 + 2/p) = 0 and ε > 0 is an arbitrary small constant. For large p we have 1j ~ 10.48. On the other hand, employing p simultancous swaps, Kanungo et al. [2004] achieve an approximation factor of E² + where 5² (2 + 2/p) - (3 + 2/p) = O. If we set p R I/c this yields a (9 + ((e)- approximation. In the same paper, they prove that 9-approximation is indeed the best possible for p-swap local search, if p is constant (see Theorem 3.1 in Kanungo et al.] (2004)). )) They showed that 9 the right locality gap for local search, but they matched it with a very slow algorithm. To achieve (9 + (0++e)-app --aapoooximttion they discretize the space reducing to O(ne-a) candidate centers and perform an exhaustive search over all size-(1/e) subsets of candidates every step. As saw in the related work section, it is possible to combine techniques from coreset and dimensionality reduction to reduce the number of points to n' = k poly(( and the number of dimensions d' log ks E
This reduces the complexity of Kanungo et al. 2004] to kO(E
Plogs" ')
In this section, we leverage techniques from Cohen-Addad et al. 2021 to achieve (9- e)-
approximation faster In particular, we obtain the following.
Theorem 12. Given set of n points in Ie with aspect ratio A. there exists algorithm that computes a 9 + c-approximation to k-means in time o(= log Ole
(0)='((()) ((lee))()) 2 .
(A)
poly(e="
Notice that, besides being asymptotically slower, the pipeline obtained combining known techniques is highly impractical and thus it did not make for an experimental test-bed. Moreover, it is not obvious how to simplify such an ensemble of complex techniques to obtain a practical algorithm.
Limitations of MSLS. The barrier we need to overcome in order to match the bound in Kanungo et al. [2004] is that, while only consider points in P as candidate centers, the discretization they employ considers also points in IRd \ P. In the analysis of MSLSS we show that we sample each point 4 from core(O)) or equivalently that q: E B(o,, (1 + <)p:), where P3 is such that Of would have the same cost w.r.t. o; if all its points were moved on a sphere of radius p; centered in og. This allows us to use a Markov's 's inequality kind of argument and conelude that there must be (1O points in O n (o.., (1 + e)p:). However, we have no guarantee that there is any point all in O n B(o,, (1 =)p.). Indeed, all points in OF might lie on aB(o, pi). The fact that potentially all our candidate centers q are at distance at least p; from 0; yields (by Lemma []) cost(O) 2 Scost(O; o,). which causes the zero-degree term in g² (2+2/p) 23+2/p)
from Kanungo et al. 2004 to become 4 in our analysis.
Improving MSLS by taking averages. First, we notice that, in order to achieve (9 + e)-
approximation we need to set p ((1/s). The main hurdle to achieve (9 19+ (-p=--appoxiiiaatonn is that we need to replace the 9; in MSLS with a better approximation of o;, We design a subroutine that computes, with constant probability, an e-approximation o; of o, (namely, cost(O), ó,) < (1 + #)cost(O;; o.)). The key idea is that, if sample uniformly O(1/E) points from O and define 6, to be the average of our samples then cost(O), 6,) I (1 + s)cost(O;, o.)
Though, we do not know Of, so sampling uniformly from it is non-trivial. To achieve that, for each 4 we identify a set N of nice candidate points in P such that poly(s)/k /k fraction of them are from
O We sample (11E)) points uniformly from N and thus with probability (e/k) O(1/)) we sample only points from O Thus far, we sampled (1//)) points uniformly from N n O7 What about
"The complexity in Theorem 12 can be improved by applying the same preprocessing techniques using coresets and dimensionality reduction, similar to what can be used to speed up the approach of Kanungo et al. 20044. Our algorithm hence becomes asymptotically faster.
7

---

the points in O; \ N? We can define N so that all points in O; \ N are either very close to some of the (4)) or they are very far from qi- The points that are very close to points (g)); are easy to treat. Indeed, we can approximately locate them and we just need to guess their mass, which is matters only when 2 poLy and we pay only a log logo(//)(1//
multiplicative overhead to
guess the mass close to qj for j = 1...P P = (1/e). As for a point f that is very far from qi (say, Ilf 1-- q:ll > pi) notice that, although f's contribution to cost(O;, 0s) may be large, have cost(f. o) s cost(f.o: for each o E P;) S B(o, (2 + E)p;) assuming q; E core(o.)
Experiments
In this section, we show that our new algorithm using multi- swap local search can be employed to design an efficient seeding algorithm for Lloyd's which outperforms both the classical k-means++ seeding and the angle-swap local search from Lattanzi and Sohler (2019).
Algorithms. The multi-swap local search algorithm that we analysed above performs very well in terms of solution quality. This empirically verifies the improved approximation factor of our algorithm, compared to the single-swap local search Lattanzi and Sohler [019]]
Motivated by practical considerations, we heuristically adapt our algorithm to make it very competitive with SSLS in terms of running time and still remain very close, in terms of solution quality, to the theoretically superior algorithm that we analyzed. The adaptation of our algorithm replaces the phase where it selects the centers to swap-out by performing an exhaustive search over (*+P) subsets of centers. Instead, we use an efficient heuristic procedure for selecting the centers to swap-out, by greedily selecting one by the centers to swap-out. Specifically, we select the first center to be the cheapest one to remove (namely, the one that increases the cost by the least amount once the points in its cluster are reassigned to the remaining centers). Then, we update all costs and select the next center iteratively. After p repetitions we are done. We perform an experimental evaluation of the "greedy" variant of our algorithm compared to the theoretically-soun algorithm from Section B and show that employing the greedy heuristic does not measurably impact performance.
The four algorithms that we evaluate are the following: 1) KM++: The k-means++ from Arthur and Vissilvitskii [207]. 2) SSLS: The Single-swap local search method from Lattanzi and Sohler [2019], 3) MSLS: The multi-swap local search from Section 3 and 4) MSLS-G: The greedy variant of multi-swap local search as described above.
We use MSLS- SGp x and MSLS-p z. to denote MSLS-G G and MSLS with p z, respectively, Notice that MSI.S-G-p = 1 is exactly SSI.S. Our experimental evaluation explores the effect of p-swap L.S, for p > 1, in terms of solution cost and running time.
Datasets. We consider the three datasets used in Lattanzi and Sohler [019 to evaluate the perfor- mance of SSLS: 1) KDD-PHY 100,, 000 points with 78 features representing quantum physic task kdd [2004] 2) RNA- - 488, 565 points with 8 features representing RNA input sequence pairs Uzilov et al. 2006]] and 3) KDD-BIO 145, 751 points with 74 features measuring the match between protein and a native sequence kdd [2004). We discuss the results for two or our datasets, namely KDD-BIO and RNA. We deffer the results KDD-PHY to the appendix and note that the results very similar to the results on RNA.
We performed preprocessing step to clean-up the datasets. We observed that the standard deviation of some features was disproportionately high. This causes all costs being concentrated in few dimensions making the problem, in some sense, lower-dimensional. Thus, we apply min-max scaling to all datasets and observed that this causes all our features' standard deviations to be comparable.
Experimental setting. All our code is written in Python. The code will be made available upon publication of this work. We did not make use of parallelization techniques. To run our experiments, used a personal computer with 8 cores, 1.8 Ghz processor, and 15.9 GiB of main memory We run all experiments 5 times and report the mean and standard deviation in our plots. All our plots report the progression of the cost either w.r.t local search steps, or Lloyd's 'ss iterations. We run experiments on all our datasets for k = 10, 25, 50. The main body of the paper reports the results for k 25, while the rest can be found in the appendix. We note that the conclusions of the experiments for k = 10, 50 are similar to those of k- = 25.
Removing centers greedily. We first we compare MSI.S-G with MSL.S To perform our experi- ment, we initialize k 25 centers using k-means++ and then run 50 iterations of local search for both
8

---

KDD-BIO, 25
RNA, 25
140
0.25
035
Cost
Cost
025
0.30
080
0.55
HOL5 Gipd
HOLS-=2
1ELSp=3
RM1
9111
LS llerations
LS lleraaioos
Figure 1: Comparison between MSLS and MSLS-G, for p = 3. for l = 25, on the datasets KDD- BIO and RNA. The y axis shows the solution cost divided by the means solution of KM+++
algorithmss, for p swaps. Due to the higher running of the MSLS we perform this experiments on I 1% uniform sample of each of our datasets. We find out that the performance of the two algorithms is comparable on our instances, while they both perform roughly 15%-27% at convergence. Figure shows the aggregate results, over repetitions of our experiment.
may happen that MSLS, which considers all possible swaps of size p at each LS iteration, performs worse than MSLS-G G-G as a sub-optimal swap at intermediate iterations may still lead to better local optimum by coincidence. Given that MSLS-G G performs very comparably MSLS, while it much faster in practice, we use MSLS-G for the rest of our experiments where we compare to baselines. This allows us to consider higher values of p, without compromising much the running time.
our first experiment we run KM+++ followed by 50 iterations of MSLS-C with 1, 4, 10 and plot the relative cost w.r.t. KM++ each iteration, for 25. The first row Figure plots the results. Our experiment shows that, after 50 iterations MSLS-G for p = 4, 7, 10 achieves improvements of roughly 10% compared to MSLS-G-p = 1 and the order of 20% 30% compared to KM++. We also report the time per iteration that each algorithm takes. For comparison, we report the running time of single iteration of Lloyd's next to the dataset's name. It is important to notice that, although MSLS-G-P
Results: Evaluating the quality and performance the algorithms.
is faster, running more iterations MSLS-G-P 1 is not sufficient to compete with MSLS-G when p > 1.
Evaluating the quality after postprocessing using Lloyd. In second experiment, we use KM++ and MSLS-G as a seeding algorithm for Lloyd's and measure how much of the performance improvement measured in the first experiment retained after running Lloyd's. First, we initialize our centers using K.M++ and the run iterations of MSLS- for = 1, 7. We measure the cost achieved by running 10 iterations of Lloyd's 's starting from the solutions found by MSLS-G as well as +++ In Figure ( (second row) we plot the results. Notice that, according to the running times from the first experiment, 15 iterations iterations of MSI.SG take less than 10 iterations of Lloyd's for = 4, 7 (and also for = 10, except on RNA). We observe that MSLS-G for p > performs at least as good SSLS from Lattanzi and Sohler 2019 and in some cases maintains non-trivial improvements.
Results:
Evaluating the quality and performance of the algorithms against a fixed deadline. In this experiment run KM++ followed by MSLS-G with p = 1, 4, 10, for of fixed amounts of time. This setting allows the versions MSLS-G with smaller swap size to perform more iterations compared the versions the algorithm with a larger swap size, smaller swap size leads to lower running time per iteration. Let 7 be the average time that Lloyd's algorithm requires complete simple iteration on specific instance. We plot the cost of the solution produced by each algorithm after running À x for each E {1.... 20} in Figure 3 Our experiment shows that MSLS-G for p 4, 7, 10 achieves improvements of more than 5% compared MSLS-G-1 = 1 even when compared against fixed running time, and of the order of 20% 30% compared to KM++
Results:
Conclusion and Future Directions
present a new algorithm for the k-means problem and we show that it outperforms theoretically and experimentally state-of-the-ar art practical algorithms with provable guarantees in terms of solution

---

KDD-BIO, k 25 (0.19 per Llayd's iter)
RNA, k = 25 (0.56 s per Lloyd's iter)
1.0
2.05
200
0.55
IX+
Cost
Cost
M615-Gpp11
0.90
MS
M515 1..5:-0p-10..0
0.80
0.25
0.00
30
20
LS Iterations
LS iterations
KDD-BIO, 25
RNA,
37000
24500
MM
5-0p21
Me Selppl
34000
VSL5-G0-7 7
NSLS-G-p-7
VSL5-G0-19
MBLS-Gio-10
22000
20000
Cost
Cost
23000
30900
22500
22006
31900
ll Iterations
Lloyd's Iterations
Figure 2: The first row compares the cost of MSLS-G, for E (1,4, 7, 10]. divided by the mean cost of KM++ at each LS step, for = 25. The legend reports also the running time of MSLS-G per LS step (in seconds). The second row compares the cost after each of the 10 iterations of Lloyd with seeding from MSLS-G, for p E {1, 4, 10) and 15 local search steps and K.M++, for k = 25.
KDD-BIO, = 25
KDD-PHY, k 25
110
1.00
0.95
191+
0.55
Cost
MSL5-G-p=1
MSL5-G-pm4
NSLS-Gpm4
0.90
M5L5-0-pp>
NSLS-G-p77
MB15-6-p-10
NSLS p
0.80
0.85
0.0
20
Timeout/ (time for 1 Lloyd's iter.)
Timeout / (time for 1 Lloyd's iter.)
Figure 3: Comparison of the cost produced by MSLS-G, for p {1, 4, 10) and k 25 the datasets KDD-BIO and KDD-PHU, divided by the mean cost KM++ after running for fixed amount of time in terms of multiplicative factors to the average time for an iteration of Lloyd's algorithm (i.e., for deadlines that are 1x...., 20x the average time of an iteration of Lloyd).
quality. A very interesting open question is to improve our local search procedure by avoiding the exhaustive search over all possible size-p subsets of centers to swap out, concretely algorithm with running time (12
This work was partially done when Lorenzo Beretta Research Student at Google Research. Moreover, Lorenzo Beretta receives funding from the European Union's Horizon 2020 research and innovation program under the Marie Sktodowska-Curie grant agreement
Acknowledgements.
No. 801 199.
10

---

References
Kdd cup. 2004.
Sara Ahmadian, Ashkan Norouzi-Fard, Ola Svensson, and Justin Ward. Better guarantees for k- means and Euclidean k-median by primal-dual algorithms. SIAM Journal on Computing, 49(4): (OCSI77-9-FFOCSIT--1 2019.
David Arthur and Sergei Vassilvitskii. K-means++ the advantages of careful seeding. In Proceedings of the eighteenth annual ACM- MMMM symposium on Discrete algorithms, pages 1027-1035, 2007.
David Arthur and Sergei Vassilvitskii. Worst-case and smoothed analysis of the ICP algorithm, with an application to the k-means method. SIAM J. Comput. 39(2)766-7882, 2009. doi: 10.1137/ 070683921. URL http:::/a...ol..rr/11........ 111///00013770006832111
Vijay Arya, Naveen Garg, Rohit Khandekar, Adam Meyerson, Kamesh Munagala, and Vinayaka Pandit. Local search heuristics for k-median and facility location problems. SIAM Comput.
33(3):544-562, 2004, doi: 10.1137/5009753900411666 URL https: org/org///// 2tgg10.1137 50097539702216402
Pranjal Awasthi, Moses Charikar, Ravishankar Krishnaswamy, and Ali Kemal Sinop. The hardness of approximation of euclidean k-means. In Lars Arge and János Pach, editors, 31st International Sym- posium Computational Geometry, SoCG 2015, June 22-25, 2015, Eindhoven, The Netherlands, volume 34 of LIPics, pages 754-767. Schloss Dagstuhl Leibniz Leibniz-Zentrum Zentrum für Informatik, 2015. doi: URL //oii/////////////////// org/ 10. 4.300222IIIcc. www.sSOCG... 2015. 1..45.330 754
Sayan Bandyapadhyay and Kasturi Varadarajan. On variants k-means clustering. In 32nd International Symposium Computational Geometry (SoCG 2016). Schloss Dagstuhl-Leibniz- Zentrum fuer Informatik, 2016.
Luca Becchetti, Mare Bury, Vincent Cohen-Addad, Fabrizio Grandoni, and Chris Schwiegelshohn. Oblivious dimension reduction for k-means: beyond subspaces and the johnson-lindenstee lemma. In Proceedings of the 51st Annual ACM SIGACT Symposium Theory of Computing, STOC 2019, Phoenix, AZ, USA, June 23-26, 2019, pages 1039-1050, 2019. doi: 10.1145/3313276 3316318. URL httph::/ //oii org/ 10. 11111111133111 3316318
Anup Bhattacharya, Dishant Goyal, Ragesh Jaiswal, and Amit Kumar. On sampling based algorithms for k-means. In Nitin Saxena and Sunil Simon, editors, 40th IARCS Annual Conference on Foundations of Software Technology and Theoretical Computer Science, FSTTCS 2020, December 14-18, 2020, BITS Pilani, K K Birla Goa Campus, Goa, India (Virtual Conference), volume 182 of LIPlcs, pages 13:1-13:17, Schloss Dagstuhl Leibniz-Zentrum für Informatik, 2020, doi: 10.
1.300LIPPIc..SSTTC.. 200.CC.000...... 13. URL https: :/ //wii org/10..//00 10..1..///IIIc... PTT.IIICCSSSTTTTTSSC FSTTCS 2020 13
Moses Charikar and Sudipto Guha. Improved combinatorial algorithmss for facility location problems.
SIAM J. Comput., 34(4):803-824, 2005. doi: 10.1137/500975337701999939 URL https
//doi.orr/// 0tte://0..11///099777970001533
Davin Choo, Christoph Grunau, Julian Portmann, and Vaclay Rozhon. k-means++: few more steps yield constant approximation. In Hal Daumé and Aarti Singh, editors, Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 1909-1917. PMLR, 2020. URL htttps:/pproceediii choo20a. html
mlr. press/v119/
Vincent Cohen-Addad. A fast approximation scheme for low-dimensional k-means. Artur Czumaj, editor, Proceedings of the Twenty-Ninth Anmual ACM-SIAM Symposium Discrete Algorithms, SODA 2018, New Orleans, LA, USA, January 7- 2018, pages 130-440 SIAM, 2018 doi: 10.1 13.1137/1.00 10.1137/1.978161199
29. WRLLhttts::///
org/ 10. 1137/1 9Tt1://1..77111117788888 29
Vincent Cohen-Addad and Karthik C. S. Inapproximabiliiy of clustering in Ip metrics. In David Zuckerman, editor, 60th IEEE Annual Symposium on Foundations of Computer Science, FOCS 2019, Baltimore, Maryland, USA, November 9- 2019, pages 519-539, IEEE Computer Society, 2019. doi: 10. 10. 00040. URL https: ///oii///oo// org/ 10. 1109/F0CS 00040
11

---

Vincent Addad and Claire Mathieu. Effectiveness of local search for geometric optimization. In 31st International Symposium on Computational Geometry, SoCG 2015, June 22-25, 2015, Eindhoven, The Netherlonds, pages 329 1999-333,22 2015. doi:
URL.
http://åx.dos: doi. org/ 4230/L1PIcssss SOCG. 2015. 329
Vincent Cohen-Addad and Chris Schwicgelshohn. On the local structure of stable clustering instances. In Chris Umans, editor, 58th IEEE Annual Symposium on Foundations of Computer Science, FOCS 2017, Berkeley, CA. USA, October 15-17, 2017. pages 49-60. IEEE Computer Society, 2017. doi: N0.1009FFCC.2201... 14. URL/httss: //doi. org 10011199FFOOC 2017. 14
Vincent Cohen-Addad, Philip N. Klein, and Claire Mathieu. Local search yields approximation schemes for k-means and k-median in euclidean and minor-free metrics, SIAM J. Comput. 48(2); 644-667, 2019. doi: 01.113//7M1127111
URL URL
1137/17M112717X
Vincent Cohen-Addad, Karthik C. S., and Euiwoong Lee. On approximability of clustering problems without candidate centers. In Dániel Marx, editor, Proceedings of the 2021 ACM SIAM Symposium on Discrete Algorithms, SODA 2021, Virtual Conference, January 10 13, 2021, pages 2635- 2648. SIAM, 2021. doi: 10.1137/1.9781111111 178111197761166115555 htRlsttps:://de org/10. 1137/1 9781611976465. 156
Vincent Coden-Addad, David Saulpic, and Chris Schwicgelshohn. Improved coresets and sublinear algorithms for power means in euclidean spaces. Advances in Neural Information Processing Systems, 34:21085-210998 2021.
Vincent Cohen-Addad, Hossein Esfandiari, Vahab S. Mirrokni, and Shyam Narayanan. Improved approximations for euclidean k-means and k-median, via nested guasi-independent sets. In Stefano Leonardi and Anupam Gupta, editors, STOC '22: 54th Annual ACM SIGACT Symposium on Theory of Computing, Rome, Italy, June 20 - 24, 2022, pages 1621-1628. ACM, 2022a. doi: 10.1 0el114///51993533555 35500111. URL (tRLittpss:///ot //cttps:://o...o. org /0/ 1149//5199//9 3520001
Vincent Adden-Addad, Kasper Green Larsen, David Saulpic, and Chris Schwiegelshohn. Towards optimal lower bounds for k-median and k-means coresets. In Stefano Leonardi and Anupam Gupta, editors, STOC '22: 54th Anmual ACM SIGACT Symposium on Theory of Computing, Rome, Italy, June 20 24, 2022, pages 1038-1051. ACM, 2022b. doi: 10.1145/55193353519996 URL
https://doi org/ 10. 1145/3519935 11444//331993335599966
Vincent Cohen-Addad, Kasper Green Larsen, David Saulpic, Chris Schwiegelshohn, and Omar Ali Sheikh-Omar. Improved coresets for euclidean k-means. In
NewrlPS, 2022c. URL http://papers. cac.ips.ccpppppppp faaar_filee files/papeess http//pppp/////aaatttee haash/ htw.:/bbbeaaaaaaaaaaa2aaaaaeee--
hhmll
Vincent Cohen-Addad, Euiwoong Lee, and Karthik C. S. Johnson coverage hypothesis: Inapproxima- bility of k-means and k-median in l, metrics. In Proceedings of the 2022 ACM-SIAM Symposium on Discrete Algorithms, SODA 2022. SIAM, 2022d.
D. Feldman and M. Langberg. A unified framework for approximating and clustering data. In STOC, pages 569-578, 2011.
Zachary Friggstad, Mohsen Rezapour, and Mohammad R. Salavatipour. Local search yields a PTAS for k-means in doubling metrics. SIAM J. Comput., 48(2):452 480, 2019. doi: 10.1137/ 17M1127181. URL https:://oot. org/ 10. wttp:/..1117/77111111 1137/17M1127181
Fabrizio Grandoni, Rafail Ostrovsky, Yuval Rabani, Leonard J. Schulman, and Rakesh Venkat. A refined approximation for euclidean k-means. Inf. Process. Lent. 176251,,,, 2022. doi: .b.0016f.pp1..222.1063
httttp:://do..rgg// 1016/j ipl hpw.pp2.1.......... 2022. 106251
Christoph Grunau, Ahmet Alper Özüdogru, Václav Rozhon, and Jakub Tetek. A nearly tight analysis of greedy k-means+++. In Nikhil Bansal and Viswanath Nagarajan, editors, Proceedings of the 2023 ACM-SIAM Symposium Discrete Algorithmss, SODA 2023, Florence, Italy, January 22- 25, 2023, pages 1012-1070. STAM, 2023. doi: 10. 1111111779911111 1977554.ch39, URL https //doi 1ttp://0..177//.77111111997977/9 ch39
12

---

Anupam Gupta and Kanat Tangwongsan. Simpler analyses of local search algorithms for facility location. CoRR, 00000009.25544 2008. URL httpp::/aaax...org//bb...s wbwgg//bb/0000/0/ 2554
Venkatesan Guruswami and Piotr Indyk. Embeddings and non-approximability of geometric problems. In Proceedings of the Fourteenth Annual ACM- -SMM Symposium on Discrete Algorithms, January 12-14, 2003, Baltimore, Maryland, USA., pages 537-538, 2003. URL http://dl acm. org/
citation.cfmiid
644198
Mary Inaba, Naoki Katoh, and Hiroshi Imai. Applications of weighted voronoi diagrams and randomization to variance-based k-clustering, In Proceedings of the tenth annual symposium on Computational geometry, pages 332-339, 1994.
Ragesh Jaiswal, Amit Kumar, and Sandeep Sen. A simple D 2-sampling based PTAS for k-means and other clustering problems. Algorithmica,, 70(1):22-46, 2014, doi: 10.1007700043--01--9339.9 URL https://doi. org/10.
Tapas Kanungo, David M. Mount, Nathan S. Netanyahu, Christine D. Piatko, Ruth Silverman, and Angela Y. Wu. A local search approximation algorithm for k-means clustering. Com- putational Geometry, 28(2):89-112 2004. ISSN 0925-7721. doi: http:://6i..or//00.000 Lcomgee.2004.000 0..0000.33.03 $0925772104000215 Special Issue on the 18th Annual Symposium on Computational Geometry SooC22002
URL https: htmps://wwww //ww. sciencedirect. wwwpeeetttcc/ccieeieeeeeeeee http///ceenneaaticell//e/e
Madhukar R. Korupolu, C. Greg Plaxton, and Rajmohan Rajaraman. Analysis of a local search heuristic for facility location problems. J. Algorithmss, 37(1); (46()):466-88)) 188, 2000. doi: 10. 10.60/m/j...
2000.1 htRt:/tp:///x doi org/ 10. /000/jaa..... 2000. 1100
Amit Kumar, Yogish Sabharwal, and Sandeep Sen. Linear-time approximation schemes for clustering problems in any dimensions. J. ACM, 57(2):5:1-5:333 2010. doi: 10.1145/1667053316660000 URL
https://doi. org/ 10. 1141/111700000 1667054
Silvio Lattanzi and Christian Sohler. A better k-means++ algorithm via local search. In Kamalika Chaudhuri and Ruslan Salakhutdinov, editors, Proceedings of the 36th International Conference on Machine Learning, volume 97 of Proceedings of Machine Learning Research, pages 3662-3671. html
PMLR, 2019. WRRLittts///proeeedingg mlr pww.ss//9//aaaaaa.a
Euiwoong Lee, Melanie Schmidt, and John Wright. Improved and simplified inapproximability for k-means. Inf. Process. Letr. 120:40 43, 2017. doi: 10.10166f..p..2016.11 11.009. URL https
//doi org/10....... 1t.11...../..................... 1p1 2016. 11.009)
SP Lloyd. Least square quantization pem. bell telephone laboratories paper. published in journal much later: Lloyd, sp: Least squares quantization in pem. IEEE Trans. Inform. Theor 1957//192), 18, 1957.
Konstantin Makarychev, Yury Makarychev, Maxim Sviridenko, and Justin Ward. A bi-criteria approximation algorithm for k-means. In Klaus Jansen, Claire Mathieu, José D. P. Rolim, and Chris Umans, editors, Approximation, Randomization, and Combinatorial Optimization. Algorithmss and Techniques, APPROX/RANDOM 2016, September 7-9, 2016, Paris, France, volume 60 of LIPIcs, pages 14:1-14:20. Schloss Dagstuhl LLibniz-Zentrum für Informatik, 2016. doi: 10.4230/LIPIcs 14. httL:htps:://dtl http://t....oo////////////// 4230/LIPICs APPROX -RANDOM
2016. 14
Konstantin Makarychev, Yury Makarychev, and Ilya P. Razenshteyn. Performance of johnson- lindenstrauss transform for k-means and k-medians clustering. In Moses Charikar and Edith Cohen, editors, Proceedings of the 51 Annual ACM SIGACT Symposium on Theory of Computing, STOC 2019, Phoenix, AZ, USA, June 23-26, 2019, pages 1027-1038. ACM, 2019. doi: 10. 10.1145/3313276 3316350. URL htttps://doii oraidoi..oog/goo
/txg//10.10 10. 1145/33132/// 3316630
Jirí Matousek. On approximate geometric k-clustering. Discrete & Computational Geometry. 24(1):61-84 2000. doi: 10. 00.100775004444000019,9 5004540010019
URL http::/dx..de doi. org
1trg/10.10 1007/
Andrew V Uzilov, Joshua M Keegan, and David H Mathews. Detection of non-coding mas on the basis of predicted secondary structure formation free energy change. BMC bioinformatics, 7(1): 1-30, 2006.
13

---

Supplementary Material
Proofs from Section 3
Lemma 8. Eper cost(p, (p]) 20PT+ ALG 2VALGVOPP
Proof.
cost(p,
PEP
cost(p, Clo.])
O PEO?
Σ
10 cost(o,, Clo.)) cost(tt
OPT +
cost cost(O"lp),com ||," [p])
PEP
OPT
cost(("(p)) Clp)) S
PEP
OPT
(110 [p] pll + 11p Clpl
PEP
2OPT + ALG +
Σ 1001p
1|p, Clpll < 2OPT + ALG 2V OPT.
PEP
The second equality is due to Lemma and the last inequality is due to Cauchy-Schwaan
Proofs from Section 4
In this section, we prove the following.
Theorem 12. Given set of n points in Rd with aspect ratio A, there exists algorithm that computes a e-approximation to k-means in time odkkOs=²) log Ole "A) 2 p-iylee
We start with key lemma showing that sample of size O(1/e) is enough approximate 1-mean.
Lemma 13 (Form Inaba et al. (1994)). Given an instance C Rd. sample m = 1/(E0) points uniformly random from P and denote the set of samples with S. Then cost(P, (()) (1 e)cost(P, P.p(P)) with probability at least 1 6.
8 we have I|p(S) p(P) 12 < ecost(P, p(P)) /IP Then, applying Lemma gives the desired result. First, we notice that ((P) is an unbiased estimator ((P). namely E ETu(S)] p(P). Then, have
Proof. We want prove that with probability
[s|
[||s - cost(P, ((P())
[llx(S) -
m Σ
m PP
where are uniform independent samples from P. Applying Markov's inequality concludes the proof.
The algorithm that verifies Theorem 12 is very similar to the MSLS algorithm from Section and we use the same notation describe it. The intuition is that in MSLS we sample Q {4) ... yp} hoping that qa core(o;) for each i; here we refine q to a better approximation û; of o; and swap the points (o:); rather than (g:); Our points are generated taking the average some sampled point, thus possibly have o; P while, on the other hand, e P.
14

---

A approxxmaani MSLS algortihm. First, we initialize our set of centers using +-mesns Then, ndkkO' 9), 2poty(s +) local search steps, where local search step works follows. Io of points from P (without updating costs).
Set p = ( We De-sample a set Q = [g]
Then, we iterate over all possible sets Out {cr ................................ cp} of p distinct elements in C' UU Q. We define the set of temporary centers T = -C Q) \ Out and run a subroutine APX-CENTERS (T) which
returns a list of poly(e=') log ((g) size-s sets In {ōr (61...6.) ...,} (where s = IQ \ Out)) We select the set I in this list such that the swap (In, Out 2) yields the maximum cost reduction.
(A)
Then select the set Out that maximizes the cost reduction obtained in this way. If (In, Out Q)
actually reduces the cost then we perform that swap.
A subroutine to approximate optimal centers. Here we describe the subroutine APX-CENTERS(T). Let Q Out = (-11..................... Recall that s S = 0 This subroutine out- puts a list of 20olv( "1 llog O(= ) (A) size-s sets In = {ō1 ... 6. Here we describe how to find a list of 20oly(e ") >-log(A) A) values for o₁. The same will apply for 02 ó. and taking the Cartesian product
((==)(()))
yields list of 2poly(="") log log (A) size-s sets. Assume wlog that the pairwise distances between points in lie in [1, A]. We iterate over all possible values of p1 E {1, (1++)))) ... (L [log] AlJ. We partition P in three sets: the set of far points F = {z E P cost(x, q)) > pi/e"), the set of close points C = {a E P \ F cost(1z T) S 5 pî) and the set of nice points N P (CU F). Then, we sample uniformly from N a set S of size ((e-!). For each (s 1)-tuple of coefficients
ao, a, E { (1 - e), (1 s)² ..... (1 (1-c8 E) E) (e")- }u(0} we output the candidate solution given by the convex combination
61 = 0 (,0......0 ..,) co0((S) + Ei= 0gii
2i 0i
so, for each value of pi, we output 20oly(s values for 01 Hence, 2polu(e 1) (A) values in total.
Analysis
The key insight in the analysis of the MSLS algorithm form Section 3 was that every q was proxy for o, because qi E core(o;), and thus q provided good center for O In the analysis of this improved version of MSLS we replace q with o which makes a better center for O Formally, fixed Out, we say that a point o; is a perfect approximation of o, when cost(O;, , (C Ú (6.)) Out) s
(1+ =)OPT, + =OPT/k. /k. We define Ô and Č as in Section 3 except that we replace 8 with E (which here is not assumed to be constant). Likewise, we build the set S of ideal multi-swapp Recall that we say that a multi-swap (In, Out) is strongly improving if cost(P, P, (C (CU (CUTT) \ Out)
in Section 3
(1 - e/k) cost(P,C) Let In = {o....o.} },} S Ô and Out = (ci...Ca c. S C, we overload the definition from Section 3 and say that the ideal multi-swap (In, Oul) is good if for every
In {ó.... } such that each ô, is a perfect approximation of o, for each = l...s the swap (In, Out) is strongly improving. We call an ideal swap bad otherwise. As in Section 3 we define the core of an optimal center; once again we replace with e, which is no longer constant. The two following lemmas our stepping stones towards Theorem 12
Lemma 14. If ALG/OPT > 9 + O(s) then, with probability k k-0(e=") polulee="" there exists Out C ECUC Q such that:
(i) If Q Out {q1 .q....q.) q.} then qi E core(o:)... E core(o,) for some ....., E O
(ii) If we define T = {01... ... .0) then (In, Out Q) is good ideal swap.
Lemma 15. If (1) from Lemma 14 holds, then with probability k k00 ?) 2- 2-polyls 1), the list returned by APX-CENTERS contains Tn = (Ø1... such that ô, is perfect approximation of o for each i 1...s
Here we prove that our improved MSLSS algorithm achieves a + (())- approximation, which is equivalent to Theorem 12 up to rescaling E. Combining Lemma 14 and
Proof of Theorem 12
15

---

Lemma 15 we obtain that, as long ALG/OPT > 9 + O(=), with probability at least k 00e -((e))) 2 polyyy the list returned by APX-CENTERS contains In {Ôi... ...).) such that (In, Out @) is strongly improving. If this happens, we call such a local step successful. Now the proof goes exactly as the proof of Theorem 3 Indeed, We show that (0(s 22oly(s local steps suffice to obtain [ log log k/E) successful local steps, and thus to obtain the desired approximation ratio, with constant probability.
(01...0.
To prove the running time bound it is sufficient notice that local search step can be performed in time nd logo(=
(())))))()))) (A) (A) 2
In the rest of this section, we prove Lemma 14 and Lemma [5
Observation 16. If we assume 5 non-constant in Lemma explicitly we obtain Pr(q E (o 2 poly(E).
then performing the computations
In order prove Lemma 14] we first prove the two lemmas. Lemma 17 is the analogous of Lemmal 10
and Lemma 18 is the analogous of Lemmal] 11 Overloading once again the definition from Section 3
we define G as the union of cores of good optimal centers in O, where an optimal center is defined to
be good least one of the ideal multi-swaps S belongs to good (exactly as Section
Lemma 17.
an ideal swap (In, Out) is bad, then have
cost(Oīm: C) S (1 + +) cost(Oin. O*) + Reassign Ceassign(In,.Ol Out) + cALG/k.
(4)
Proof. Let In = for ... J.). ..)) In {ōr ... ] such that o; is perfect approximation of o, for each i 1...8. Recall that 0i = U O then
cost Cest(onn (C (CUIn) Out) $ cost(O), (Cu Out) - (1 (1+ =)cos((O,, O*). (5)
Moreover, Reassign(In, Out) cost(P)) Oj,, C Out) - cost(P) (P) O7, C) because points
in P Cout are not affected by the swap. Therefore, cost P. (C ( \ Out) s (1
e) cost(OO 0*) Reassign(In, Out) + cost(P Ofn:C) C). Suppose by contradiction that Equa- tion does not hold, then
cost(P,C) cost P. (CUTT)) \ Out)
(est(P\Onns cost(O),,.)) - cost (P, (CU \ Out) Z eALG/k
Hence, (In, Out) is strongly improving and this holds for any choice of In, contradiction.
Lemma 18. ALG//PT OPT > 9 + O()) then cost(G.C) 2 cost(P.C) poly(E). Thus, if D2 -ammle we have Plg E Gl 2 poly(E).
Proof. First, we observe that the combined current cost of all optimal clusters in O \ Ö is at most
eALG/k = eALG. Now, we prove that the combined current cost all O such that or is bad < (1 2a) ALG.. Suppose, by contradiction, that it is not the case, then we have:
- AE) < Σ
cost(O;,C) <
Σ
w[In, Out) cost(O),...)
o,Eō
Bad OUn,Os
we[n, Out) ((1 + cooot((i O*) + Reassign(In, (ALG//)
Bad OInt
(1 + =) + + 2/p))OOPPT + /LGG EALG.
The second and last inequalities make use of Observation 7 The third inequality Lemma 17
Setting 2 ALG/OPTT we obtain the inequality 7/² - (2 2/p + 0(e))) (3 S
Hence, we obtain contradiction in the previous argument as long as 17" (2 (
16

---

O/)= > O, which holds forp ((==²) and nz² = 9 ++0 A contradiction there implies that least an E-fraction of the current cost is due to points in UGood o, O Thanks Observation I6
we have P grrcostly.C)19 E core(O;) q E O,] 2 poly(E). Therefore, we can conclude that the current cost of G U@ood (Ore(O?) is at least a poly(E)-fractin of the total current cost.
Proof of Lemma 14] Thanks to Lemma 18 we have that Plgi E EG] 2 poly(E). Whenever qi E G have that qi core(oi) for some good on. Then, for some 8 S p we can complete or with 02...00 such that In = {or- ... belongs to a good swap. Concretely, there exists Oul C C such
that (In, Out) is good swap. Since In C O we have cost(O), C) > cOPT/k for all E In, which combined with Observation 16 gives that, for each i = 22.. 8, Plg: E (ore 2 poly(e)/k. s) 2 2- 2--ool( x k00(=)), Notice, however, that (In, Oul) is s-swap and we may have p. Nevertheless. whenever we sample 9i ......9.. followed by any sequence 1s+1... it is enough to choose Oul' = Oul U (1s+1 ...9p) to obtain that (p)): Oul') is improving p-swap.
Hence, we have P(qs- E core(o,) for = 1 s]
In order to prove Lemma 15 we first need few technical lemmas. Lemma 19 (Lemma 2 from Lattanzi and Sohler (2019)) For each x, v, 2 E Rd and & > 0,
cost(z, v) S 51++))oost(n 2) + + c/s)cott(2))))
Lemma 20. Given q Rd and Z S Rd such that cost (Z. q) E then, for each E Rd
(1)) c(s))coot(Z2))) o) O S | co (1
7, + O())
Proof. To obtain the first inequality, we apply Lemma 19 to bound cost(z, o) < (1++)costi( (1 + 1/E)cost(2, 9) for each E Z. To obtain the second inequality, we bound cost(q.0) (1 (1+=)cost(e o) + (1 1/E)cost(2, q) for each & E Z.
Lemma 21. Let X (z) be weighted set of points in B such that x; has weight wa
Let f be the weighted average of X. Let i = p(an ...))) be the weighted average of where T, has weight ag If w; Sa, w1/(1 for each i 1...f, l, then if we interpret cost(X, C) as Ex,x W; cost(r,, C) we have cost(X, S (1 c(s)))os((XX,..)))
Proof. We note that js minimizes the expression cost (X, Moreover, cost(X, 2) S Ei= a cost(a; 2) S cost(X,=) /[1 e)). Since 인 minimizes the expression Σi= 03 cost(ri, z) it must be cost(X,, pi) S cost(X, 14) /(1 e).
Adopting the same proof strategy, we obtain the following.
Observation 22. Thanks to Lemma 20 we can assume that the points in Z are concentrated in for the purpose of computing a (1 +1++0 O(E)) -aproximation to the 1-means problem on Z, whenever an
additive error I is tolerable. Indeed, moving all points in Z to q introduces a +(E) multiplicative error on cos (ost(Z..) and a (e)) additive error.
The next lemma shows that point 2 that is far from center experiences small variation of cost(2, o) when the position of o is slightly perturbed.
Lemma 23. Given E Rd such that l|o =|| 2 r/e we have that for every of E B(o, r), cost(2,o') (1 + (E)))
Proof. It is enough to prove it for all o' that lie on the line L passing through o and 2, any other point in o" E r) admits a point o' E B(o, r) O L with I1o - =|| = ||o" - = is enough to acoet(..)
compute the derivative of cost(z,)) -) with respect to the direction of L and see that ar. B(orr) = (1=))r))r/s Thus, cost(2,00 o') = cost(2, o) + (1 + 0(E))r2/E = (1 + ((E)) o).
Proof of Lemma 75 Here we prove that for each 0 ..... o, there exist coefficients ag as (1) [, (1 2). ....... (1 E) [og 1e"11} U (0} such that the convex combination ô; = à, (i) ... .....") ()))
is perfect approximation of o₁, with probability --0 p--oov(s="), Wlog, show this
17

---

for 01 only. Concretely, we want to show that, with probability k- O- 2 2--ool((= there exist coefficients ag .... a, such that ôN = ő₁ (ag.... ...)..).)))))).. satisfies cost(O (CU {ó1}) Out) v
(1 + ((s))OPT, + ((e)OPT/k. Taking the joint probability of these events for each i = 1...8
we obtain the success probability k- 2 poly(e""), Note that we are supposed to prove that cost(Oj, (C U {or)) \ Out) < (1++))PPT + <OPT/k, however we prove weaker version where E is replaced by O(e). which is in fact equivalent up to rescaling E.
Similarly to Cl:] and [+[-] define T[p] as the closest center to in T. Denote with C1, Fi and N1 the intersections of 0f with C, F and N respectively. In what follows we define the values of cg og..--0, that define ô1 = ô1 (ag...... 0s)) and show an assignment of points in Of to centers in (CU {or }) \ Out
with cost (1 + (s)))TPT O(e)OPT/k. Recall that we assume that qi E core(o:) for each i i= 1. s.
In what follows, we assign values to the coefficients (a:);. It is understood that if the final value we choose for O; is v then we rather set a; to the smallest power of (1 E) which larger than v, if u > =². Else, set a; to O. We will see in the end that this restrictions on the values of o; do not impact our approximation.
In what follows, we will assign the points in Of to C Out, if this can be done inexpensively. it cannot, then we will assign points to 01. In order to compute good value for ôj we need an estimate of the average of points assigned to ô1. For points in N1, computing this average is doable (leveraging Lemmal 13) while for points Of \ N1 we show that either their contribution is negligible or we can collapse them so to coincide with some q Q without affecting our approximation. The coefficients (0x))2) represent the fraction of points in O₂ which is collapsed to th. Oig represents the fraction of points in Of which average we estimate as u(S). Thus, Equation (3) defines o; as the weighted average of points qi where the weights are the (approximate) fractions of points collapsed onto q, together with the the average (SS) and its associated weight ag.
Points in Ct. All points p E Ci such that Tip] g Q can be assigned Typ C \ Out incurring total cost of at most e O by the definition of C1. Given a point E C1 with Tip E Q we might have Tip & C \ Out and thus we cannot assign p to Tpp.. Denote with W the set of points p with Tip Q. Our goal is now to approximate (WW). ). In order to do that, we will move each p E W to coincide with 9: = T(p). We can partition W into Wi... W, so that for each 2 W; T qi If E Z; then we have IP - 9:112 < e²pî. Hence, thanks to Observation 22 we can consider points in W; as if they were concentrated in 4 while losing at most an additive factor OP)) and multiplicative factor (1 + E) on their cost. For i = 1.... set a, + |Wl//O0]) In this way, En a; x / Σi= (g is an approximates solution to 1-mean on W up to a multiplicative factor (1 (1+e) e) and an additive factor OPs)
Points in No. Consider the two cases: (i) cost(N1, T) OOPT/K; (M) cost(N1, T) s OOPTTKk
Case (i). We show that in this case (SS) is (1 + ='aappooxinttion for 1-mean on No. with probability --0(e-)). 2- -ooly(="),, First, notice that if we condition on S S N₁ then Lemma 13 gives that (uS) is (1 11++6 appppoxxxintaooonn for 1-mean on N₁ with constant probability. Thus, we are left to prove that S C N₁ with probability k-O(e-"). ..te"")22-00 2-poly(=-)) We have that the P PprvcoetppTT) (p E Ni P E N] > 2²/k, however the costs w.r.t. T of points in N varies of at most a factor poly(s-1) ), thus Pp-Unirlp E Ni P E N] Z poly(=)/k. The probability of S C Ni is thus (poly(=)/k)
22-pooll==)) In this case, we set ap + |Nil/1Of because p(S) approximates the mean of the entire set N1.
k --
Case (ii). Here we give up estimating the mean of N1 and set orp + 0. The point N1 such that T[x] & Q can be assigned to T[x] incurring a combined cost of 22OPT /h. We partition the remaining points in N1 into Zı U... ... Ze where each point is placed in Z, if T[x] = q-- Now, collapse the points in Z, so as to coincide with q, and show that this does not worsen our approximation factor. In terms of coefficients (ax):. this translates into the updates a; <- a; Z for each 1 1...8.
Indeed, using Observation 22 we can move all points in Z, to g incurring an additive combined cost of :OPT/k and a multiplicative cost of 1 + O(e).
18

---

Points in Fj. Points in Fi are very far from 91 and thus far from 01, hence even if their contribution
cost(Oi.o) 0 might be large, we have cost FF 01) (1 (())) o') for all of in a ball
of radius p1 /s centered in o,, thanks to Lemma 23
Let H be the set of points that have not been assigned to centers in C') Out. In particular, H WUN if points in N1 satisfy case (i) and H- = W U Z1 Ze if points in N₁ satisfy case (ii). We consider
two cases.
If [lx(H) - 21ll S p/e, then [l((H) - 01 S p(1 + + 1/E) because q1 E core(o)). Since for each f E Fi have Ilf -0-1 01ll 2 fl n1 (1+ E)p (((p==) then cost(f.d') = 01) for each of in ball of radius (p/E)) centered in 01, and so in particular for o' = ((H)). Thus in this case we can simply disregard all points in F1 and computing 01 according to the (az); defined above yields a perfect approximation of o;-
(1 (()))
p/e, similar argument applies to show that cost (H, = (1 + E)cost(H,,) for each of in ball of radius (p)) centered in 0,. Indeed, we can rewrite cost cost(H,o) o') as [M] cost(u(!l)) od) + cost(u(H), ), .f.. If [Ix(H)) - qı || < p/z the first term varies of at most factor (1 +)) and the second term is constant. Thus in this case 0 = qi is a perfect approximation of and we simply set 02 and a, 0 for & 1. In other words, here n(N1 U H) is too far from 91 (and thus 01) to significantlyt influence the position of oN and the same holds for any point in F1. This works, of course, because we assumed 91 core(o1).
Else, if []x(H) all
Discussing the limitations on the coefficients values. we were allowed to set a; to exactly the values discussed above, representing the fractions of points from O, captured by different G.S. However, to make the algorithm efficient we limit ourselves to values in [. (1 - -).... (1 - [o01 (e]1 01} {0). Lemma 21 shows that as long as the values of (as); estimate the frequencies described above up to a factor (e)) then the approximation error is within a multiplicative factor + (s))
The proof above would work smoothly
We are left to take care of the case in which a; is set to value e². We set a when dealing
with points in C₁ U N1 and for each x E C₁ U Ni we have, for each o' E B(q,, (1 + s)p). cost(ar o') < 2cost(g1 o') + 2-2cost(2x q1) = ((pp==6), Thus, if we simply set a; + 0 whenever we have a; < E" then the combined cost of points in Of with respect to o' varies by e² 101] Pie (s)) Effectively, ignoring these points does not significantly impact the cost. hence solving 1-mean ignoring these points finds a (1 + O -apprxxmmae solution to the original problem.
Additional Experimental Evaluation
In this section we report additional experiments which presentation did not fit in the main body. In particular, we run experiments on the dataset KDD-PHY and for k = =10,50 50.
In Figure we compare MSLS-G G with MSLS. To perform our experiment, we initialize k = 25 centers using KM++ and then run 50 iterations of local search for both algorithms, for p 3) swaps. We repeat each experiment 5 times. For ease of comparison, we repeat the plot for the KDD-BIO and RNA datasets that we present in the main body of the paper. Due to the higher running of the MSLS we perform this experiments on 1% uniform sample of each of our datasets. We find out that the performance of the two algorithms is comparable on all our instances, while they both perform roughly 15%-27% better than k-means+++ at convergence.
In Figure 5 we run KM++ followed by 50 iterations of MSLS-G with = 1, 4, 10 and k 10, 25, 50 (expcluding the degenerate case p = k = 10) and plot the relative cost w.r.t. KM++ at each iteration. The results for k = 25 KDD-BIO and RNA can be found in Figure 2 We repeat each experiment& times. Our experiment shows that, after 50 iterations MSLS-G G for p = 4, 7, 10 achieves improvements of roughly 5- 10% compared to MSLS-G-p = 1 and of the order of 20% - 40% compared to KM+++ These improvements are more prominent for = 25, 50. We also report the time per iteration that each algorithm takes. For comparison, we report the running time of a single iteration of Lloyd's s next to the dataset's name. Notice that the experiment on RNA for k 50 is performed on 10% uniform sample of the original dataset, due to the high running time.
In Figure 6 we use KM++ and MSI.S- as seeding algorithm for Lloyd's and measure how much of the performance improvement measured is retained after running Lloyd's. First, we initialize
19

---

KDD-BIO, k
KDD-BIO, k= 25
1005
1.00
Cost
Cost
0.92
090
0.50
AS1S-Cpp2
085
MSLS0pp2
MSLS-p22
MSLSp=3
104
106-1
LS Iterations
LS lterations
KDD-PHY, k 25
KDD-PHY, = 25
1.00
HSLS p-2
MGLS
1.05
PSLS
Gp
12111
8h
1.50
140
0.25
Cost
0.00
0.0
0.75
19
LS Iterations
LS lterations
RNA,, = 25
RNA.A 25
100
239
Cost
2.92
Cost
285
875
MP0
0.7
134+-
194+
LS Iterations
LS Iterations
Figure 4: Comparison between MSLS and MSLS-G, for 2 (left column) and = 3 (right column), for k 25, on the datasets KDD-BIO (first row). KDD-PHY (second row) and RNA (third row). The y axis shows the mean solution cost, over the 5 repetitions of the experiment, divided by the means solution cost of KM+++
20

---

KDD BIO, 10 (5.35
Lloyd's iter)
KDD BIO, 50 (24.57 per Lloyd's iter)
1.000
0.92
082
0,900
Cost
0.90
890
0.900
MGLS-Gi-e,, 4.21
18++
0.050
MS.S.6-p-. 12.07
26.55-ppm4, 18.65
M065-65m4, 2.35
S.Gapm3. 22.00
MG.S-Gxp- 13.80
4:6-p-11 23.08
0.75
LS Iterations
LS Iterations
RNA, k = 10 (17.29 per Lloyd's iter)
RNA, = 50 (10.29 per Lloyd's iter)
HSLS
Nbub 5./09 60399
HSL5-66p-44 2
4
15.41
23:6
Cost
0.40
0.65
PSLL5 0-0-1 9485
MbL5 o p-1. 1/105
MBL5:Crpp-) 8.09
HBL5-Crp-é 6
0.7
Iterations
Iterations
KDD-PHY, 10 (5.97 s per Llayd's iter)
KDD-PHY, = 25 (14.69 s per Lloyd's iter)
GM
MoL5 lo p- 5:22
M5L5 o p-1, 3:02
G-L5-G-p-6 1.315
Cost
MSL5--0 o pw?: 13.15
0/ 2.77
3/33300 9 0-10) 25955
MSLS O pei, 3255
0.7
20 Iterations
Iterations
KDD PHY, 50 (25.51 per Lloyd's iter)
100
HG15-Ge1)
WSL5-C-p44
MR15-p-7 25.59
0.5
M515-Gpp-00 17.4
0.90
Cost
0.7
LS Iterations
Figure 5: We compare the cost of MSLS-G, GG, for p E (1,4, 10), divided by the mean cost of K.+++ each LS step, for l: {10, 25, 50), excluding the degenerate case also the running time of MSLS-G per LS step (in seconds). The experiments were run on all datasets: KDD-BIO, RNA and KDD-PHY, excluding the case of l 25 for KDD-BIO and RNA which reported in the main body of the paper.
10. The legend reports
21

---

our centers using KM++ and the run 15 iterations of MSLS-G for 1,4,7. We measure the cost achieved by running 10 iterations of Lloyd's starting from the solutions found by MSLS-G as well as ++++ We run experiments for k = 10, 25, 50 and we repeat each experiment 5 times. We observe that for l = 25, MSL.S- for p > performs at least as good as SSL. from Lattanzi and Sohler [2019] and in some cases maintains non-trivial improvements. These improvements are not noticeable for 10; however, given how Lloyd's behave for k 10 we conjecture that & 10 might be an "unnatural" number of clusters for our datasets.
22

---

KDD-BIO, k 10
KDD-BIO, k == 50
40000
MM++ HRl MBL5-6-9-7
MitS-ppp NSLSOp-?
Cost
42500
12500
$U500
Lloyd's Iterations
Lloyd's Iterations
RNA,
RNA, 50
17500
W++++
M+++
MSLS-Orpz
MSL5-0pp-
12000
HBL5-0-0-7
40000
Cost
003000
36000
15300
32003
Lloyd's Iterations
Lloyd's Iterations
KDD PHY, k 10
KDD PHY. k 25
N5L5-G-11
155000
700000
NBLS-Opea
4515.p=7
152500
HSLS-0-2m44
He2.5-0-p->
105000
Cost
147502
142503
185000
190000
137500
Lloyd's Iterations
Lloyd's Iterations
KDD PHY, k = 50
128000
124030
Cost
118000
134000
Lloyd's Iterations
Figure 6: We compare the cost after each of the 10 iterations of Lloyd with seeding from MSLS-G
for p E 41,4,77 7, 10} and 15 local search steps and KM++, for E {10, 25, 50). We excluded the degenerate case p le 10, and the experiments reported in the main body of the paper.
23


Reviews:
[{"review_id": "o22mQ97TwV", "reviewer": "Reviewer_swvB", "summary": "The following results are given in the paper: \n1. A tighter analysis of the local search algorithm of Lattanzi and Sohler. The paper shows a constant approximation guarantee for their algorithm.\n2. The paper extends the approach of Lattanzi and Sohler to multi-swap local search (where more than one point is swapped in every iteration). Using this approach, an approximation guarantee of 10.48 is obtained.\n3. The paper also gives a (9+\\eps) approximation. Such a 9-approximation algorithm using local search is already known (https://www.sciencedirect.com/science/article/pii/S0925772104000215). The algorithm given in this paper has a slight running time advantage.\n\nThe algorithm is a local search algorithm. The justification of the title is that the swap-in points in the local search strategy are obtained by D^2 sampling (i.e., the sampling strategy used in k-means++).", "questions": "- Some of the other queries are mentioned above within the other fields.", "limitations": "This work is mostly theoretical. There are no negative societal impacts.", "rating": 7, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 4, "strengths": "1. Knowing a better analysis for an already published algorithm is good.\n2. Extension to multi-swap and its analysis giving better approximation guarantee is interesting.\n3. Matching the best possible approximation (possible using local search) using a slightly better running time is also interesting.\n4. The experimental section shows that the local search iterations do better than Lloyd's. This is interesting and can be of practical use. It would have been better if the running time comparison was also given to make the tradeoffs more visible.", "weaknesses": "1. The paper looks good from the point of view of extensions to the Lattanzi and Sohler work. However, from the point of view of developments in local search algorithms (which is essentially what the paper is), the improvement over (https://www.sciencedirect.com/science/article/pii/S0925772104000215) seems incremental. Even the proof techniques are similar.\n2. There are some discrepancies in the experimental section. For example, Figure 2 says \"first row...\", but there is none. I did not find the running time for the local search iterations (even MSLS-G).", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their time and effort put into reviewing the paper. \n\nWe acknowledge that the proof technique is similar to the one by Kanungo et al. on a high level. However, the exhaustive local search technique of their paper leads to their algorithm being prohibitive for practical applications. The inefficiency of the Kanungo et al. paper serves as the main motivation of our work that sharply improves the running time to $O(n\\cdot d \\cdot poly(k))$ for an algorithm with approximation ~10.48 for any constant number of swaps, which is further implementable and gives rise to a very efficient heuristic. Our more elaborate $9+\\epsilon$ approximation algorithm essentially shows that one can match the approximation factor of Kanungo et al. by leveraging coresets and obtaining an algorithm that is faster for a wide range of values of k.\n\nBy \"first row\", we refer to the first row of figures, as the two top figures. We kindly note that the running times of the algorithms (per iteration) are reported in the legends of Figure 2 following the algorithm names. However, we will clarify both these points in the final version of the paper."}, {"review_id": "4M7rFJDDPR", "reviewer": "Reviewer_FVuS", "summary": "This paper proposes a new k-means algorithm: multi-swap local search (MSLS) which combines local search and k-means++, to achieve a constant approximation guarantee with efficient time complexity. Specifically, the local search framework includes a step that selects alternative centers for optimizing the cost function, and this paper selects those centers using the $ D^2 $-sample from k-means++  (without updating costs). \n\nThe authors prove that the random candidate centers obtained in this way lead to an improvement over the cost. By iterating on this step, MSLS will return a constant approximation result.\n\nThe author also leverage ideas from coreset and  dimensionality reducion to propose a more efficent $ 9+\\varepsilon $-approximation algorithm within the local search framework.", "questions": "- Both k-means++ and local search are not limited to the Euclidean space. Can thm3 be generalized to the metric space? I notice that the analysis are influenced by [1] and there exists mean computing ($ \\mu(\\cdot) $) that prevents such generalization. However, [2] gives a constant approximation for k-means in the metric space. Therefore, I am interested in what the current analysis method lacks for the generalization and whether such limitation can be overcomed? \n\n\n\n[1] A Local Search Approximation Algorithm for k-Means Clustering\n\n\n[2] Simpler Analyses of Local Search Algorithms for Facility Location", "limitations": "-     ", "rating": 7, "confidence": 4, "soundness": 4, "presentation": 3, "contribution": 4, "strengths": "- The paper is well written. All sources of motivations and techniques involved in this work are clearly presented, making it easy to follow the key ideas of the algoirthm and the related works. It does make a novel contribution to the classic problem.\n- The idea behind the MSLS algorithm is natural and the proofs are very clear. Although the author refer techniques from the previous work, the approximation analysis is non-trivial and results in very tight ratios.\n- The MSLS algorithm is simple and easy to implement. The experiments are comprehensive. This work is valuable in both theory and practice.", "weaknesses": "- the time complexity is $ \\tilde{O} (nd k^p) $, which grows exponentially w.r.t. $p$.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their time and effort put into reviewing the paper.\n\nWe would like to clarify that some exponential dependence on p is likely needed for any local-search based algorithm. In particular, if one could devise an exhaustive local search algorithm with polynomial dependence on $p$, then setting $p=k$ one would get a polynomial-time algorithm to solve this NP-hard problem that is known to be hard to approximate within a factor of 1.06 [1]. Notice that our algorithm is a $(1+\\epsilon)$ factor away from the approximation achieved by the exhaustive local search by Kanungo et al. which is analyzed for a fixed number of swaps. This means that if we were able to devise an algorithm with sub-exponential dependence on $p$, for $p=k$, we would hope to achieve a $(1+\\epsilon)$ factor approximation to the k-means problem, which is NP-hard. However, it is a very interesting and open question to achieve a dependence on $c^{f(p)}$, for a constant $c$.\n\nIt is indeed an interesting open question on whether one can achieve a similar result in general metric spaces. We tried to extend the analysis from [2] to $k$-means but we find one main obstacle: their projection lemma (Lemma 2.4 in their arxiv version) does not hold for $k$-means because distances squared do not satisfy triangle inequality. The corresponding lemma in our work is Lemma 8 whose proof crucially relies on Lemma 1, which in turn only works for euclidean spaces.\n\n[1] Vincent Cohen-Addad and Karthik C. S. Inapproximability of clustering in lp metrics. In FOCS'19.\n[2] Gupta, Anupam, and Kanat Tangwongsan. \"Simpler analyses of local search algorithms for facility location.\" arXiv preprint arXiv:0809.2554 (2008)."}, {"review_id": "bJ7eSNJQ6y", "reviewer": "Reviewer_2sgw", "summary": "This paper studies the standard $k$-means problem. Given a $k$-means instance $(P,k)$, the goal is to find a set $C$ of centers with size at most $k$ such that the sum of the squared  distances from $P$ to $C$ is minimized. For the $k$-means problem, Lattanzi and Sohler (ICML 2019) proposed an elegant combination of the local search and the $k$-means++ seeding methods. Instead of enumerating all the swap pairs constructed between the data points and the current centers, they showed that $k$-means++ sampling can be used for sampling a data point to serve as the candidate center for swapping in such that a $509$-approximation can be achieved in expectation in time $O(ndk^2loglogk)$. This paper extends the single-swap strategy to multi-swap strategy. In a single local search step, this paper shows that using $k$-means++ to sample $t$ data points induces a good $t$-swap for the local search process such that the clustering cost can be reduced significantly by at least a $\\Omega(1-1/k)$ fraction with certain probability. After $O(k^{t-1})$ local search steps, an improved approximation can be achieved related to the swap size $t$, where the total running time can be bounded by $\\tilde{O}(ndk^{2t})$. When $t$ is large enough, the approximation could be smaller than 10.48. By combining the techniques from Cohen Addad et al. (NeurIPS 2021), the papers shows that by using more local search steps (with exponential dependence on $\\epsilon$), the approximation guarantee can nearly match the lower bound of the $9+\\epsilon$ of the standard multi-swap local search method.", "questions": "1. During the construction of ideal multi-swaps, can the construction process described in the paper guarantee that all the optimal centers in $\\tilde{O}$ is used at least once for constructing ideal multi-swaps (see the example of weakness 2). If some of the optimal clustering centers   are not used for constructing ideal multi-swaps, will the bound for summation of clustering cost for all \"bad'\" optimal clusters in Lemma 11 still hold?\n\n2. In Lemma 2, this paper shows that if the current clustering cost for some optimal cluster $O_i^*$ is larger than $(2+3\\delta)$ times its optimal clustering cost, the clustering cost of data points close to $o_i^*$ should take a large fraction (related to $\\delta$) of the whole clustering cost induced by $O_i^*$. Hence, in each local search step, the success probability for performing a good swap should be related to parameter $\\delta$, for example $\\Omega(1/\\delta)$. However, the running time of the proposed multi-swap local search algorithm is independent of $\\delta$. I think there should be at least an $\\frac{1}{\\delta}$ term in the running time, or does it mean that $\\delta$ is a constant given as the input, please explain.\n\n3. In previous multi-swap local search approximation schemes, the approximation ratio is given as a function related to the swap size $t$ (i.e., in [1], the approximation is given as $(3+2/t)^2$). However, in this paper, the approximation ratio is determined by an inequality related to $\\eta$ and $t$. Can the author provide the approximation ratio in the form similar to that of [1] such that one could figure out how large that swap size is enough for obtaining an approximation guarantee of 10.48.\n\n4. Why are all the instances in the experimental setting small? In the recent results for clustering algorithms, instances of size over 1 million [2] or even 100 million [3] have been considered. How does the proposed multi-swap local search method scale on large datasets?\n\n5. What are the performances of different local search algorithms within a certain time limitation (large enough for reaching the convergence)?\n\n[1] Kanungo T, Mount D M, Netanyahu N S, et al. A local search approximation algorithm for k-means clustering[C]//Proceedings of the eighteenth annual symposium on Computational geometry. 2002: 10-18.\n\n[2] Ren J, Hua K, Cao Y. Global Optimal K-Medoids Clustering of One Million Samples[J]. Advances in Neural Information Processing Systems, 2022, 35: 982-994.\n\n[3] Matsui Y, Ogaki K, Yamasaki T, et al. Pqk-means: Billion-scale clustering for product-quantized codes[C]//Proceedings of the 25th ACM international conference on Multimedia. 2017: 1725-1733.\n\n", "limitations": "Since this is a theoretical paper, I don't think there is potential negative societal impact of this work. ", "rating": 3, "confidence": 4, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "The strengths of this paper can be summarized as follows:\n1. This paper proposes a fast multi-swap local search method for solving the standard $k$-means problem, which runs in linear time in the data size. The approximation ratio significantly improves the previous one with linear running time in the data size (i.e., 509-approximation)\n2. This paper gives another multi-swap local search method which better approximates the optimal clustering centers during the local search swaps using a subroutine called APX-CENTERS, where an improved approximation guarantee, i.e., $9+\\epsilon$, can be obtained with running time exponentially dependent on $poly(\\epsilon^{-1})$. The approximation ratio matches the approximation lower bound of local search methods for the $k$-means problem and significantly improves the running time of previous work based on direct enumeration of swap pairs constructed between the whole dataset and the current centers opened.\n3. In experiments, this paper proposes a heuristic method which avoids the exhaustive searching process for determining which subset of the sampled set of data points should be swapped in. The experiments show that the multi-swap local search method achieves better performance on clustering quality with fixed local search iterations compared with single-swap local search methods.", "weaknesses": "1. Although the authors show that sampling-based multi-swap local search can improve the approximation to a very small constant and even matches the lower bound of local search methods, the core idea behind is to use the successive $k$-means++ sampling strategy to construct candidate set of centers for swapping in, which is an extension of the sampling-based local search method (denoted as LS++ method )proposed by Lattanzi and Sohler. The claimed improvement in ratio seems to be minor as LS++ already guarantees constant approximation (in [1] the authors claimed that they did not attempt to optimize the constants and this paper verifies that LS++ can yield approximation ratio smaller than 26.4). The key idea behind the analysis for multi-swap local search method is as follows: (1) the authors further divide the optimal clusters into different groups according to their current clustering cost (i.e., $cost(P_h^*,C) \\ge \\delta ALG/k$ or $cost(P_h^*,C) < \\delta ALG/k$) (2) the authors give weights for swap pairs such that each current center can be used for at most $(1+1/t)$ times when performing a summation of all bad swap pairs. The analysis for probability lower bound that induces a good swap is similar to that of LS++ method.\n\n2. Some theoretical details is unclear to me. During the construction of ideal weighted multi-swaps, this paper partitions the optimal clustering centers and current clustering centers based on $\\tilde{O}=${\n$o_i \\ | \\ cost(O_i^*,C) \\ge \\delta ALG/k$} and $\\tilde{C}=$ $C\\backslash$ {$C[o_i] \\ | \\ o_i \\in O^* \\backslash \\tilde{O}$}. Then, in the construction of ideal multi-swaps, each swap is consist of $In$ and $Out$ such that $In \\subseteq \\tilde{O}$ and $Out \\subseteq \\tilde{C}$. If my understanding is not wrong, some of the optimal centers in $\\tilde{O}$ is not used for construction of ideal multi-swaps. If not all optimal centers is used at least once for constructing ideal multi-swaps, then the bound for summation of clustering cost for all \"bad'\" optimal clusters may not hold in Lemma 11. The ideal multi-swaps construction considers the two caes: (1) for each center $c_i \\in \\tilde{C}$ that is neither busy nor lonely, form $In$ with the optimal centers captured by $c_i$ (multi-swap) and borrow $|In|-1$ lonely centers from $L$ (2) for each center in $\\tilde{C}$ that is busy, form single swap using each of optimal centers captured by $c_i$ and the lonely centers in $L$. However, there may exist some case that for some optimal clusering center $o_i^* \\in \\tilde{O}$, the current center $c_j$ that captures $o_i^*$ is not in $\\tilde{C}$ and hence $o_i^*$ will no longer be used for constructing ideal multi-swaps. Here is an example where $o_i^*, o_2^*, o_3^*$ is the optimal clustering centers, and $c_1$, $c_2$, $c_3$ is the current centers with swap size $p=2$. In this instances, it holds that $c_1$ captures $o_1^*$ and $o_2^*$, $c_3$ captures $o_3^*$, $cost(O_1^*,C) \\ge \\delta ALG/k$, $cost(O_2^*,C)<\\delta ALG/k$ and $cost(O_3^*,C) \\ge \\delta ALG/k$. Hence, we have $\\tilde{O}=${$o_1^*,o_3^*$} and  $\\tilde{C} =${$c_2,c_3$}. During the construction of ideal multi-swaps, since $c_1$ is not in $\\tilde{C}$, $o_1^*$ will never be used during the construction. \n\n3. Another main weakness of this paper is the experiments. (1) As a local search algorithm, it’s important to validate the performance of the proposed algorithm within a certain time limitation, which is not presented in the experimental results. This paper mainly compares the experimental performances with fixed local search steps. It is unfair since multi-swap local search method takes more time than single-swap local search methods (i.e., the LS++ method). (2) In experiments, the maximum local search step is set to be 50, which is much smaller than the theoretical bounds (100000kloglogk) for obtaining good performance for LS++ algorithm. For this setting, even single-swap local search method may not converge to a good local optimal solution. The authors should present more evaluation for larger local search steps (i.e., several hundreds or thousands of local search steps) and fixed running time that is large enough for each algorithm to reach the convergence. (3) The size of the tested datasets is rather small (i.e., no larger than 500,000). However, in the recent results for clustering algorithms, instances of size over 1 million [1] or even 100 million [2] have been considered. Since multi-swap local search is much slower than the single-swap local search method (even with heuristic acceleration by avoiding the enumeration of all subsets of data points for swapping in), it is unclear that whether the proposed algorithm can scale well on large-scale datasets compared with single-swap local search method.\n\n[1] Ren J, Hua K, Cao Y. Global Optimal K-Medoids Clustering of One Million Samples[J]. Advances in Neural Information Processing Systems, 2022, 35: 982-994.\n[2] Matsui Y, Ogaki K, Yamasaki T, et al. Pqk-means: Billion-scale clustering for product-quantized codes[C]//Proceedings of the 25th ACM international conference on Multimedia. 2017: 1725-1733.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their time and effort put into reviewing the paper. We start by addressing the weaknesses suggested by the reviewer and then we answer the specific questions.\n\nAnswers to weaknesses:\n\nWeakness 1. Although the algorithm is an extension of the one in Lattanzi and Sohler, the analysis as carried out in [1] does not yield any low constant. We incorporated more elements into the analysis (some of which are borrowed from [1]) to obtain a better constant. We note that in the study of approximation algorithms, once the first milestone of a constant factor approximation is achieved, the goal becomes to minimize the gap between the lower and upper bound. As an example, a large body of work is devoted to improving the constant approximation factor for popular clustering objectives like $k$-median and $k$-means [3, 4, 5].\nOverall, we agree that the techniques used in the proofs in the main body of this paper are similar to those in [2]. However, we believe that (i) improving the constant for an interesting constant-factor approximation algorithm is per se interesting; (ii) the main technical novelty is in the 9+eps algorithm, which is presented extensively in the appendix.\n\nWeakness 2. \nWe thank the reviewer for finding this minor issue in our proof. We submit in the rebuttal an updated version of the paper (we couldn't attach it to our response to the specific reviewer) where we implement a simple change. In particular, adding lines 192-194 and 218 to our manuscript fixes this issue, which are highlighted in blue to make it easily verifiable. The reviewer rightfully pointed out that not all optimal centers in $\\tilde O$ were assigned to a swap. Indeed, this happens whenever $o \\in \\tilde O$ and $C[o] \\not\\in \\tilde C$. We fixed this by creating 1-to-1 swaps between these “orphaned” centers and the leftover lonely centers. We have many enough of such lonely centers because by $|\\tilde C| \\geq |\\tilde O|$. \n\nWeakness 3.\n(1) The reviewer raises a good point. We omitted this experiment given that the running times only differ by a small constant and the times reported suggest that our algorithms still outperform the 1-swap algorithm when run against a deadline. We will add an experiment where we compare the different algorithms against a set of deadlines that is set in the range [1x, …, 20x] the time for a single iteration of Lloyd's algorithm. In our revised version attached in the rebuttal we include this experiment for two of our datasets (KDD-BIO and KDD-PHY) in Figure 3. The experiment confirms our observation that our algorithm performing swaps of size greater than one still achieves better performance compared to the single-swap version of the algorithm despite performing fewer local search steps.\n\n(2) We note that this is a standard setting as in the case of the single-swap local search paper by Lattanzi and Sohler. Even for 50 iterations the algorithms appear to reach a point of diminishing returns in terms of quality improvements. \n\n(3) The main point of our work is to obtain a practical algorithm with good theoretical guarantees. We didn't try to optimize our implementation and we didn't employ any parallelization, which would lead to faster implementations. We believe the current setting of our experiment aligns well with the goal of our paper. \n\nWe further note that the scalability comparison to the two papers mentioned by reviewer is not valid:\nThe paper by Ren et al. only obtains reasonable times for datasets with millions of points for the case of $k=3$ and by applying parallelization using over a thousand cores, and still these times are in the order of hours.\nThe paper by Matsui et al. considers only a heuristic approach, and shows that it is over an order of magnitude faster than plain $k$-means++. We iterate that our goal is to devise an efficient algorithm with provable approximation guarantees.\n\nFinally, we recall that even for the case of the 10-swap version of our algorithm, each iteration is faster than each iteration of Lloyd which is the standard post-processing algorithm applied in practice.\n\nAnswers to questions:\n\n(1) See comment to Weakness 2 above.\n\n(2) Notice that the success probability indeed depends on $\\delta$, as indicated by the notation $\\Omega_{\\delta}(1)$, which is introduced in Line 108. This notation is used to hide the depende on $\\delta$ to hide unnecessary complexity in the claims, and is assumed to be a constant as stated in Theorem 3.\n\n(3) We agree that a simple formula to determine the approximation factor would be ideal. We note that solving equation (1) for $\\eta$ would derive $\\eta = \\sqrt{5 p^2 + 4 p + 1}/p + 1/p + 1$, and hence the approximation factor in terms of $p$ and $\\delta$ becomes $\\left(  \\sqrt{5 p^2 + 4 p + 1}/p + 1/p + 1 \\right)^2 + \\delta$, which is arguably not a very helpful formula to state in our main theorems. For this reason we preferred the simplicity of stating the theorems in terms of $\\eta$, and just state that $p$ should be large enough. In the next version of our paper we can explicitly state the solution to equation (1) in a form of a footnote to facilitate the reader.\n\n(4) See comment on weakness 3 (2).\n\n(5) See comment on weakness 3 (1).\n\n[1] Kanungo T, Mount D M, Netanyahu N S, et al. A local search approximation algorithm for k-means clustering[C]//Proceedings of the eighteenth annual symposium on Computational geometry. 2002: 10-18.\n\n[2] Lattanzi, Silvio, and Christian Sohler. \"A better k-means++ algorithm via local search.\" ICML 2019.\n\n[3] Jain, Kamal, and Vijay V. Vazirani. \"Approximation algorithms for metric facility location and k-median problems using the primal-dual schema and Lagrangian relaxation.\" JACM 2001.\n\n[4] Ahmadian, Sara, et al. \"Better guarantees for k-means and euclidean k-median by primal-dual algorithms.\" SICOMP 2019. FOCS17.\n\n[5] Cohen-Addad, Vincent, et al. \"Improved approximations for Euclidean k-means and k-median, via nested quasi-independent sets\". STOC 2022."}, {"review_id": "2Ilc3KhNkV", "reviewer": "Reviewer_coFq", "summary": "This paper studies local-search algorithms for k-means clustering. The goal here is to obtain a local-search algorithm which (1) give a close to 9-approximation ratio (which is best possible for local search algorithms) and (2) is practical. In the past literature, there has been many local search algorithms developed for the k-means problem, but all of them either have a relatively big approximation ratio or have a prohibitive running time in practice due to the complexity of the local moves at each step. \n\nIn this paper, the authors give a simple local search algorithm which is practical and guarantees a 10.48-approximation. They also give a more complex local search which matches the approximation lower bound of 9, but with a significantly faster running time than in Kanungo et al. [2004]. However, this algorithm is significantly more complex and is unlikely to be used in practice. They also give experiments that show that their first algorithm and a slight variant of it perform well in practice.\n", "questions": "It is unclear to me if the approximation ratio of 10.48 is tight for the algorithm. As well as the 26.64 for the single swap version.\n\nIn experiment section, I do not see any discussion on the running time of the compared algorithms? Maybe it would be nice just to mention it.", "limitations": "Yes.", "rating": 8, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 4, "strengths": "In my opinion, the strengths are as follows.\n1) The paper is clearly written.\n2) The first algorithm is a nice result since it combines practicality with an almost best possible theoretical guarantee. As far as I know (not being an expert), this is the first time the D^2-sampling is used and analyzed in combination with multi-swap local search. Overall, the techniques in the paper are nice.", "weaknesses": "I do not see an obvious weakness in this work, the results are tight or close. Overall, I think this is a nice paper that fits well in the context of NeurIPS. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for their time and effort put into reviewing the paper.\n\nAnswer to the questions:\n\nWe did not manage to prove the tightness of our analysis, neither for $p = 1$ (26.64..), nor for $p$ approaching infinity (10.48..). This is an interesting open question. However, it seems to us that this should be the case. Indeed, if we could consider the optimal centers as candidate centers to swap in, then the term $2 + \\delta$ in Eq. (2) would turn into a $1 + \\delta$. If that term is $1 + \\delta$, then the analysis gives a $9 + \\epsilon$ approximation ratio for p sufficiently large. This suggests that, allowing only data points is actually the only factor holding our algorithm back from achieving optimal approximation ratio. Moreover, the factor $2 + \\delta$ in Eq (2) is best possible: think of a cluster where all data points lie on a sphere centered in their optimal center. Thus, we should expect the analysis to be tight, at least for $p$ approaching infinity. All this is a bit speculative, and we have no formal proof that our analysis is actually tight.\n\n\nWe kindly note that the running times of the algorithms (per iteration) are reported in the legends of Figure 2 following the algorithm names. We realized this was confusing, and we will report the running times more prominently in the next version of our paper."}]

Author rebuttal:
While we respond to each reviewer individually, we attach here a revised version of our paper where we incorporate things discussed with reviewer 2sgw.

We note that while this attached version is slightly longer than 9 pages, we only use it to refer to it in this rebuttal phase, and will be made to abide to all formatting guidelines in the final version.

Meta-review: This paper is a difficult case as one very thoughtful reviewer has raised valid concerns that seem to be largely addressed, though the authors did not follow the correct format in uploading a rebuttal version of the changes. These concerns are in a sense outweighed by the positive consensus among other reviewers, who relatively provide very little information behind their assessments.
Decision: Accept (poster)