Linear Time Algorithms for k-means with Multi-Swap Local Search
Authors: Junyu Huang
Qilong Feng
Ziyun Huang
Jinhui Xu
Jianxin Wang
Venue: NeurIPS 2023 poster (2023)
Keywords: Approximation Algorithms
k-means Clustering
Local Search

Abstract: The local search methods have been widely used to solve the clustering problems. In practice, local search algorithms for clustering problems mainly adapt the single-swap strategy, which enables them to handle large-scale datasets and achieve linear running time in the data size. However, compared with multi-swap local search algorithms, there is a considerable gap on the approximation ratios of the single-swap local search algorithms. Although the current multi-swap local search algorithms provide small constant approximation, the proposed algorithms tend to have large polynomial running time, which cannot be used to handle large-scale datasets. In this paper, we propose a multi-swap local search algorithm for the $k$-means problem with linear running time in the data size. Given a swap size $t$, our proposed algorithm can achieve a $(50(1+\frac{1}{t})+\epsilon)$-approximation, which improves the current best result 509 (ICML 2019) with linear running time in the data size. Our proposed method, compared with previous multi-swap local search algorithms, is the first one to achieve linear running time in the data size. To obtain a more practical algorithm for the problem with better clustering quality and running time, we propose a sampling-based method which accelerates the process of clustering cost update during swaps. Besides, a recombination mechanism is proposed to find potentially better solutions. Empirical experiments show that our proposed algorithms achieve better performances compared with branch and bound solver (NeurIPS 2022) and other existing state-of-the-art local search algorithms on both small and large datasets.

Linear Time Algorithms for k-mmans with Multi-Swap Local Search
Junyu Huang 1.2 Qilong Fengl.2. Ziyun Huang", Jinhui Xu Jianxin Wang 12,5. School of Computer Science and Engineering. Central South University, Changsha 410083, China 2Xiangjiang Laboratory, Changsha 410205, China
Department of Computer Science and Software Engineering, Penn State Erie, The Behrend College
*Department of Computer Science and Engineering, State University of New York Buffalo. NY, USA
The Hunan Provincial Key Lab Bioinformatics, Central South University, Changsha 410083, China
j junyuhuangceu@@ooinaa com, csufengemail, csu edu en, zxh201 @psuu. edu, iinuiinuufa falo. edu, jxwang@mail, csu. edu. cn
Abstract
The local search methods have been widely used to solve the clustering problems. In practice, local search algorithms for clustering problems mainly adapt the single- swap strategy, which enables them to handle large-scale datasets and achieve linear running time in the data size. However, compared with multi-swap local search algorithms, there is considerable gap on the approximation ratios of the single-swap local search algorithms. Although the current multi-swap local search algorithms provide small constant approximation, the proposed algorithms tend to have large polynomial rumning time, which cannot be used to handle large-scale datasets. In this paper, we propose multi-swap local search algorithm for the k-means problem with linear running time in the data size. Given a swap size t, our proposed algorithm the current best result 509 (ICML 2019) with linear running time the data size. Our proposed method, compared with previous multi-swap local search algorithms, is the first one to achieve linear running time the data size. To obtain a more practical algorithm for the problem with better clustering quality and running time, we propose a sampling-based method which accelerates the process of clustering update during swaps. Besides, recombination mechanism is proposed find potentially better solutions. Empirical experiments show that our proposed algorithms achieve better performances compared with branch and bound solver (NeurIPS 2022) and other existing state-of- state-of-the-ar the--ar local search algorithms on both small and large datasets.
achieve (50(1 + !) + >)-aprroximation which improves
Introduction
Clustering is fundamental problem in the field of machine learning with many real-world applica-
tions. The goal of clustering partition a given set of data points into different clusters according to their similarity such that data points within the same cluster share high similarity as much as possible. Among different objective functions, the k-means clustering aims to minimize the sum of the squared distances between data points to their closest centers. More formally, given a set P C Bd
"Correpooning Authors
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

of data points in d-dimensional Euclidean space, the goal of the k-means clustering is to find set C' C IR of size at most l with the following objective: minCCE* Eperminnec: min llp ap--cl cll².
For the k-means problem, Lloyd's algorithm [2] is one of the most widely used heuristic in practice. However, there is no theoretical guarantee for Lloyd-type method unless certain data distribution assumptions are introduced. is known that there several constant approximation schemes based The current best approximation ratio in polynomial time is 5.912 [3] which is based on primal-dual method and nested quasi-independent set. For fixed dimension or the number of clusters k; several (1 + d-aaproximmaioon algorithms
on primal-dual and randomized rounding techniques [ g
were proposed 6 7
The k-mcans++ algorithm proposed by Arthur and Vassilvitskii [3 is a good seeding method that runs in linear time in the data size with Olog k) -aapproximation. It is also known that the +-means++ algorithm gives a constant approximation by opening O(k) centers I 13 16. Lattanzi and Sohler [1I showed a combination of k-means++ seeding and the local search algorithm (named as LS++ algorithm), which yields a constant approximation with running time Ondkk² log log k). In each round, LS++ algorithm samples data point using k-means++ seeding and enumerates possible swap pairs to make improvements on clustering cost. They proved that after Olk log log k) rounds of sampling and swaps, one can obtain -0p-appro solution in expectation. Choo al. 4 proved that one can achieve an O(1/e* >apprrximatiin using reduced ((ek) rounds of LS++ algorithm. Under the assumption that each optimal cluster has size (2(n/k), Huang et al. gave an improved approximation algorithm with ratio (100 + e) by random sampling methods.
However, there are still several issues for local search methods. Although current local search methods with multi-swap strategy can achieve good theoretical guarantee, the running time of them have polynomial dependence
the data size, which are hard to be used to handle large-scale datasets. Compared with the (9 + e >approximmation multi-swap local search algorithm given in [10), the approximation ratio 509 of LS++ is a large constant since it can only apply the single-swap strategy. Numerical experiments [9] showed that LS++ could easily fall into a poor local optimum when handling real-world datasets. An immediate idea is apply the multi-swap strategy to LS++ algorithm for improvements, However, the swapping process of LS++ relies heavily on the one-to-one matched swap pairs defined in (II). Thus, it is challenging apply the multi-swap strategy solve the k-means problem while maintaining linear dependence of data size on the running time.
Secondly, in the process of local search swaps, the time and space complexities for clustering cost update during swaps have great impact on the efficiency of the algorithms. A direct way in [ is to maintain the nearest and the second nearest centers for each data point such that picking the best swap pair and updating the clustering cost can be implemented in time (mdd) and (ndlk), respectively. Maintaining the distances from data points their centers requires an (nd)). To obtain faster implementation, as pointed out in 4 one can use binary search trees to store the distances from each data point to each of the clustering centers. By using this data structure, each local search step can be implemented in time O(nd log k). However, the space complexity becomes O(ndk). To get much practical algorithms for the k-means problem, it is necessary to further improve the time and space complexities for updating the clustering cost during swaps.
space complexity of
1.1 Our Contributions
In order to further narrow the gap between theory and practice, in this paper, we propose the first multi-swap local search algorithm for the k-means problem with linear running time the data size. A common feature for the existing multi-swap local search methods is that (((n))) candidate swaps should be enumerated for finding clustering cost improvements in a single local search step, which leads to at least quadratic running time. To overcome this challenge, our idea is to use sampling-base based strategy to construct a candidate set of centers that are close to the optimal clustering centers in linear time, which serves as the set of potentially good centers for swapping in. Based on the candidate set of centers constructed, enumerations on the current set of centers opened suffice to determine good swap. Hence, the number of possible swaps can be reduced from O((nk))) to @(k'), where each local search step can be conducted in linear-time in the data size.
Sampling- strategy has been used in |I for designing single-swap local search algorithms. However, the theoretical analysis relies heavily on the one-to-one matched swap pairs. Thus, the approximation ratio is large constant since there may exist some optimal clusters that cannot be well

---

approximated by performing matched swaps. In this paper, to obtain better approximation guarantee, we extend the notions of swap pairs and propose new consecutive sampling method to construct candidate centers for swap such that data points close to a subset of optimal clustering centers can be swapped in simultaneously. A key challenge here that there may exist some optimal clusters whose clustering costs only take tiny fraction of the total clustering cost such that sampling methods may fail. To overcome this challenge, we propose new structures that divide optimal clusters into different groups for establishing lower bound for the success probability of sampling.
By using the proposed multi-swap local search method (denoted as MLS algorithm for short), given a swap size f, an improved (50(1 + 3) + e)-approximation can be obtained in time O(ddk.2++1 log(c-" log k)). To benefit more from the proposed multi-swap local search method when handling large-scale datasets, we propose a sampling-basee gorithm. The proposed algorithm (denoted as MLSP algorithm) accelerates the updating of clustering cost during the swaps such that each local search iteration in MLS algorithm can be implemented in time O + poly(k)d (k)) with extra space complexity of Ö(d)). In order to obtaín better clustering quality, we develop recombination mechanism in MLSP which combines sampling and scoring methods to help the local search algorithm find better solutions when the search falls into a poor local optimum. By picking the top-k data points with the highest scores as new initialization, the chance to get out of the local optimal solutions becomes large. Numerical experiments show that our proposed method achieves better performances compared with branch and bound solver and other local search algorithms. The main contributions of this paper are as follows.
method for accelerating the ML.S al-
We propose the first multi-swap local search algorithm (MLS algorithm) with running time linearly dependent on the data size. Given a swap size t with 2 2. our MLS algorithm achieves a (50(1 + )) apprroxxiintio in time O(ndk22++1 log(c log k)), which improves the current best approximation ratio 509 with linear running time in the data size.
We give practical heuristic algorithm (MLSP algorithm) for better implementation of the proposed MLS algorithm, which accelerates the process of clustering cost update during swaps and provides better scalability to large-scale datasets. Besides, recombination mechanism proposed to prevent the local search algorithm falling into poor local optimum too early.
Table 1 shows a detailed comparison of our results with the state-of ofite--f-theaarr ones.
Result
Approximation Guarantee
Method
Assumption
Running Time
Ollog k)
k-meuns++
O(nak)
[3 2)
Multi-Swap Local Search
O ig d'd A)
509
Sampling Single-Swap Local Search
O OludkPioglogk)
((1)
Sampling Single-Swap Local Search
O(ndk log k).
100 +
Sampling + Single-Swap Local Search
PX >
O(ndka log -
This Paper
50(1 + !) +
Sampling + Multi Multi-Swap Local Search
O(ddi:2 log(e- log k:)
Table 1: Comparison with related results on k-means clustering, where the size of the given dataset, d is the dimension, k is the number of clusters opened, t is the parameter representing the swap size of local search methods, and A the aspect ratio (aspect ratio is defined as the maximum pairwise distance of the given instance divided by the minimum pairwise distance).
2 Preliminaries
We use P C IRdd and l to denote the given dataset and the number of clusters, respectively. For any two points p. q E P. we use d(p, 9) = 1p ql12 to denote the squared distance between them. Given two sets A, B C P, let (AA, B)= Epe A minge d(p, q) denote the sum of the squared distances from data points in A to their closest points in B. Let C* = {ct. cg. ci.,... } be an optimal solution. We use P(C = (Pr, P Pi..... PH} to denote the corresponding optimal clusters by assigning data points in P to their closest centers in C*. Let Opt be the cost of the optimal solution. For subset Q S PCCC* of optimal clusters, we use (QQ) = UP: EQ P; to denote the set of data points in clusters of Q. Denote Z'(Q)= = {ci= : P% E Q] as the set of optimal centers of clusters in Q. Given a subset S C* let J(S)= = [Pr of E S) be the collection of optimal clusters whose clustering centers
are G. For an integer I, let [1] = 21,2,.. .,.,, 1). The following lemma is a folklore for the k- k-means problem.
3

---

Lemma 1
3 Let P S R be a set of data points, and ju() = EpEp denote the center of gravity. For any data point c E Rd, we can get A(P, (c)) = Pld(u(P), + A(P, {μ(P))).
Theorem 3
Algorithm returns an (log -approxiiaeee solution time (ndk)).
Algorithm 1 k-means++
Input: An instance (P, k) of the k-means problem.
Output: A set C _dd of centers with size at most k.
1: Randomly sample point P and set C = {p)-
2: for i= 1 to k 1 do
3: Pick point p E P with probability ((p)) ()A(P.C)) and add p to C.
4: return C.
3 Linear Time Local Search Algorithm with Multi-Swap Strategy
The general idea solving the k-means problem with multi-swap local search is that we propose a new consecutive sampling method to construct candidate centers for swap such that data points close to subset of optimal clustering centers can be swapped in simultaneously. The multi-swap local search algorithm (denoted as MLS) is given in Algorithm 2. There are mainly two stages in each round of the MLS algorithm. Given a swap size i, in the first stage (steps 3 to 5), a candidate set of centers with size will be sampled using the k-means++ + method. This avoids enumerating all the data points for constructing the candidate sets for swapping in. In the second stage (steps 6 to 7), the algorithm enumerates all subsets (with size most t) of the candidate set of centers opened for swapping out. By extending the notions of swap pairs to swap set and carefully analyzing the structures of local optimal solutions, we prove that the clustering cost can be reduced significantly with certain probability in each iteration of the MLS Algorithm. The following is the main result of this paper.
Theorem In the i-th iteration of Algorithm 2. let C" be the set of centers obtained in step 7. If the current clustering cost A(P, C) is larger than 50(1 + O)tt,, then with probability at least ((k-"), we have A(P,C") S (1 TOok (PP, C). After ((k((t) log(e" lo) iterations, we ger an approximate solution with ratio (50(1 + 7) + E) in expectation
3.1 Analysis
In this subsection, we analyze our proposed Algorithm 2, where a candidate set of centers for swap is constructed by independent sampling steps in each iteration. Our objective is to show that the clustering cost can be reduced with certain probability in each iteration. Due to space limit, all the detailed proof' are given in Appendix A.
In the following, we consider a single iteration of the proposed MLS algorithm. Assume A(P, C) > 50(1 Opt holds within single iteration in Algorithm 2. Otherwise, C is already a 50(1 ))- approximate solution for P. Let C = {c1, C2 .... Ck denote the set of centers before the swap (steps 6-7) of Algorithm 2. Let P(C) {P1, P2.... P&) be the corresponding partition of clusters induced by C. For an optimal cluster Pž let oj be its clustering center. For each cluster PA E IP(C), let Ch C denote its clustering center. Following the ones in [1], extend the definition of good clusters with respect to C" as follows.
Definition 1 Good single cluster. A cluster P P(CC") is called good with pair of points (ch, cj) such that cj & C and A(Pa, C) - C(P.C, Cr. - (cth) Tho A(P, C), where <(P, C, = A(P', C. Clfc,)) AP C)) is the reassignment cost by swapping out. Otherwise we say that Pr bad single cluster with (c),
Definition 2 Good I-Clusters. Given an integer with I. 2 2, for a collection of optimal clusters Q C (CC')) with [Q] < Q is called good with a pair of sets (Z' (Z))), V) such that V C C. IV = (2)) and (ZZ(Q), C*) - ((P, C, ('(Q),
(Q) C*) Took A(P, C),, where = A(P) C. ((P) C) is the reassignment cost by swapping the points in V out. Otherwise say that Q is a of bad t-clusters with (Z'(Q), V).
C(P.C,Z'(Q))
Note that the authors did not try optimize the constants

---

Algorithm 2 MLS
Input: An instance (P, k) of the k-means problem, parameters T and t. Output: A set C _ Rd of centers with size at most k.
1: Initialize C = k-meanst+(P, k).
2: for i = 1 to T do
3: 1= 0.
4: for j j= 1-110 to t do
5:
Pick a point p E with probability A(fp), C))((PP and add to
6:
if U C I and V _ C s.t. [U] = [V] and A(P, (VP,)) < (1- 100%
1
) C*)
then
7:
C = C.V
8:
return C.
The above definitions estimate the changes of clustering cost by replacing a set of centers V C with a set Q' of data points that are close to a subset of optimal clustering centers, where new clustering is constructed by reassigning data points in Z (Q) to Q' and data points in (IZ(Q) to C.V. The main idea behind Algorithm 2 is to iteratively find data points close enough to the optimal clustering centers for swap to make reduction on clustering cost by at least (1 - + In the following, we will prove that with probability at least (2k-"), the sampling and swap process induces a significant reduction on clustering cost in each iteration.
C
We start by dividing the optimal clusters into several groups to give an upper bound of reassignment costs. Given a swap size f, for each optimal center ci, E C we define (cc,, as a mapping function that maps of to its closest center in C. For simplicity, we say that cj is captured by V(cj,). For a center Cj E C, let (cg) be the set of optimal centers captured by cj. If (c)) = 0, then oj is called a lonely center. If [ -w'(c,,) [c,,]] = 1, let o7 = (c,). Then it is called that (c), eÿ) forms type-1 matched swap pair. If 1 < (e,), S t, let σn be arbitrary set of unused lonely centers with 0nl | (c,)l 1. Let A = v (c3)) and A' = {cg} Uog. Then, it is called that (A, A') forms a type-1 matched swap set. Let M1 = { ea) : (c), is a type-1 matched swap pair] and M2 = {(A, AA) (A, A') is type-1 matched swap set) be the collections of type-1 matched swap pairs and type- I matched swap sets, respectively. If v "c,)) > t, then find each lonely center Cy E C that has not been used for constructing type-1 matched swap pairs sets. For each lonely center Cy and each cf E (c)), (ch, G2y forms type-2 matched swap pair. We use M3- = [(c), (g)= : (c%, (g) is a type-2 matched swap pair) to denote the set of type-2 matched swap pairs. For a set V C of centers, let X(V) Ph be the set of data points in P whose closest centers are in V. Given a subset S C C' of clustering centers, we also use J(S) = {Ph ca E G) to denote the set of clusters whose centers are in S. The following lemma gives upper bounds of reassigmment cost by matched swap pairs or sets.
Lemma 2 Given a rype-1 or type-2 matched swap pair (ch, cj). it holds that C(P, C, ch. 9) S 24 C*) + AA(, C). Given a type-l matched swap set (Q, V), it holds that 24 (AA((()) + +
Let H4 = fch: cg) E 3 be the set of optimal centers that participate in constructing type-1 matched swap pair. Let Hz, = {A (A, A') M2} be the collection of the subsets of optimal centers that participate in constructing type-1 matched swap set. Let L = {c, : (c,, c2) E M2} be the set of optimal centers that participate in constructing type-2 matched swap pair. Let H3 = {cy (c), Cy) E } be the set of centers in C' that participate in constructing type-I matched swap pair. Let H2 = {A' (A, A') E .22 be the collection of the subsets of centers in C that participate in constructing type-1 matched swap set. Let L' = {c₂ (ch, Cq) E M2} denote the set of centers in C that participate in constructing 2 matched swap pair.
During the sampling and swap process in steps 4-7 of Algorithm 2, there are two cases that may happen: (1) 3 cj E L. such that PV is a good single cluster with type-2 matched swap pair (ch, cg) E M3: (2) oj E L. P is a bad single cluster with any type-2 matched swap pair (ch, C.) E M3. We will discuss the two cases separately in the following. If case (1) happens, we first show that with probability at least 2 the clustering cost can be reduced at least by 1 - O Let of E T. be an optimal center such that P is good single cluster with a type-2 matched swap pair (ch, cy) M3. By the definition of good single cluster, we have A(P,, C²) 2 9A(P,, (c]]).
5

---

as the set of data points in PK that are Otherwise, the clustering cost of data points P (P)) is at least A(P) (g)) using c/ as center, which contradicts with APPAXx(PM) (G)) < A(P, (c.))-
Define x(P) E E Pk d(p, c)) S 1.A(P, PR close to the optimal center ch. Observe that (x(P 2
3 I
Next, we will show that whenever an optimal cluster Px has large clustering cost with respect C, i.c. A(P C)= = b {c(,) for a real number b 2 3, it suffices to use data points in x(P) to approximate the clustering cost of Ph.
Lemma 3 Let P7 be an optimal cluster with A C) 64( {c))) for a real number b 2 3. Then, A (x 2 2uo 1) (G.))
We now argue that the clustering cost of each good single cluster takes certain fraction of the total clustering cost. Then, by sampling according to the squared distances, with good probability, data points close to the center of good single cluster can be sampled. According to the definition
of good single cluster, we can assume that A(P; C) = 6A(P% {c))) for a real number b 2 9. b-1 A(P. (ct,)) 2006
By Lemma 3, we know that A(x(P)), 2 200
6-1 A(Pi, C) 2 300
1 A(P)
By the definition of good single cluster, we also have A(P C) Z TOOE A(P, C), which implies that Axx(P)), Z S0000R A(P, Thus, in each round of the sampling process in steps 4-5 of Algorithm 2, with probability at least (2))), we can sample a point q E (4%) for good single cluster P" with type-2 matched swap pair (cf, Cq) E M3. The following lemma shows that, by swapping 9 with Cy and assigning all the data points in P to q, the clustering cost can be reduced at least by 1 1 -()).
Lemma 4 By swapping with Cap the clustering cost of (P,, C*) can be reduced at least by 1 (
We have shown that if case (1) happens, we can sample a data point g close to the center of good single cluster P7 to make the clustering cost reduced significantly. Next, we assume that case (I) never happens and case (2) happens. In case (2), the idea behind is to sample data points close to the optimal centers in H, or sets of optimal centers in H2 to reduce the clustering cost. We first bound the clustering cost of optimal clusters in (L).
Lemma 5 If case (2) happens, then A(X(L), S (1+ ()) + 24A(((L).C') + (A(X C) + 100 A(P,
Now, consider set of centers Q H2. We will divide the optimal clusters in ((Q) into two groups. Let Qr = {P% E J(Q): A(Pr, ((P,C)) > sonwooga A (P,, C}} be the set of optimal clusters in J()) with large clustering cost with respect
and Qs = {Px E J(Q) : A(P C) < 30000 A(P.C) ])} be the set of optimal clusters in (() with small clustering cost with respect to C. respectively. We define Q's [PP E Qs A(PR C) < 3A { as the set of optimal clusters in Qs whose clustering centers are close to one of the centers in C. Let Q's = Rs'O's and Qr Q2 U Q²2 respectively. We first show that, it suffices to only consider optimal clusters in Qr-
Lemma 6 Let Q E H2 be set of centers in H2, where (Q) is a set of good t-clusters with type-1 matched swap set (Q, A') E M2. Define V = A';h where oj is the center in A' with (c,) 1. Let U S P be the set of data points with IU [V] such that U n (x()))) + 0 holds for each PA E Qr. Then, A(P.CVV S (1
1 A(P, C).
Took
Note that there are t sampling iterations in each step of Algorithm 2. Let H₂ = (P% : PA E QT, Q E H2} and H, = ((H)) U Hj. We will define mapping function m which maps each Pr E H, to an integer m(P;) E [t]. For each P E H, such that of E H,, we define m (P:) = 1. Then, consider each Q E H2. For each Pr E Qr, we define m (PP)) i such that i E [[QT]] and m(P:) # m(()) for any two optimal clusters PH, P3 E Qr- For each Q E H2, let (Q) be the success probability that for each PA E Qr, a data point is sampled from x(P)) in the m(p)) -th iteration of Algorithm 2. For each o E 11,, let P. (cf) be the success probability that a data point is sampled from X (PP)) ) in the first iteration (note that in this case m P7) 1) of step 4. Let H² (c), E H, P is a good single cluster with the type-1 matched swap pair (ct, 03) E Mj } be the set of centers in H1 whose optimal clusters are good single clusters with type-1 matched swap pairs in M.. Let H2 = [Q E H2 J(Q) is a set of good 1-clusters with the type-1 I matched swap set (Q, V) E M2} ) be the collection of the sets of centers in H2 whose corresponding optimal clusters are good

---

t-clusters with type-1 I matched swap sets in M2. Define = E Po(c)) + EgEnz (Q) as the summation of the success probability. Since all the optimal clusters in (L) belong to bad single cluster, and the events defined related to Ps or Pa(c)) for each Q E H2 and cj E H1 are mutually exclusive, p gives a lower bound success probability to get a (1 1000 reduction on the clustering cost in each iteration of steps 2-7. In the following, we will show how to obtain lower bound for p!.
Consider a set Q E H2 of optimal centers such that JQ) is a set of good I-clusters with a type-1 matched swap set (Q, A'). There are two subcases that may happen: (1) PA E J()) such that Pk is a good single cluster with a swap pair (c), Gj). where ey a lonely center used for constructing type-1 matched swap set; (2) V P E J(Q). P% is a bad single cluster with any swap pair (cf, e)), where Cj is a lonely center used to construct type-1 matched swap set. In subcase (1), since is usually constant and could be much smaller than k with probability at least S2(k= a data point E x(P) can be sampled in the first iteration of step 4 in Algorithm 2 for swap to make the clustering cost reduced at least by (1 T according to Lemma 4. Next, we assume that subcase (1) never happens and subcase (2) happens. In subcase (2), our objective is to sample a set V of data points such that V contains at least one point from x x(P)) for each P% E Qr. For each optimal cluster PV E QL. with probability at least ((k--),
)) we can sample a data point E x(P)) in the m( --th iteration of step 4 in Algorithm 2. Thus, the probability can be bounded by 1211-11.).. Then, only need to consider optimal clusters in Qs- By Lemma it suffices to consider optimal clusters in Q²z. The following lemma gives an upper bound of the failure probability of not sampling a data point from (4%) for each P7 E Q² in the m m(p)) )-th iteration in step of Algorithm 2.
Lemma? 7 Given a set Q E H2 of centers such that J()) is a set of good t-clusters, the probability that step 4 of Algorithm fails to sample a data point q from x(P)) for each PT. E Q3 in the for
m (P)) )-th iteration is at most (1 + 30000k )c -(Z(())) (1008))33333/P.C))
Now we can bound the success probabilities of (Q()) and (c(c for each Q E H2 and of E H1. We will divide the optimal clusters in J(() into two different groups Hg and Hg. Firstly, we consider the centers in H1. For center ch H1, if Pj is good single cluster with the type-1 matched swap pair (ck, cÿ) E M1,, then add Pž to Hg. Otherwise add P to Hg. For a set Q E H₂ of optimal clustering centers, if [Q) is set of good (-clusters with the type-1 matched swap set (Q, V) E M2, then add each P E Q'g to Ha and each P7 E (Q))Q" to Hn. If J(Q) is a set of bad I-clusters with the type-1 matched swap set (Q, E M2, then add each P% E ((Q) to Hn. The following lemma argues that the summation clustering cost of the optimal clusters in He is large.
Lemma 8 For the optimal clusters in HG, we have A (Z C)
A(P.C).
100
For a good single cluster PV with the type-1 matched swap pair ()) E M, where cj E Hi, define py (cf)
ps (ch) as the probability that the first iteration in step 4 of Algorithm 2 fails to sample a data point q E (PP)). Then, we have pf(C))) S (.(( S - /)..0)/3330//CC by Lemma 3. For a set Q E H& of optimal centers, we have that pf(Q%) is the failure probability of not sampling data point from x(P:.) for each Pr E Q's in the m(PF)-th iteration of step 4 in Algorithm 2. Recall that pd = -Ee (c() + +EQEH3 (Q) is the summation of success probability. Then, we have pf 2 1(k (( Ge Po(ch) + EqEng (QEng)o(O')) (O)))) Let P. = EGGEnE Po(c)) +
Eqeng Pa (Q%). Observe that pod 2 1 - Notenyy (1 - po(())) II
-Po())) IIQE (1 (Q2)). Since
there most good t-clusters, by Lemma 7, we have p!" 2 1 - (1 + 30000t )& S00aTP.CT 2(2(gg) 3000(P.C)
1 1 In(1+ n00000
2 1- 1-e soico soico 2 1 soino 2 60001 where the third inequality
1
follows from Lemma 8 and In(1 + +x) S Then, it holds that r! > =2(k=)pp!' = 2(k-'),, ") which indicates that with probability at least ((k-)), we can sample data points close to a set of good t-clusters good single cluster for swap to make the clustering cost reduced at least by (1 1ook according to Lemma 6. Putting all things together, Theorem 2 can be proved (Detailed proof of Theorem 2 is given in Appendix A).
Running Time Analysis. By Theorem 2, in order to obtain (50(1 + :) + c) approxxiee solution, the iteration rounds for Algorithm should O(k log(e log k)). In each iteration, it takes ((udk) time to update the distances between data points to their closest centers. During the sampling process, I data points are sampled according to the D2-Sampling distribution to serve as the candidate set of centers for swapping in, which takes time ((t)) if the distances from data points to their
7

---

centers are already known. It takes Okk) time to enumerate each subset with size at most of the set of current centers opened. It takes (Indk") time to recalculate the clustering cost after each swap nearest centers of each data point are maintained during the whole process. Thus, the total running time of Algorithm is Ondd log( - logk)). ))).
3.2 Accelerating Multi-Swap Local Search for k-means
In this subsection, we provide more practical algorithm for the k-means problem to accelerate the proposed multi-swap local search process. The algorithm is given in Algorithm 3. The main idea behind is to use sampling-base
methods to obtain fast clustering cost updating during swaps. In step 7 of Algorithm 3, a small sample set S of size log is randomly taken from P. For clustering cost updating. instead of calculating the clustering cost of all the data points in P. we use the clustering cost of S _ P as an estimation. This reduces the time for picking the best swap pair during single local search iteration from (ndd) to (ppoly(k)d). Then, in steps 16-18, we design a recombination method to find better initialization which prevents the local search algorithm from falling into poor local optimum too early. step 16 of Algorithm 3, we randomly take a set D C P centers with size ( log k)) Let C1 C D be the set of the new candidate centers. For each center e₁ C1. we add a score of ((P,. C)/ C) to it. Similarly, for each center C₁ E C, we also add a score APP. C)/ ((P, C) to it. Then, by giving each center in C1 and C' a score weight of 0.75 and 0.25, respectively, we pick the top-k data points with the highest scores as a new initialization of clustering centers find potentially better clustering costs until convergence.
Algorithm 3 MLSP
Input: An instance (P, k) of the k-means problem, parameters T, 6, R', and Output: A set C C Rd of centers with size at most k.
I: Initialize C = k-means++ k), (, = 0, C = 0.
2: while r < R' do
3: for = 1 to T do
4:
-= 0.
5: for j to do
6:
Pick a point p E P with probability ((p)), ; (P/A(P,C), and add p to
7:
Randomly sample a set S from of size log n²
8:
Let (U,V) V) be a swap set such that U I.V C C, |0| = [V] and UU) is
minimized.
9:
if A CVVU) (1 - 100% T (P, C.) then
I
C UV. For each center e E C, find the 50-nearest neighbors in P to C' for improvements on clustering cost by swapping e with one neighbor until convergence.
11:
12:
if (A(P.C) < A(P, Cf) then
C,
14:
else
r + 1, randomly sample set D from with size f log $ and set C₁ CUD.
For each e C. calculate S'(cx) A(P,.C) and add score of 0.25.5 (cr) to Cj-
16:
A(P.C)
17:
For each Ch E C1 calculate S'(cx)) = A(P,.C))
A(P.C) and add a score of 0.755' (c)) to
18: Reset C' as data points in C₁ with top & scores.
19:
return Cj.
4
Experiments
In this section, we compare our proposed algorithms with the branch and bound solver and other local search methods. For hardware, all the experiments are conducted on 72 Intel Xeon Gold 6230 CPUs with 500GB memory.
Datasets We evaluate the performance of our algorithms on 8 datasets used in [5] with sizes over 50,000, two datasets SUSY (5,000,000 x 17) and HIGGS (11,000,000 x 27) from the UCI Machine
8

---

Learning R and one dataset SIFT (100,000,000 X 128) with size 100,000,000 used in 14] In Appendix B, we give experiments on small datasets with sizes smaller than 50, 000 used in 15].
Experimental Setup In our experiments, we choose centers from the datasets such that the problem can be solved by branch and bound solver (denoted as BB for short), which can serve as references of optimal solutions. We compare our MLSP and MLS algorithms with the BB method and other local search algorithms. Following the settings in [5], the centers obtained by different algorithms are projected to their closest data points in P. and the number of clusters k is set to be 3. 5 and 10. Each algorithm is executed for 10 times, and the average results with deviation, the best results and the average running time are given. We test the performance of different choices of parameters in Appendix B. We also conduct the experiments with fixed running time in Appendix B.
Algorithms In our experiments, we consider six algorithms. The first is the BB method in [5]. which is the sfate-of-the-are solver for handling large-scale datasets. For local search algorithms, we compare our MLS and MLSP algorithms with the LS++ algorithm in [1] and FLS algorithm in [9]. using the Lloyd's algorithm [12] as a baseline. For MLS algorithm. we use the sampling-base based method designed for MLSP to accelerate the swapping process. For MLSP and MLS algorithms, we set the sampling rounds T 400 with a swap size 1 2. For MLSP algorithm, the failure upper bound is set to be Rº = 5. For parameters e and n, we fix them as 0.5. For fair comparison, we also set the number of sampling rounds as 400 for LS++. Following the settings in [.] for all local search algorithms, the Lloyd's 's algorithm [2] is used as the final step to adjust the centers.
Motied Dataei
SamgleSee
BBICost)
Bos Clei
Asarge Clos
Tanete)
Delect
SuxpleSior
BBrCned
Bea Cial
Anange Cill
Titelo
5501F4000
PLS
13L.34
136.1922211
21.25
251220 eicc 1.83403
000
132 71
15808-334
0.87
SINIT 1D
C0M
0.63
MLS
131.75
136.718277
3.13
2.50965405
257788 4lec 7.36403
25,39
561610144:11733
$1.5
4.15588447
11 1154l
20.94
163156-34
261856.27 189155
1157.1
1lledd
KBOG
55.41-223
1554E 107
5.15568.10 un
0.38T0E (0::2113300 122388-.00 105
0.21
L000000012
COM
582761991
ST1151 23:37955 96
0.50
MLS5
568582.36
68.11
247377
517.52
F1.
34856-355
120025408
11665--00999 9.dEi63
1iced
1M0000012
24665
2566420066
0.00
USC_1900
2 24407T-44::100000
48.32
MI.5
24680
24599-257
20005500
210636-00--2240
60916
1,5e
1417E-11
357-31
1.3155-
MLSPP
14
1.3058--44:2.022
201.11
0.47
50003001177
3462.57
1.5r
56720
1570
30.75
5834E 100 1.4805
2434.97
1kcad
55546
00505-413 13
0.41
30000331
1 55CAE (00:) 1.2540
171.66
MLSP
175.55
1.87738--88
21V224
LSt+
13
25-10
121353
MI
Llead
SIFT
100000.0000 125
OOM
133500-
596.9
Table 2: Comparison results on clustering costs and running time with k 10 on datasets with sizes larger than 50,000, where OOM is short for out of memory
Results Table 2 shows the results on 11 datasets with sizes over 50,000 using A = 10. For clustering cost, on each dataset, the best clustering cost returned by our MLSP algorithm is smaller than BB method and other local search methods. For each dataset except for RNG RNG_AGR, the average clustering cost returned by our MLSP algorithm outperforms the result of BB method. It can be seen that, as the sizes of datasets grow, BB method requires larger space complexity. When the size of dataset over 360,000, it requires a memory of over 500GB, which is not practical for handling large-scale datasets. For L.S++ method, it is more difficult to find high-quality solutions as the sizes of datasets grow. As for FL.S, although the performance of FLS is better than that of LSS+, our MLSP algorithm improves the performance of FLS on clustering cost with smaller deviation than FLS. By calculating the average values over all datasets, for MLSP algorithm, we can get that the clustering cost is reduced by 1.9%, 1.3% and 2.4% compared with LS++ FLS and Lloyd's 's algorithm, respectively. For MLS algorithm, on average, the clustering cost is reduced by 0.6% and 1.1 % compared with LS++ and Lloyd's algorithm, respectively. As for running time, it can be seen that proposed algorithm scales well as the sizes of datasets grow. By calculating the average values over all datasets, our proposed MLS algorithm is at least 1.79 times faster than LS++ algorithm.
9

---

The experimental results on the performances with varying T and R' (Appendix B.1) show that larger T and R' will not influence the results too much. In general, larger sampling rounds and larger failure upper bound can result in potentially better solutions with higher running time. For parameters and 17, the results (Appendix B.1) show that the performances of our proposed algorithms are almost the same for different choices of and v). Tables 3, 4 and 5 show the comparison results of different algorithms with varying iteration rounds on dataset rds, KEGG and Urban_ 10, respectively. The results show that, our proposed MLSP algorithm always achieves the best clustering cost compared with LS++ and MLS algorithms. Tables 6, 7, and 8 show the results of MLSP algorithm with varying number of failure upper bound R' for fixed с 0.5, 17 0.5 and T = 400. It can be seen that larger failure upper bound will lead to smaller deviation, and the running time becomes higher. Tables 9, 10, and 11 show the results of MLSP algorithm with varying parameters and 17 for fixed T = 400 and R' = 5 on dataset rds, KEGG, and Urban_ 10, respectively. It can be seen that smaller values of e and 1) result in better performances clustering cost with smaller deviation, and the running time becomes higher. The experimental results on small datasets (Appendix B.2) suggest that the proposed MLSP method not only outperforms other algorithms in terms of clustering quality but also runs much faster than the BB solver. The experimental results on the performances with different values of k (Appendix B.3) show that our proposed MLSP algorithm can still achieve the best clustering quality for smaller values of k. The experimental results on the performances with fixed time limit (Appendix B.4) show that our proposed MLSP algorithm achieves the best clustering quality within any given time constraints.
Tables 12 and 13 show the results on 18 datasets with sizes smaller than 50, 000 using k 10. It can be seen that, for each dataset, the best clustering cost returned by our MLSP algorithm is smaller than BB method and other local search methods. For each dataset, the average clustering cost returned by our MLSP algorithm nearly matches the result of the BB method. As for FLS, although the performance of FLS is better than that of LS++, our proposed MLSP algorithm improves the performance of FLS on clustering cost with smaller deviation on most datasets. As for running time, there is no significant difference on running time among different local search algorithms for small datasets. Tables 14 and 15 present the performances of different local search algorithms with k on small datasets. Table 18 presents the performances of different local search algorithms with k- = 3 on large datasets. Tables 16 and 17 present the performances of different local search algorithms with & 5 on small datasets. Table 19 presents the performances of different local search algorithms with k = 5 on large datasets. It can be seen that, our proposed MLSP algorithm achieves the best clustering performance on most datasets with different values of k. On each dataset, the best clustering cost returned by our MLSP algorithm matches the clustering cost of BB method. On each dataset, the average clustering cost returned by our MLSP algorithm nearly matches the result of the BB method. As for running time, there is no significant difference running time among all local search algorithms on small datasets. However, as the data sizes grow, our proposed MLS algorithm becomes much faster than other local search algorithms.
5 Conclusion
In this paper, we propose fast local search algorithms for the k-means problem with multi-swap strategy, which runs in linear time in the data size. We develop new sampling techniques, which accelerate the process of clustering cost update during swaps. By proposing a recombination mechanism, the proposed algorithm can find potentially better solutions. Experimental results show that our algorithms achieve better performance on both small and large datasets compared with the state-of-the-ee the- algorithms. An interesting future direction is how to design fast local search approximation algorithms for handling high dimensional clustering datasets.
Acknowledgments
This work was supported by National Natural Science Foundation of China (62172446, 62350004, 62332020), Open Project of Xiangjiang Laboratory (22XJ02002, 22XJ03005), and Central South University Research Programme of Advanced Interdisciplinary Studies (2023QYJC023). This work was also carried out part using computing resources at the High Performance Computing Center of Central South University.
10

---

References
[1] Ankit Aggarwal, Amit Deshpande, and Ravi Kannan. Adaptive sampling for k-means clustering. In Apprximmation, Randomization, and Combinatorial Optimization. Algorithms and Techniques, pages 15-28. 2009.
[2] 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 Computing, 49(4):17-97, 2019.
[3] David Arthur and Sergei Vassilvitskii. k-means++: The advantages of careful seeding. In Proceedings of the 18th Anmual ACM-SIAM Symposium on Discrete Algorithms, pages 1027- 1035, 2007.
[4] Davin Choo, Christoph Grunau, Julian Portmann, and Václav Rozhon. k-means++:: Few more steps yield constant approximation. In Proceedings of the 37th International Conference on Machine Learning, volume I 199, pages 1909-1917. PMLR, 2020.
[5] Vincent Cohen-Addad, Hossein Esfandiari, Vahab S. Mirrokni, and Shyam Narayanan. Im- proved approximations for euclidean k-means and k-median, via nested quasi independdd sets. In Proceedings of the 54th Annual ACM SIGACT Symposium on Theory of Computing, pages 1621-1628, 2022.
[6] 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 Journal on Computing, 48(2):644-667, 2019.
[7] Zachary Friggstad, Mohsen Rezapour, and Mohammad R Salavatipour. Local search yields plas for k-means in doubling metrics. SIAM Journal on Computing, 48(2):452 -800, 2019.
[8] Fabrizio Grandoni, Rafail Ostrovsky, Yuval Rabani, Leonard Schulman, and Rakesh Venkat. A refined approximation for euclidean k-means, Information Processing Letters, 176:10625,,
2022.
[9] Junyu Huang, Qilong Feng. Ziyun Huang, Jinhui Xu, and Jianxin Wang. FLS: new local search algorithm for K-means with smaller search space. In Proceedings of the 31st International Joint Conference on Artificial Intelligence, pages 3092-3098, 2022.
[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. In Proceedings of the 18th Annual Symposium on Computational Geometry, pages 10-18, 2002.
[11] Silvio Lattanzi and Christian Sohler. A better k-means++ algorithm via local search. In Proceedings of the 36th International Conference on Machine Learning, volume 97, pages 3662-3671. PMLR, 2019.
[12] Stuart Lloyd. Least squares quantization in PCM. IEEE Transactions on Information Theory, 1292)11 1982.
[13] Konstantin Makarychev, Aravind Reddy, and Liren Shan. Improved guarantees for k-means++ and k-means++ parallel. In Proceedings the 34th Conference Neural Information Process- ing Systems, pages 16142-16152, 2020.
[14] Yusuke Matsui, Keisuke Ogaki, Toshihiko Yamasaki, and Kiyoharu Aizawa. Pqk-means: Billion- scale clustering for product-quannized codes. In Proceedings of the 25th ACM International Conference on Multimedia, pages 1725-1733, 2017.
[15] Jiayang Ren, Kaixun Hua, and Yankai Cao. Global optimal K-medoids clustering of one million samples. Proceedings of the 36th Conference on Neural Information Processing Systems, pages 982-994, 2022.
[16] Dennis Wei. A constant-factor bi-criteria approximation guarantee for k-means++. In Pro- ceedings of the 30th Conference on Neural Information Processing Systems, pages 604-612, 2016.
11

---

A Missing Proofs in Section 3
Lemma 2. Given type-1 rype-2 matched swap pair (g), h. it holds that <(P, C, Cj: 1 2A (*) + AA(P, C). Given type-1 matched swap set (Q, V), it holds that C(P. (,
(*)
Proof We first consider a type-1 matched swap set (Q, V). For a data point E P, we use Sp to denote its closest center in C. Let op be its closest optimal center in C*. Observe that
E
d(p. Sop) dpp.
(EX(V)(2J((I)
2 (VlD. + /(lpp 300)
E
ddp, 8p)
PEX(VII((IO)
(WuD.oo) + +2o)++/d(@e vd(oms Sp)
2 d(p, sp)
E
(EX(VIZ((IO)
2 d(p,sp) sp)
E PEX(VIZ((ION
2(d(p/o))) + v vd(p,sp)
2 od) +2 7 22dd(p, oy)d(p, sp)
E
(EX(1)Z((1I)
E
(4+ dd(p op) + +2d(pp sp).
PEX(VIZ((IO)
where the first inequality follows from the fact that So, is still in CVV for each data point p e (JQQ) according to the definition of the swap set defined, the second and fourth steps follow from the triangle inequality, the third step follows from the fact that So is the nearest point to op in C, and the last step follows from Cauchy Inequality. Let À to- Then <(P, C, Q, V) s 24)(X(V), C*)) + & (VX((V) The cases for matched swap pair are similar to the cases for matched swap set. For a matched swap pair (ch, Cj), we can get that C(P, C, c, S 24A (P,, C*) AA(P, C).
Lemma 3. Let PK be an optimal cluster with A(P,, b) (ct))) real number 2 3. Then, A (xx(P), C) 2 200 (b (P% (G,,)-
Proof Let 8c;, be the closest center in to cj. Observe that d(cr, sct) >> (0-1)A(P%.(GH)) Otherwise
[PA
A(Pr, {8et }) < bA(Pi) (ct)) holds according to Lemma 1. Consider an arbitrary data point
p E x(P)). By triangle inequality, we have
A(Pr, ((r)))
VA(tp).C) 2 dei)) vdlp,
(Vb 1V-1-- -V1.)-
Since 16 22, we have
3
(b () (GH))
a((p).022
1 1-V V
PA
Together with the fact that Ix(P)) 2 } P. we have ((x(P)), Z proves the Lemma.
1-1
A(P), (cf)). which
200
Lemma 4. By swapping with Cp the clustering cost of A(P. C) can be reduced at least by 1 - +
12

---

Proof Observe that A(P) (4)) < ( ( holds by Lemma 1. Then, we can get that
A(P,CUfc)UUal) = -
(A(P) +) + C) A(P), (Mgg (a))
& 1 A(P.C)
c.) (
S (1 100k
where the second to the last inequality follows from the fact that A < <9A(P% (.)))) and the last inequality follows from the definition of good single cluster.
Lemma 5. If case (2) happens, then ((X(L), S (1 () (4 (((()))))) + 100 A(P.C))
Proof By the condition of case (2), it holds that all the optimal clusters in (LL belong to bad single cluster. For any bad single cluster P% J[),, let (c), c% E M3 be the type-2 matched swap pair with minimum reassignment cost, i.e., (c) % arg (o.,o) C(P, C, c Cn). According to the definition of bad single cluster, we have
APP 1 A(P.C) (G3) ((P.C.c,.C))
100k
Let (MM)) be the set of lonely centers used for constructing type-2 matched swap pair. For an optimal
center of E C*, denote Seg the center in C closest to c Let L1 Uegez By the definitions of ((Ms) and L, for each Sej, E L1, we can find set (scc;) S x(Ms) with size ! 1p-11 ((s;))
such that z(a) n z(b) = 0 for any a, E L1. For each Seg L1. since (s.; } > t, it holds
that N - )[()((s;))))))) < 1 + 1. By taking summation over all the centers in L, we have
[L] EegEL 1 Σey EL, (s.;() Ee, EL, [=(se;)((1 + 7) s |s(MM))((1 !). Then, by considering all the optimal centers in L and taking summation, we have
A(X(L),C)= Σ ((P.C)) Σ A(P.C) + + (2t1)(//)))))))))()) + C(P,C,
100k GEL
EL
EnEN(M)) S(P,C, (k(M3)
< -XP0))
1
Σ
100
SEL
s (1+ (ONNZ(J)))CCC)+2)XXXXX
1 + -(X()))))) 100
1 A(P.)))
where the second the last inequality follows from the fact that cf is the center with minimum reassignment cost, and the last inequality follows from Lemma 2 and the fact that [|)]s <
1 1+t
Lemma 6. Let Q E H2 be set of centers in H2, where J(Q) is a set of good t-clusters with
type-1 matched swap set (Q, A') E M2. Define V AA" {c; where oj the center in A' with
[ j (c;) 1. Let U C P be the set of data points with IU V such that U n (x((P)))))))))))))))))())) # 0 holds
for each P E QT. Then, U(P.C'VUU) I (1 Took ((P.C).
Proof Consider an arbitrary optimal cluster P7 E Q's- We know that A(P C) < 3A(PH (ci)) by the definition of Q'g. Let Sc; be the closest center in C' to ch. We first argue that Se; is close to cj, with d(c )e7)) < SA(P; (c..)) For data point E P let 8p be its closest center in C. According
13

---

to the triangle inequality, we have
Σ pe sp)
d(c Scc)) -
Σ
+ (1+ 1 x (d(p
1
PP
PEPA
(1 + () (0.)) + (1+ ()A(P,
P%
(A R. (2,())))
PR
where the last inequality follows from A(P), < 3A(P) (c),) and feeding = v3. By assigning all the data points in Pž Seg, we have A(P) (8c;)) }) S 9A( (ct,) by Lemma 1. Then, by swapping U with V, we have that
A(P,C/VUU) A - -AP.0O
< A(P.C) - -
S - -
1 AA(PCC)
S 100k
where the first inequality follows from the fact that se; & V for each P7 E Q's, the second inequality follows from ((Ph) {sez }) < 9 h {er, }) for each PK E Q's, and the last inequality follows from the definition of good t-clusters and V [sag) = A'.
Lemma 7. Given set Q E Hg of centers such that JQ)) is set of good t-clusters, the probability that step 4 of Algorithm 2 fails to sample a data point from x(P)) for each P E Q's in the
mL P iteration is at most (1 + 30000k ) A112///2.CC////00PPP...
Proof We first consider an optimal cluster P% E Q'g. Observe that APP, C) < 30000k2t-r A(P,C)
by the definition of Q'g. Define p as the failure probability of not sampling data points from x(P)) in the m( -tth iteration of step 4 in Algorithm 2. Then, we have that
p = A(x(P)). A(P.C)
A(P.C)) A(P, C) 2 1
1
2 21-
30000k-'-11
Let 1/ 1 1--1-1 p be the success probability. Then, we have
p
30000
30000k 1) 1)
k-2
30000k2
30000
where the last inequality follows from 30000 k. Let p (Q') denote the probability of not sampling data point from (Pt) for each PV E Q' in the ml --h iteration of step 4 in Algorithm 2. For each P E Q'x. define X₁,
the m (Px))tt iteration of step in Algorithm samples a data point g E x(P)), and X₁ 0 if the m m(P))-th iteration fails to sample a data point q E x(P)). Let E, be the event that Xa = for each P7 E Q'g, and Eg be the event that there are at least two clusters P7. PJ E Q² such that Xn, V X I and X A X,
0. Observe that there are at most 2' 2 subcases if event Fe happens. For each subcase 8, let n(s)
(P% : PK E Q², X = 1} and (s) = (Pr P" E EQ', Xn = 0}. respectively. Define Pr(s) as the probability that subcase happens. Since the sampling iterations in step 4 of Algorithm 2 are mutually independent, we have (s)) = III: p Ip; (s) P/. By applying the inequality that pt < 30000k 1- P/) we
have
(s(s) S 300000 - 1)'
II
S
P
II
P
30000A - 1)
PIEQU
PIEQU
for each subcase of Fez, where the last inequality follows from 300004(2-1111 1. Then, we have (x() < 30000k Irgegg p/ by taking a probability summation over all 2 -2 subcases of event
14

---

Eg. Observe that the failure (robab)lityyy((9)) = P.(E))+PP((E)) Thus, can get that
where (,(E1) = II Pg eqg p/.
Py(('s) S (1 +1+ 30000k
II p) = (1
I
(1 (1)
30000k
17 E4'S
17 E4%
(1 A(x(P,)
0)
II (1 APP, C)
= (1 + 30000k
II PREQU
S (1+ 30000k
A(P.))
300A(P,C)
PIE
S (1 30000k
where the second to last inequality follows from Lemma and the last incquality follows from the inequality 1 6
Lemma 8. For the optimal clusters in HG, we have ((Z(H)), TOO A(P.C).
Proof Given any type-1 matched swap set (Q, V), if (QQ) is set of good t-clusters with type-1 matched swap set (Q: V) E M2, according to the assumption subcase (2), we know that each optimal cluster P% in Qu is bad single cluster with swap pair (c%,, ci), where is lonely center V. Since J(Q) is a set of good t-clusters, not all optimal clusters J(Q) are bad single clusters. Hence, for each optimal cluster P Qr. can find a lonely center ((c)) E V such that P7 is bad single cluster with the swap pair (c, (())). By the definition of bad single cluster, we have
A(P C)
+9 ((PC,
100k
For each optimal cluster P% Q's, by the definition of Q's holds that A(P; C) 3A(P%, (c))). ((Q) a set of bad t-clusters, by the definition of bad t-clusters, we also have
C) 100k
A(P.C) (A. + ((P.C.O.V)
Define H = (cr, E H1 P% is a bad single cluster with the swap pair (ch, cj) E M1). H3 (Q H2 J(Q) is a set of bad t-clusters with the swap set (Q, V/) E M2) and H3 (Q's UQ: QE ((Q) a of good t-clusters with the swap set (Q, V) E M2). respectively. Putting all things together and taking summation over all clusters in HB, we can get that
A(Z(HD),C)= Σ A(P..C)
Σ A(P,, C)+ Σ A(Z(J()).C) Σ A(Z(Q)).C)
QEB
Q'EH?
Σ
A(P.C) 9 (ct) (AA(
C*) ((((H))))))
100k
Pf
+ 244(X (H) + A(X(H)).))
24A(((())) 1 A(((H))) {
< 100
+ 24)(X(H)) CC*) + {
15

---

where the first inequality follows from Lemma 3. Then, together with Lemma 5, we can get that
(Z) A(P,)) (Z(([()), (()()))
2 A(PC) (1 -
+ 1 A(P.C)+
1 - C
100
100
+ ((((H) ()) + 24 (1) A(X(H)) C))
5
2 A(P,C) - (1 )
(P.()) -(P.C) 330pt)
50
33
33
2 A(PCC A (P.C)
(P.CC)
100
50
1 A(P.C),
N 100
where the second inequality follows from Lemma 5, the third inequality follows from the fact that Hi n H2 n L' = and Z ((H)) n ((L) = 0. and the fourth inequality follows from 2 2 and the assumption that (PP, C) 2 50(1 + p)pll.
Theorem In the i-th iteration of Algorithm let C' be the set of centers obtained in step the current clustering cost A(P, C) is larger than 50(1 + +)tpt, then with probability at least O(k-"), ") we have A(P, C²) S (1 1 (P,, C). After O(k((t) log(e- log iterations, we get an approximate solution with ratio (50(1 + ?)) + c) in expectation.
Proof Let T = 5k++1 1og(24c- log k), where ö is sufficient large constant, Following the work in
[]]. we define another random process X with initial clustering cost APP, for set C" returned
by the k-mmans++ algorithm such that for T" iterations of sampling and swaps, it reduces the value of
A(P,C" by least (1 Toog with probability Ak: - and increases the final value of A(P, C') by 50(1 + OT,, where is a constant with X < 100 is easy to see that E[A(P, < E[X].
Then, have
E = 50(1 + 10ppt
+ A((P,C' Σ
(11
1
1T' (1
AET
AA
100k
A(P.C' 100 ST +"550 ++5011+ + 10ptt
t
eA(P, + 50(1 + 1ptt.
24 k
This implies that E [[A(P,))CC
cA(P.C' + 50((1 + Opt. Then, we can get that ll
E[A(P.C) Σ EELIP.
P,(C')( <A(P.C")
24 50(1 + )10)t)
(F(((P,C)) + 50(1 + 1 24 log k
Since the k-means++ algorithm provides an approximation ratio of 8(log + 2) in expectation by Theorem 1, we have E[A(P. C)] < (50(1 + 2)) + O)tpt.
Running Time Analysis. By Theorem 2, in order to obtain (50(1 + :) + c) approxxiie solution, the iteration rounds for Algorithm should O(k*+1 log(e -1 log k)). In each iteration, it takes ((ndk) time to update the distances between data points to their closest centers. During the sampling
16

---

process, t data points are sampled according to the D'-Sampling distribution to serve as the candidate set of centers for swapping in, which takes time O(t) the distances from data points to their centers are already known. It takes O(k*) time to enumerate each subset with size at most i of the set of current centers opened. takes (mndk") time to recalculate the clustering cost after each swap nearest centers of each data point are maintained during the whole process. Thus, the total running time of Algorithm is (nddk2++ log(e" log k)).
B
Complementary Experiments
B.1 Experiments with Different Parameter Settings
In this section, we present the experiments on the performances of our proposed MLSP algorithm with different parameters T, R', and r). We also compare MLS and MLSP algorithms with LS++ algorithm using different iteration rounds to show the parameter robustness different algorithms.
Tables 3, 4 and 5 show the comparison results of different algorithms with varying iteration rounds on dataset rds, KEGG and Urban_ 10, respectively. The results show that, our proposed MLSP algorithm always achieves the best clustering cost compared with LS++ and MLS algorithms. It can be seen that the number of iteration rounds has little impact on our MLSP algorithm. For MLS algorithm, performs better on clustering cost compared with LS++ algorithm for most cases, which indicates that larger iteration rounds will not influence the performance of our algorithms.
Tables 6, 7, and 8 show the results of MLSP algorithm with varying number of failure upper bound R' for fixed c 0.5, 17 = 0.5 and T 400 on dataset rds, KEGG and Urban_10 respectively. can be seen that larger failure upper bound will lead to smaller deviation, and the running time becomes higher. In general, the failure upper bound does not influence the performance of our proposed MLSP algorithm. Tables 9, 10, and 11 show the results of MLSP algorithm with varying parameters and y for fixed 7 400 and R' 5 on dataset rds, KEGG, and Urban_10, 10, respectively. can be seen that smaller values of and y result in better performances on clustering costs with smaller deviation, and the running time becomes higher. There is no significant difference between different choices of parameters and 17, which indicates that our proposed MLSP algorithm ensures the parameter robustness.
Method Rounds
Cost
Average
Time(s)
Method
Cost
Average Cost
Time(s)
LS++ MLS
133.4074
139.672413.5172
2.56
LS++
133.1909
135824-2.0054
3.81
100
133.9806
139.6108++I.0053
1.12
MLS
200
132.5915
138.6551#2.3461
2.27
MLSP
131.7410
132.9932:10.9546
29.02
MLSP
131.7410
/33.1258:222266
30.25
LS++
132.0117
137.2619::2.6 2544
5.49
LS++
132.1966
138.25261:2.089
MLS
300
133.5915
138.6551++2..05
3.55
MLS
400
131,8264
136.9465 1638
4.48
MLSP
131.7410
132.3278:10.866
42.67
MLSP
131.7410
132.0323-10.4511
58.61
LS++
133.6918
140.7157442.0
0771177-2200022 .0066
12.26
LS++
131.7681
137,7633- 3.026
17.24
MI.S
500
131.8248
138.3528 +-5,2943
MISS
600
137.1133+
MISP
131.7410
132.85599+1..3888
66.73
MLSP
131.7410
132.1285
72.64
LS++
131.8519
137.8405 +3.31 80
22.83
LS++
131 131.8338
139.067442.5677
24.23
MLS
700
131.8493
136.1385
7.66
MLS
800
131.7709
135.7753--..085
7.98
MLSP
131.7410
132. 0.8551
81.88
MLSP
131.7410
132.2103 0.73338
91.52
Table 3: Comparison results on dataset rds with varying iteration rounds
17

---

Method
Rounds
Best Cost
Average Cost
Time's)
Method
Roumds
Best Cost
Average Cost
Time(s)
LS++
6.18718407
6,2712E-0711.2E400
2.13
LS-+
6.1619E+077
6.3341E+07-1 068
MLS
100
6.1557E407
7:25:1-40+111E4077
2,11
MLS
200
6.1684E+07
6.1154E4772222-44
4,00
MLSP
6.1534E+07
6.1599++0::49--000
22.22
MLSP
6.1534E+07
6.1534++07=0
LS++
6. 614424407
6.2732E-07+
10.41
LS-+
6,1547E-07
6:3276407=1
49
16.58
300
6.5379E=-7-7722Fxw0
400
6.1547F+07
6.3255E;07
+.8
MLSP
6,1534E+07
6.15599-447::449----0
46.18
MLSP
6,1534E+07
6.1546E+07-3.7E
57.13
6.1547E+07
6.28911447711..220
6.1547E+07
6,2878E+07
MLS
6. 1661
6,3595F:-07 1.7:446
10.54
MI.S
6,1534E
12.95
MLSP
6,1534+07
63.61
MLSP
6.1534E+07
6.1572E+07
5.7E::44
77.49
LS++
6.1534E+07
6.2874E-071
43.03
LS-+
6,1547E+007
6.2261E:07 5.99E005
54.28
MLS
700
6.1987E+07
6:1113E-07334E446
13,97
MLS
800
6.1953E+07
6.5863E407-3.54440
16.22
MLSP
6.1534E+07
6.1547E+00=337B440
86.74
MLSP
6.1534E+07
6.1559E+07=49E-0
96.73
Table 4: Comparison results on dataset KEGG with varying iteration rounds
Method
Rounds
Best Cost
Average Cost
Time(s)
Method
Rounds
Best Cost
Average Cost
Timo(s)
LS++
24659.3733
25255.4895 4895 1613.506
3.08
LS++
24660.9511
25024.9944 356.9319
8.36
MLS
100
24787.2318
30022.6100+334.100 1652
1.88
MLS
200
24659.389
25022 5201-332--777
3.34
24726.3142
MI.SP
LS++
24659.3733
25594.2717:19 12..222222
15.99
LS++
24659.0733
25383 49144997777777 577.4748
25.79
300
400
45166.5432--526.599
MLSP
24659.0733
24761.5254
$4,50
MLSP
24689.0733
24793.0067:22996877
58.74
LS++
33639.3732
25217.6713:--46 9049
37.94
LS++
24659.0733
25282 12282_5551=727-0 8617
MLS
500
24639. 3732
24835, 88 16-+-111 3455
7.79
MLS
24639. 3732
9.18
MLSP
24659.0733
24766-426++1188666888
73.87
MLSP
24659,0733
24742.35133:54:4211
76.58
LS++
34339
24852, 2682.1 446. 1449
69.19
LS++
24639 0844
2561 1.111,.111.-111 891.9473
.4733
88.
MLS
700
4008
25073,7342=295 8309
10.61
MLS MI.SP
800
24659 5079
25001 2214-223433555 35 15511.00 8913
12.01
24659.0733
134.1987
Table 5: Comparison results on dataset Urban 10 with varying iteration rounds
Failure
Best Cost
Average Cost
Time(s)
131.7408
133.8048+2.5306
39.49
131.7408
132.6853++1.5759
48.01
131.7408
132.8304+11.6403
59.41
6
131.7408
132.7752++1.1013
55.65
131.7408
132.2285-+0.8758
60.03
8
131.7408
132.4954400.6819
69.60
131.7408
132.30644106885
72.11
10
131.7408
132.2246:10.6323
85.36
Table Results of MLSP Algorithm on dataset rds with varying number of failure upper bound R'
Failure
Best Cost
Average Cost
Time(s)
3
6.1534E+07
6.1559E+07+4.9E544
44.52
4
6.1534E+07
6.1559E+07+4.9E:44
59.61
6.1534E+07 5
6.155999407144444444
64.45
6.1534E+07 6
6.1559E=0714499-444
70.69
6.1534E+07 7
6.1534B+07:10
76.37
8
6.1534E+07
6.1534E+0740
91.05
9
6.1534E+07
6.1534E+07+00
100.82
10
6.1534E+07
6.1534E+0740
105.89
Table 7: Results of MLSP Algorithm on dataset KEGG with varying number of failure upper bound
R!
18

---

Failure
Best Cost
Average Cost
Time(s)
3
24659.0733
24824.3035+214.2367
41.44
4
24659.0733
2481900528:238.56888
56.65
5
24659.073)
24755.9209+132.11
59.97
6
24659.072 07339
14806.0814+155
66.19
7
24659.075
1900000
80.24
8
24339.0733
44679.49471:88.7733
88.09
9
04659.0733
24703. 7859:15611331
92.57
10
24659.0733
24684.7586151.1044
95.37
Table Results of MLSP Algorithm on dataset Uran 10 with varying number of failure upper
bound R'
eta/epsilon
Best Cost
Average Cost
Time(s)
0.25/0.25
131.7408
1.2.687551.9399
61.57
0.5/0.25
131.7408
133.5867+2..0229
58.97
0.75/0.25
131.7408
133.0508++1.224
63.59
0.25/0.5
131.7408
132.330
50.81
0.5/0.5
131.7408
132.3244+1.5983
51.40
0.75/0.5
131.7408
+32.6404
50.50
0.25/0.75
131.7693
134.4831+:22116
52.86
0.5/0.75
131.7408
132.3092::0.9371
61.99
0.75/0.75
131.7408
133.209342.1465
60.73
Table 9: Results of MLSP Algorithm dataset rds with varying parameters and !]
eta/epsilon
Best Cost
Average Cost
Time(s)
0.25/0.25
6.1534E+07
6.1572E-07:5.7E-000
85.99
0.5/0.25
6.1534E+07
6.1540E+07+33775:44
78.57
0.75/0.25
6.1534E+07
6.1572E+0745.77--44
72.15
0.25/0.5
6.1534E+07
6.1572E+07+5577--44
68.67
0.5/0.5
6.1534E+07
6.1540E+07+3.77-04
65.05
0.75/0.5
6.1534E+07
6.1559E-07:4:99E-0
62.09
0.25/0.75
6.1534E+07
6.1559E+07+4.9E--0
65.15
0.5/0.75
6.1534E+07
6.1546E-00443777:-4
63.55
0.75/0.75
6.1534E+07
6.1546E+0743.7E+44
64.35
Table 10: Results of ML.PP algorithm dataset KEGG with varying parameters and
eta/epsilon
Best Cost
Average Cost
Time(s)
0.25/0.25
24659.0733
24681.3945:240.9771
66.64
0.5/0.25
24659.0733
24837.5858:223.333
72.86
0.75/0.25
24659.0733
81786.8942-1168.0088
59.95
0.25/0.5
24659.0733
24671.9425+3883418
60.31
0.5/0.5
24659.0733
24756.0949-133.5952
63.27
0.75/0.5
24659.0733
24746.6868+1144
56.24
0,25/0.75
24659.0733
1171
57.11
0.5/0.75
24659.0733
24804.6754
54.95
0.75/0.75
24659.0733
44753.6144-204.220
56.56
Table 11: Results of MLSP algorithm dataset Urban 10 with varying parameters and
19

---

B.2 Experiments on Small Datasets
In this section, we present the experiments of the performances of our proposed MLSP algorithm on other datasets with sizes smaller than 50, 000 used in [15].
Table 12 and Table 13 show the results on 18 datasets with sizes smaller than 50, 000 using k 10. can be seen that, on each dataset, the best clustering cost returned by our MLSP algorithm is smaller than BB method (note that BB method already guarantees a gap smaller than 0.1% to the optimal solutions on small datasets) and other local search methods. On each dataset, the average clustering cost returned by our MLSP algorithm nearly matches the result of the BB method. As for FLS, although the performance of FLS is better than that of LS++ our proposed MLSP algorithm improves the performance of FLS clustering cost with smaller deviation for most datasets. As for running time, there is no significant difference on the running time among different local search algorithms for small datasets.
20

---

Method
Size
BB(Cost)
Best
Mean
Time(s)
LS++
30.01
30.45+0.53
0.04
FLS
29.79
29.91=0.14
0.15
LLOYD
Iris(150*4)
29.79(735s)
30.32
30.90:10.4
0.01
MLSP
29.74
29.76:10.33
6.61
MLS
29.93
30.21
0.17
LS++
216.36
219.43+
0.37
FLS
214.52
216.68+2.88
1.78
LLOYD SEEDS(210*7)
218.49(448s)
215.65
219.99-2.93
0.12
MLSP
214.52
215.08
8.07
MLS
214.95
2.822
0.46
LS++
251.86
254.17
0.37
FLS
251.86
252.05
1.61
LLOYD GLASS(21449)
251.86(2566s)
253.25
259.94+3.99
0.11
MLSP
251.86
251.99
7.25
MLS
253.29
254.00=
0.46
LS+
377753
392980=14437
0.37
FLS
375974
377496-1934
1.71
LLOYD
BM(249*6)
375974(1204s)
376590
384475+5560
0.07
MLSP
375974
376276+265
6.81
MLS
378649
384982+4131
0.37
LS++
29.57
30.16
0.38
FLS
29.27
29.3
1.73
LLOYD
UK(258*5)
29.28(4h)
29.94
0..23
0.08
MLSP
29.27
29.29
7.14
MLS
29.41
30.00
0.39
LS++
6.96E+10
1.01
0.39
FLS
6.96E+10
7.06E+104331E:+009
1.77
LLOYD
HF(299*12)
6.96E+10(4h)
6.96E+10
6.97E-10+22.E--08
0.09
MLSP
6.96E+10
6.96E+10110
8.56
MLS
6.96E+10
7001++011.3E+008
0.4
LS+
1.36E+10
3.4.E+0016..E4408
0.4
FLS
3.36E+10
3.415+1014.55408
2.01
LLOYD
WHO(440*8)
3.40E+10(4h)
3.44E+10
3.49E=00-
0.11
MLSP
3.36E+10
3.37E+1014.6E+07
9.17
MLS
3,42E+10
3.52E-1045.6E400
0.4
LS++
11311E++06
L.454E++06-8.7+003
0.45
FLS
11331E++06
L.412E+00:224E+04
1.81
LLOYD
HCV(572*12)
1.1315E006(h))
L.1329E+06
L.1538E+0+-2.E--00
0.15
MLSP
1.1311E+06
I.1410E+00=222E--44
11.05
MLS
1.1505E+06
1.1135E00+442E+44
0.9
LS++
1.0786E+06
1.0933E+0+11.1E-04
0.51
FLS
1,0786++06
1.0816E40625.35+03
2.11
LLOYD
Abs(740*21)
1.0786E+06(4h)
L,0824E+06
1.209E++0+22.1E-44
0.12
ML.P
1.0786E+06
1.0786E+06+00
13.42
MLS
1.0789E+06
1.0979E+05+1.1E:000
0.47
Table 12: Comparison results on clustering costs and running time with k = 10 on datasets with
sizes ranging from 150 to 740
21

---

Method
Size
BB(Cost)
Best
Mean
Time(s)
LS++
776.04
793.7447.3
0.51
FLS
762.16
766.3844.9
2.14
LLOYD
TR(980+10)
772.47(4h)
770.72
776.4815.0
0.12
MLSP
762.16
764.71
11.89
MLS
785.01
805.89+
0.56
LS++
1.1734E+00
1.2124E
0.58
FLS
1.1734E+08
1.2092
2.35
LLOYD
2)C(1000*21)
1.1742E+08(4h)
1.1734E+08
L.1749E+08=11777000
0.09
MLSP
1.1734E+08
1.1752E+08:55114404
13.45
MLS
1.1735E+08
1.1846E+00:9.66-000
0.62
LS++
2.7116B+06
2.7657E+06144444444
0.58
FLS
2.7073E+06
2.8519E-00=1.224-04
2.51
LLOYD
HEMI(1995*7)
2.742IE+06(hh)
2.7144E+06
6.734
0.18
MLSP
+.7070E+06
2.123E+001667E-00
13.39
MLS
2.7292E+06
2.7886E#0616.9E-04
0.59
LS++
5.3689E+09
5.4772E+0911.0--08
0.58
FLS
5.3610E+09
5,4256E+0945.66-07
2.65
LLOYD
pr2392(2392*2))
5.3578E+09(4h)
5.3599E+09
5.3624E--09+:33E--0
0.14
MLSP
5.3578E+09
5.3668E+004423BB--7
10.76
MLS
5.3629E+09
5.4203E+09+6.2E---7
0.58
LS++
1.3902E+05
L4155E+05=11.55-03
0.62
FLS
1.3829E+05
1.4018E+05+1.1--03
3.43
LLOYD
TRR(5456
1.3796E+05(4h)
1.3868E+05
1.1177E405+2.15-003
0.18
MLSP
1.3796E+05
1.3829E
25.48
MLS
1,4591E+05
1.93999405
0.64
LS++
1167.6
1184.01
0.58
FLS
1163.7
11.4.1+13.06
5.06
LLOYD
AC(7195*22)
1181.7(4h)
1169
1172.1
0.36
MLSP
1163.7
1168.2
22.32
MLS
1165.5
1181.7
0.64
LS+
1.6104E+06
1.6508E+00-12444-000
0.24
FLS
166099E+06
1.61966-+00+773364444
2.26
LLOYD
cds_cnt(10000=4))
1.6119E+06(4h)
166009E+06
1.6247E-0061.44--05
0.11
MLSP
1.6009E+06
L61105E-4067.222-00
18.67
MLS
1.6146E+06
L.6520E+0012.555-05
0.42
LS++
1.8286++07
L8421E+07+111E+04
0.28
FLS
18829E++77
1.804E+00711.E-004
0.94
LLOYD
HTRU2((1898*8)
1.8723++074hh)
1.8266E+07
1.8380E+0742.06-00
0.15
MLSP
1.8266E+07
L.312E+07+42EE-00
24.43
MLS
1.8275E+07
L.869E+00:11.77-04
0.3
LS++
900239++06
9,1784E+064993E-04
7.06
FLS
8.9904E+06
9.0384E+0615.44-44
19.87
LLOYD
G1(36733*11)
8.9909E+06(hh)
9.0001E+06
9.1299E+0617.7E-04
0.33
ML.P
8.9904F+00
9,0131E--06+3..6---0
47.51
MLS
0.0239F+06
9.17844+06+9.3E+44
6.29
Table 13: Comparison results on clustering costs and running time with k = 10 on datasets with
sizes ranging from 980 36733
22

---

B.3 Experiments with Different Values of k
In this section, we present the experiments of the performances of our proposed MLSP and MLS algorithms on different datasets with different values of k. Tables 14 and 15 present the performances of different local search algorithms with k on small datasets. Table 18 presents the performances of different local search algorithms with k 3 on large datasets. Tables 16 and 17 present the performances of different local search algorithms with k on small datasets. Table 19 presents the performances of different local search algorithms with k = 5 on large datasets.
can be seen that, compared with other local search methods, our proposed MLSP algorithm achieves the best clustering performance on most datasets with different values of k. On each dataset, the best clustering cost returned by our MLSP algorithm matches the clustering cost of BB method (note that BB method already guarantees gap of 0.1% to the optimal solution). On each dataset, the average clustering cost returned by our MLSP algorithm nearly matches the result of the BB method. As for running time, there is no significant difference on running time among all local search algorithms on small datasets, However, as the data sizes grow, our proposed MI S algorithm becomes much faster than other local search algorithms.
23

---

Method
Size
BB(Cost)
Best
Mean
Time(s)
LS++
83.96
84.82++1.05
0.26
FLS
83,96
83.96:10
0.5
LLOYD
Iris(15004)
83.96(93s)
84.68
84.68:E0
0.09
MLSP
83,96
83.96:10
3,46
MLS
83.96
84.52:1:0.73
0.43
LS+
598.29
600.83++5.59
0.27
FLS
598.29
598.29--0
0.51
LLOYD SEEDS(210*7)
598.29(84s)
598.29
598.30 598.30:0..69
0..07
MLSP
598.29
598.29-:0
3.34
MLS
598.29
600.70:15.6
0.67
LS+
629.02
629.51::0.73
0.27
FLS
629.02
629.02:10
0.62
LLOYD GLASS(214*9)
692.02(107s)
629.02
629.02=00
0.05
MLSP
629.02
629.02
3.37
MLS
629.02
630.72::.088
0.68
LS+
8.63E+05
4.6
0.27
FLS
8.63E+05
8.63E+051:0
0.57
LLOYD
BM((49*6)
8.63E+05(86s)
8.65E+05
8.65E+05++0
0.03
MLSP
8.63E+05
8.63E+05+:0
3.69
MLS
8.63E+05
8.74
0.67
LS+
50.92
51.95
0.29
FLS
50.77
50.83+0.12
0.54
LLOYD
5K(258
50.77(89s)
51.29
51.31
0.04
MLSP
50.77
50.77+00
3.67
ML.S
50.77
51.78++0.81
0.67
LS++
7.83E+11
7.905+11+6.5E+10
0.27
FLS
7,83E+11
7.84E+11+2.0E+10
0.56
LLOYD
HF(299*12)
7.838+11(1075)
7.83E+11
7.83E-1116.3E+07
0.05
MLSP
7.83E+11
7.83E+11:33E+07
4.02
MLS
7.83++11
7.86E-1113.4E+09
0.74
LS+
8.33E+10
8.43E+10-1..E400
0.29
FLS
8.33E+10
8.33E+10:10
0.58
LLOYD
WH0(440*8)
8.33E+10(113))
8,40E+10
8.40E+ 10-L0
0.06
MLSP
8.33E+10
8.33E+10-:0
3.74
MLS
8.33E+10
8.46E+10-88.1 1E+08
0.71
LS++
2.75E+06
2.75E+06=10
0.32
FLS
2.75E+06
0.61
LLOYD
HCV(572*12)
2.75E+06(215s)
2.75E+06
2.795+0614.77+04
0.08
MLSP
2.75E+06
2.75E+06=10
4.04
MLS
2.75E+06
2.79
0.78
LS++
2.63E+06
2470E400+5S.EE+44
0.33
FLS
2.62E+066
2.62E+06=00
0.69
LLOYD
Abs(740*21) 2.62E+06(1 19s)
2.62E+06
2.62E+06-:0
0.07
SL
2.62E+06
2.62
5.4
MILS
2.62E+06
1.66
0.81
Table 14: Comparison results on clustering costs and running time with k = 3 on datasets with sizes
ranging from 150 740
24

---

Method
Size
BB(Cost)
Best
Mean
Time(s)
LS++
1145.77
1159.73=10.48
0.33
FLS
1134.45
1134.45:10
0.56
LLOYD
TR(980+10)
1134.45(126s)
1 1136.93
1138.09-10.94
0.09
MLSP
1134.45
1134.45:10
4.75
MLS
1134.45
1154.49: 14.61
0.78
LS++
1.28E+09
1.28
0.35
FLS
1.28E+09
1,28
0.82
LLOYD
2)C(1000*21)
1.28E+09(140s)
1.28E+09
1.28E+091:0
0.13
MLSP
1.28E+09
1,28E+09:0
5.18
MLS
1.28E+09
1.28E+09#1.1E+06
0.81
LS++
9.91E+06
9.91E=06:22.55+03
0.53
FLS
9.91E+00
9.91E+06:10
0.92
LLOYD
HEMI(1995*7)
9.91E+06(97s)
9.91E+06
9.91E+06=0
0.09
MLSP
9.91E+06
9.91E+06+0
5.44
MLS
9.91E+06
0.91E+06+2.3E++33
0.79
LS++
2.13E+09
1.
0.48
FLS
2.13E+09
2. 144
0.89
LLOYD
pr2392((392*2))
2.13E+09(123s)
2.13E+09
2.13E+09=00
0.09
MLSP
2.13E+09
2.13E+09::0
5.92
MLS
2.13E+09
2.15
0.76
LS++
1.96E+05
L.98E+05#1.7E+03
0.49
FLS
1.96E+05
1.97E+05#6.8E+02
0.97
LLOYD
TRR(5456
1.96E+05(325s)
1.96E+05
2.005-0553.334403
0.09
MLSP
1.96E+05
1.96E:005:3..8-+02
6.06
MLS
1.96E+05
1.98E+05+1..74003
0.78
LS++
2199.1
2234.68++4191
0.92
FLS
2199.1
2201.90+8.42
2.02
LLOYD
AC(7195*22)
2199.10(222s)
2227.18
2227.18:11
0.16
MLSP
2199.1
2205.75
6.43
MLS
2199.1
2212.69
1.06
LS+
1.49E+07
499++0712694
1.1
FLS
1.49E+077
1.49E+07+10
1.91
LLOYD
ads_cn((10000=))
1.49E+07(203s)
1.49E+07
1,99:+07+331122
0.03
MLSP
1.49E+07
1.49E+0710
4.96
MLS
1.49E+07
1.49+07+226.1
0.91
LS++
8.21E+07
8.22E+07+1..E+005
1.65
FLS
8.21E+07
8.21E+0710
3.71
LLOYD
HTRU2((7898*8)
8.215+07(1555s)
8.21E+07
8.21E+07+7.994+03
0.97
MLSP
8.21E+07
8,21E+07
8.44
MLS
8.21E+07
8.22E+07+2.0E+05
1.21
LS++
1.95E+07
1.97E-07+4.EE+05
4.38
FLS
1.95E+07
1.95E+07+6.0E+04
7.28
LLOYD
G1(36733*11)
1.95E+07(1936s)
1.95E+07
1.99E-07+5.4E+05
1.18
ML.P
1.95E+07
1.95E-407773B8+03
12.83
MLS
1.95E+077
1.97E-07+4.7E+-5
3.82
Table 15: Comparison results on clustering costs and running time with k = 3 on datasets with sizes
ranging from 980 36733
25

---

Method
Size
BB(Cost)
Best
Mean
(sme(s)
LS++
50.97
52.75=1.74
0.37
FLS
50.97
51.70=1.10
0.78
LLOYD
Iris(150*4)
50.92(355s)
51.19
51.19.10
0.1
MLSP
50.97
50.9
3.53
MLS
50.97
11.78++1.11
0.81
LS++
401.21
404.58++3.00
0.3
FLS
401.21
402.46+11255
0.81
LLOYD SEEDS(210*7)
401.21(376s)
401.21
406.69+2.23
0.11
MLSP
401.21
401.46110.55
3.62
MLS
401.21
409.6519.40
0.83
LS++
437.88
443.63:1.54
0.29
FLS
437.73
439.82.16.29
0.89
LLOYD GLASS(214*9)
437.73(592s)
437.88
457.22
0.07
MLSP
437.73
437.73+0
3,43
MLS
437.73
439.37+2249
0.86
LS++
6.0249E+05
6.22135+05+
0.29
FLS
6.0249E+05
6.1293E-05=5522--03
0.83
LLOYD
BM(249*6)
6.0249E+05(3895)
6.0249E+05
6.1293E+05+4.9E-003
0.04
MLSP
6.0249E+05
6.0249E+05:0
3.71
MLS
6.0249E+05
6.2244+-05+2.0E-00
0.83
LS++
40.17
40.97=00.3
0.29
FLS
40.17
40.28=0.25
0.82
LLOYD
UK(258*5)
40.17(457s)
40.24
41.16
0.06
MLSP
40.17
40.25 0.25
3.72
MLS
40.34
41.5
0.83
LS++
3.0998E+11
3.1257E++1+2.5E++10
0.31
FLS
3.0998E+11
3.1083E+11+111B=+10
0.89
LLOYD
HF(299*12)
3.0998E+ (0998E+11(1733))
3.0998E+11
3.001E+11=1.2E--00
0.06
MLSP
3.0998E+11
3.1021E+11+7.1E-008
4.65
MLS
3.0998E+11
3.1499++1155.6+-09
0.92
LS+
5.5982E+10
56642E+1043366--08
0.33
FLS
5.5914E+10
16333E++1+22.11400
0.89
LLOYD
WHO(440*8)
(591445+0(1880s))
5.5914E+10
5.6219E:+10+1..5550000
0.09
MLSP
5.5914E+10
5.6034E-10+1..55--00
4.84
MLS
5.5973E+10
566460E+1043.55--0
0.88
LS++
1.9716E+06
L.9722E+06-622.78
0.34
FLS
1.9716E+06
2.0044E-06-22.6+04
0.91
LLOYD
HCV(572*12)
199716E+0(122688))
1.9716E+06
1.9995E-0623333--00
0.13
MLSP
1.9716E+06
1.9716E+06:00
5.02
MLS
1.9716E+06
20113E+00+3.1E-00
0.93
LS++
1.7472E+06
1.7630E+0643.8E-44
0.36
FLS
1.7472E+06
1.7472E+0620
0.98
LLOYD
Abs(740*21)
1.7472E+06(410s)
1.7472E+06
17502E-06+1.5E+03
0.09
MI.P
1.7472E+06
1.7472E+06+0
6.03
MLS
1.7472E+06
1.7676E+00+3.5E-00
1.16
Table 16: Comparison results on clustering costs and running time with k = 5 on datasets with sizes
ranging from 150 740
26

---

Method
Size
BB(Cost)
Best
Mean
(ime(s)
LS++
966.64
990,43+14.69
0.35
FLS
953.39
956.16=:3.43
1.13
LLOYD
TR(980+10)
953.39(824s)
965.74
169.1045.65
0.14
MLSP
953.39
954.62:12.48
6.28
MLS
983.22
1000.86:1
0.93
LS++
4.6926E+08
1722688+08+2.8833
0.39
FLS
4.6919E+08
4.7197E+08+22.
1.06
LLOYD
S)C(1000*21)
4.6919E+08(220S)
4.6919E+08
4.6949E-08=2.EE44
0.09
MLSP
4.6919E+08
1.00
7.7
MLS
4.7204E+08
4.7529E+08-2.74-06
1.22
LS++
5.3811E+06
5,4497E+06:5.6640
0.45
FLS
5.3811E+06
5,4090E+00=3344--00
1.24
LLOYD
HEMI(1995*7)
(.3811E+06(9305)
5.3811E+06
544262E40=11.5EH44
0.13
MLSP
5.3811E+06
5.8880E+0622.1E--44
7.04
MLS
5.3811E+06
5.478144061.44544
0.86
LS++
1,1621E+09
1.7700++0966.3E407
0.44
FLS
1.1619E+99
L.1636E+0922.3E407
1.31
LLOYD
pr2392((392*2))
L.6619+09(6255)
11620E+009
L.1620E++0922.29--05
0.14
MLSP
1.1619E+99
1.1619E+09:00
7.46
MLS
1,1227E+09
L.1759E+0=11.4E-07
0.81
LS++
1.6920E+05
1.7055E+05::666.4
0.54
FLS
16894E+05
1.7012E+05+529.7
1.92
LLOYD
TRR(5456
L68700++0(30944)
1,7866++05
1.7250E+05
0.15
MLSP
1.877E+00
L.6903E+55+4733
9.57
MLS
1 1.691EE+05
0.7027E=05
0.89
LS++
1636.1
1636.93+0.99
0.67
FLS
1636.1
1636.14::0
2.6
LLOYD
AC(7195*22)
(63552s
1542.6
1642.56
0.24
MLSP
1636.1
1636.14110
7.68
MLS
1636.1
1637.25
1.57
LS+
5.3728E+00
5.3747E++0-22226--00
1.29
FLS
5.3725E++06
5.3725E+06::103
3.11
LLOYD
<ds_cn((10000=4))
5.3725E+06(7171s)
5.3744E+06
5.3751E+06:19887
0.05
MLSP
5.3725E+06
5.3728E+0615554
7.53
MLS
5.3728E+06
537422+06=2.11540
1.58
LS++
4,2176E+07
4.2200=40=44.EHH0
2.31
FLS
4.2152E+07
42160E-407+93E5403
5.03
LLOYD
HTRUU((1898*8)
4.2154E+07(4h)
4.2171E+07
4.173E+40-2299--03
0.16
MLSP
4.2152E+07
42153E+07-229E--0
11.01
MLS
4,2268E+07
42413E407=1.44444
1.56
LS++
1.3351E+07
1.3643E+07+3.86405
5.34
FLS
1.3351E+07
13438E+00+22.6E405
10.67
LLOYD
G1(36733*11)
1.3351E+07(4h)
1.3351E+07
1.330E+00733:44400
0.19
MI.P
1.3351E+07
1.3351=+07+0
16.61
MLS
1.3358E+07
1.3673E+07+3.9FH00
4.93
Table 17: Comparison results on clustering costs and running time with k = 5 on datasets with sizes
ranging from 980 36733
27

---

Method
Size
BB(Cost)
Best
Meun
Timets)
LS++
476.88
488.32
2.93
FLS
476.88
88+0
10.44
LLOYD
r)s(500000033)
476.79(811s)
476.88
507.69++31.5
0.18
MLSP
476.79
476.8=0.07
15,52
MLS
482.77
LS++
4.9412E+08
4.9412E:08::0
1.04
FLS
4.94121:+00
4.9412E:+08 +0
7.27
LLOYD
KEGGGS3,413133)
4:9412B+08(2901))
4.9412E+08
4.9412E=08:10
0.26
MLSP
4,9412E+08
4.9412E=08:10
12.89
MLS
4.94121x+08
4.9412E:+08:4
1.97
LS++
1.15E+05
1.1E++0=1..E:4
3.26
FLS
1.15E+05-:229
LLOYD
Urban [aca__101100000010 (.1SE+0S6833s)
1.15E+05
1,23E+05
1.14+04
0.18
MLSP
1.158+05
1.15
26.14 14
MLS
1.15E+05
1.16E+05 1,6E+03
LS++
L6396E+15
L.399E++1513.8E.
8.27
FLS
16339E++5
1.6397=+15
1.33+11
58.78
LLOYD
(NG_AGR(19))8477)
L64E+15(4h)
1.639EE+15
L.6398E-15+3.2E-11
MLSP
16396E+15
1.336E+15:0
53.25
MLS
1,6396E=+5
1.63981+15+ 133++11
7.73
LS++
414253
417302:14669
9.86
FLS
414733--956
LLOYD
Urban Urban_GB(360. 1773
OOM
414253
423509 12613
0.32
MLSP
414253
314
71.46
MLS
414253
417063316123
8.87
LS++
2.2779E+07
2.2781E-07 1824
13.93
FLS
2.2779E+07
2.2780E-07 777
142.15
LLOYD
SPNET_ 3PNET_3DX43487493)
OOM
2.2781E407
2.27722B-0
MLSP
2.2779E+07
2.27779::47
121.03
MLS
2.2779E+07
2.2781E+07++1022
10.33
LS++
3.8720E+06
388722=406::10
37.1
FLS
368720E+06
3.8721E=06::44,
338.19
LLOYD
5ym(1,000,00002)
OOM
3687200:006
187224:-06+2266-
0.69
MLSP
3,8720E+06
3.8721E=06:332.9
223.97
MLS
3.8720E+06
3.87224406:11992
22.93
LS+
6,9094E+08
08
286.41
FLS
6,9094E+08
6.9094E+08:10
3244.76
LLOYD USC 19956)2955886685666)
OOM
6,9094E+08
6.9094E:08-00
18.57
ML.P
6.9094E+08
6,9094F:008
2443.28
MLS
6.9094E+08
6,9094E+08 08-00
227.59
LS++
4.6299E=07
4.78918+07-
66F406
305.27
FLS
4.6299E+07
4.8013E+07:11.E-000
3894.45
LLOYD
SUSY(5,000.000
OOM
4.6299E:07
29.2
MLSP
4.6299E+07
4.7349E-0777.774404
2610.89
MLS
4.6299E+07
4.76:2E+07=1.6E:006
262.61
LS++
2.2614E+08
2.2726E+08:t11. IE+06
614.07
FLS
2.2614E+08
2.2720E-08-99.2 2005
8728.13
LLOYD
MIGGS(11,00,,00002
OOM
MLSP
2.2614E+08
6492.59
MLS
2.261 48:+08
2.27281:408
+1.2
522.15
LS++
1.5774E+13
1.6098E+13::1.6..M
31994
FLS
1.5649E+13
1.58
18306
LLOYD
SIFTT
OOM
1,5718B+13
1.5793E-11=:1.5E-11
401
MLSP
>48h
MLS
1.5718B+13 1.5848B413+111E4M
12738
Table 18: Comparison results on clustering costs and running time with k = 3 on datasets with sizes
larger than 50, 000, where 00M is short for out of memory
28

---

Method
Size
BB(Cost)
Best
Meun
Timets)
LS++
284.27
287.39:4:2.78
3.34
FLS
290.15
290.15:40
LLOYD
rds(50,000*3)
282.65(4h)
284,22
287.50+ 8.667
0.86
MLSP
283.13
283.27:10,21
58.61
MI.S
287.23+
LS++
L9900E=008
1.920SE+08:22
ELS
1.9201 E+008
1,92051=+08
1.61:+04
10,92
LLOYD
KEGGI(3,,13-23)
(.928+08(300s)
1.9201=+08
1,9754E+08:
1.6E+07
0.26
MLSP
1.9201E+08
1.9201E+08:10
24.83
MLS
E9900E=008
1.92045+08=
2.33--
1,88
LS++
5.62E+04
56281176
5,33
5.62E+04
56248--49
LLOYD
Urban Uiban_100000000 10(100,0002)
56231(4h)
5.62E+04
57768=2762
0.14
MLSP
5.62E+04
56246-:36
40.97
MLS
5.62E+04
$6315:83
LS++
5.0721E+14
5.0721E+1443.3E=111
15.65
FLS
5.0721E+14
5.0739E+14
122.09
LLOYD
5.0721E+14(hh)
5,0721E+14
4.0781E+14=3.2E11
MLSP
5.0721E+14
3.07366+14:311E-11
73.2
ML.S
5.0721E+14
5.0759F+14
14.24
LS+
201878
202554-954
16.48
LLOYD
Urban GB(360.17792)
OOM
201878
204943:17816
0.53
MLSP
201873
201874++2.03
103.94
MLS
201878
2021661629
LS+
8.8280E+06
18309E+0011.4E-40
26.85
FLS
8,8280E+06
8.8289E+06= 1.1
331.24
LLOYD
SPNET_3D (DNET_30(33487493
OOM
8.8280E+06
8.8315E406= 1.64403
MI.SP
8.8273E++06
8.82738-00--0117
180.12
MLS
8.8277E+07
8.8313E+06= 1.04403
22.94
LS++
88662EE066
1 1.869EE+0---30 3,,2
66.78
FLS
188755++06
1.86744+06
1.8E+03
$27.22
LLOYD
5y0(1.000.000"2)
OOM
1 86661E+06
1 1,95791+006 1.1E+03
0.85
MLSP
L8652E=006
1.8699E+0=11.. 1.1E+03
391.86
MLS
18858E+06
L.8944E-0633.14400
31.62
LS+
3,9017E+08
3.9375E+08::1 1
439.74
FLS
3,9017E-08
3.9017E+08:10
4802.6
LLOYD
USC_ 199002,458.685668)
OOM
3,9017E+08
3. 3861E+08==3612 3.64407
20.11
MI.SP
3,9017E:008
3,9017E++0811
3612.99
MLS
3,9017E+08
3.9732E408= 1.4
324.97
LS++
4. 4.06238+07
4,1035F+07
1.550-00 1.51:066
451.28
FLS
3,9880E+07
4.0997E+07=4.5E--44
5329.33
LLOYD
SUSY(5.000,000117
OOM
3.98488+07
4.0648E=+07 -,4
31.67
MLSP
3,8720E+07
3.2116++07:1.4E--55
4137.19
MLS
3,9904E+07
4.8852E-07:3.6E406
297.43
LS++
2.0829E+08
2. 1016E+08=1.3E-440
1316.59
FLS
2.0829E+08
2. 1012E+08=1.2222
18429.3
LLOYD
HIGGS(11,000000022))
OOM
MLSP
2.0829E+08
10100-00:1111E--4
12713.2
MISS
2.08299:+08
2.1018F408- 1,41:+06
1019.6
LS++
1.5035E+13
1.5.948++1+8..84410
41376
FLS
L4778E+13
1.9994+
27908
LLOYD MLSP
S0FT( 128)
OOM
14400S++33
1.5021E+13=8.00-10
531
>48h
MLS
1.4778 1.4921-13-8.9---00
12931
Table 19: Comparison results on clustering costs and running time with k 5 on datasets with sizes larger than 50, 000, where O0M is short for out of memory
29

---

B.4 Experiments on Different Local Search Algorithms with Fixed Time Limits
In this section, we conduct some additional experiments on different datasets with fixed time limits. The experimental results in Table 20 demonstrate that our proposed MLSP algorithm achieves the best clustering quality with fixed time limits.
Dataset
Mothod
Time Point
Cost
Time Point
Cost
Time Point
Cost
Time Point
Cost
137.96
137.93
FLS MLSP MLSS
138. 98
137.26
134,68
134,62
rds
20s
133,62
40s
60s
132.28
80s
132.14
140.19
138.23
137.63
137.63
62855194
62855194
62793066
62793065
FLS MLS
61633499
616,34999
61633499
61608744
KEGG
20s
61534479
408
61534479
60s
61534479
80s
615344779
63419827
63419827
634198827
6341 63419827
LS++
2501 62
666
2501 1.62
FTS
25 13103.8 83
25066653
2,5058.2 2.1
24973 42
10
20s
24767.85
408
608
80s
MLS
24920.38
24885.88 80
241855.86 80
24885 SU
LS++
1.3833B+14
1.3833E+14
1. 1.3833E+14
1. 138833+14
FLS
1.3822B+ 14
1,381 3E-
38008- 14
I 1.3797E 14
RNG_AGR
60s
120s
180
240x
1.3689++
MLS
1.38261:+1
L.3826E+14
1. 1.3826E:+1
1 1,3826E:14
LS++
89234,23
89234.23
892.4423
89234.23
FLS
89691.46
89691.44 .46
89688.77 77
89406.39
Urhan GB
60s
89314.46
120k
89311.78
180
240x
88407.99
MI.S
88955.80
88955.80
88955.80
88955.80
LS++
25763966
2576396
2575784
2575784
FLS
2573606
2572637
2571369
2570776
SPNET_A0 3D
MLSP
200s
2570513
400s
2570292
600s
2570160
800s
2569484
MLS
2574266
2574266
2574266
2574266
LS
56 1444.2
56 1444.2
56 444.2
561 444.2
FLS MLSP
561914.9
56 1914,9
561914.9
5619144.
300s
560502.5
600s
560502.4
900s
560502.4
1200s
560502.4
MLS
562385.8
562385.8 8
562385.8
562385.8 8
LS++
270910336
2709 (366
2709 10366
270910360
FLS
27091 10366
770910366
270825611 16
2708256
USC_1990
MLSP
2000s
270738773
4000s
270735773
6000s
270735773
8000s
270735773
MLS
270825616
270825616
270825616
270825616
LS++
32437956
32437956
32437956
32437956
FLS
3169312
31693
31693 28
31693128
SUSY
ML
2500s
31756413
5000S
31650548
7500s
31650548
10000S
31629757
MLS
32548680
32548680
32548680
32548680
LS++
189640888866 6
188373729.5
188373729.5
187098096.4 4
189186932.5
86932.5
189186932.5
189186932.5
HIGGS
MLSP
12000s
1840841229.5
24000s
184084129.55
36000s
184084129.5
48000s
184084129, 29.5
ML.S
186589935.5
186689935.5
186689935.5
86689935.5
LS++
1,82E+13
1.77E++13
1.67E+13
1,58E+13
13
1,52E:13
SIFT
MLSP
20000s
1.85E+13
40000s
1.85E+13
60000s
1.85++13
80000s
1 1885+13
MI.S
1.38E+13
1.38E+13
L338
1,38
Table 20: Comparison results on clustering costs with h = 10 and fixed time limits
30


Reviews:
[{"review_id": "wunfyiUN5S", "reviewer": "Reviewer_BGQc", "summary": "This paper proposes a multi-swap local search algorithm for the k-means problem with linear running time in the data size, while also achieves a better approximation ratio when compared with other local search algorithms that adopt single-swap strategy. To benefit more from such algorithm when handling large-scale datasets, the authors further propose a sampling-based method to accelerate the updating of clustering cost during the swaps. The extensive empirical experiments validate the good numerical performance of the proposed algorithms.", "questions": "Line 11-12, \"which improves the current best result 509\": It might be true for local search methods that take single-swap strategy, but It is unfair to call this the best result, since the current best approximation ratio is 6.357.\n\nLine 90, \"MLS algorithm\": Please refer it and write out the full name when mentioning this abbreviation for the first time. Please consider move the \"multi-swap local search algorithm\" at line 100 to here. \n\nAlgorithm 1: \"Output: ... of at most k centers\". I do not understand why you say \"at most\". Since in each iteration point p can be selected if and only if point p is not in C, isn't it guaranteed that C has size exactly k? Similar question to Algorithm 2.\n\nLine 3, Algorithm 1: $\\Delta$ is defined over two sets, here p is a point it should be {p}, and {C} should be C. This appears in other places as well, for instance Algorithm 2.\n\nIn the numerical experiment section, in Table 2 why is MLS able to outperform the single-swap local search methods LS++ and FLS by a big margin? It does not make sense to me since MLS algorithm is clearly more complex than LS++ in [10] and it has bigger runtime complexity.  \n\nLine 345-346, \"the average clustering ... nearly match the results of the BB method\": As far as I can see, except on one dataset, isn't MLSP algorithm always outperform BB method in terms of average clustering costs?", "limitations": "The authors did not address the limitation of their work. ", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 2, "strengths": "This paper studies the classic k-means clustering problem and gives a first linear time local search algorithm in literature that adopts the multi-swap strategy. Such algorithm enjoys the two advantage from both worlds: smaller approximation ratio and linear run-time in data size. The organization of this paper is clear and well-structured, the mathematical proof is very detailed and easy to follow, and the numerical experiments are also quite extensive. ", "weaknesses": "For Lloyd algorithm, according to your experiments it is able to save up to 99% of the run time (0.07 vs 5.12), while only having around 1% more average cost. Therefore, the merit carried by the work presented in this paper does not seem to be too impressive from the practical point of view. One related question is, considering the huge amount of time being saved from Lloyd algorithm, have the authors tried running Lloyd algorithm for multiple times (if Lloyd algorithm takes 0.07s for each initialization, then you can run 5.12/0.07 times of Lloyd algorithm) with different initialization and then simply take the best solution? How does that compare with MLS when they are set to have the same run time?", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the positive rating and the thoughtful comments. In the following we address the concerns.\n\t\n**Question 1: Line 11-12, \"which improves the current best result 509\": It might be true for local search methods that take single-swap strategy, but It is unfair to call this the best result, since the current best approximation ratio is 6.357**\n\t\nResponse: We thank the reviewer for pointing out this  issue. In the revised version of the paper, we will rewrite the statement as \"which improves the current best result 509 with a linear running time in the data size.\"\n\t\t\n**Question 2: Line 90, \"MLS algorithm\": Please refer it and write out the full name when mentioning this abbreviation for the first time. Please consider move the \"multi-swap local search algorithm\" at line 100 to here.**\n\t\nResponse: We thank the reviewer for the kind suggestion. In the revised version, we will give the full name of \"MLS algorithm\" when mentioning this abbreviation for the first time.\n\t\n**Question 3: Algorithm 1: \"Output: ... of at most k centers\". I do not understand why you say \"at most\". Since in each iteration point p can be selected if and only if point p is not in C, isn't it guaranteed that C has size exactly k? Similar question to Algorithm 2**\n\t\nResponse: We thank the reviewer for raising this subtle issue. \n\nIn the standard definition of the $k$-means problem, the goal is to output a set $C \\subseteq \\mathbb{R}^d$ with size at most $k$ such that the clustering cost is minimized. To be consistent with the literature, we phrase our output requirement as \"a set $C \\subseteq \\mathbb{R}^d$ with size at most $k$\". Since local search methods for the $k$-means problem in literature all return a solution of size exact $k$, we also follow the tradition and require our algorithms to return a solution of size exact $k$.\n\t\n**Question 4: Line 3, Algorithm 1: $\\Delta$ is defined over two sets, here $p$ is a point it should be $\\lbrace p \\rbrace$, and $\\lbrace C \\rbrace$ should be $C$. This appears in other places as well, for instance Algorithm 2**\n\t\nResponse:  Thanks for pointing this out. The $\\Delta(p,\\lbrace C \\rbrace)$ used in line 3 of Algorithm 1 is actually $\\Delta(\\lbrace{p\\rbrace},C)$, which is a typo. Sorry for the confusion.\n\t\n**Question 5: In the numerical experiment section, in Table 2 why is MLS able to outperform the single-swap local search methods LS++ and FLS by a big margin? It does not make sense to me since MLS algorithm is clearly more complex than LS++ in [10] and it has bigger runtime complexity**\n\t\nResponse: We thank the reviewer for raising this important question. Theoretically, the running time of the MLS algorithm is larger than that of LS++ and FLS. The key point for MLS to outperform  LS++ and FLS in experiments is that the sampling-based method is applied to accelerate the swapping and clustering cost updating processes (as mentioned in line 337).\n\t\n**Question 6: Line 345-346, \"the average clustering ... nearly match the results of the BB method\": As far as I can see, except on one dataset, isn't MLSP algorithm always outperform BB method in terms of average clustering costs?**\n\t\nResponse: Thanks for raising this question. Since the average clustering cost of our MLSP algorithm does not outperform that of the BB method on dataset RNG\\_AGR, for the sake of rigor, we did not state that our proposed MLSP algorithm outperform the results of the BB method. In the revised version, we will make a clearer statement as \"the average clustering costs outperform the results of the BB method except on one dataset RNG\\_AGR\".\n"}, {"review_id": "2IynEuiudw", "reviewer": "Reviewer_dZuQ", "summary": "Local search algorithm is a well studied technique for clustering problems. In the k-means problem, the simple swap heuristic states that start with an intiial chosen set of k centers  and then at each local search step, check if swapping an existing center with a new one leads to decrease in cost. It is known that this gives a constant factor approximation algorithm (with the constant being 25). The $t$-swap heuristic generalises this local search algorithm by swapping a subset of $t$ centers at each step. It is clear that the running time of such an algorithm will be exponential in $t$, though the approximation factor converges to 9 as $t$ increases. \n\nA related line of work on the k-means problem deals with the Lloyd's algorithm. In this algorithm, we choose an initial set of k seed centers and then run a different local search algorithm on it. A lot of work has happened on how to choose these initial centers. The k-means++ heuristic says that we choose the initial centers as follows: the first center is picked uniformly at random. After having picked $i$ initial centers, a point is chosen as the $(i+1)^{th}$ one with probability proportional to square of the distance from the $i$ centers chosen already. It is known that if we choose a set of $k$ centers in this manner, then it can lead to $O(\\log k)$-approximation algorithm for the k-means objective (and this bound is tight). Subsequent works have focussed on whether we can improve this initial set in a small amount of time. For example, reference [10] showed that if we run the single swap heuristic on the $k$ centers chosen by this random sampling procedure for about $O(k)$ steps, then we get a 509-approximation algorithm. \n\nThis paper is an extension of this last result by giving a fast implementation of the $t$-swap heuristic on an initial set of centers chosen using the random sampling procure described above. Each local search step chooses a set of $t$ new centers using a random sampling procedure and checks if it can replace an existing set of $t$ centers in the current solution. The paper shows that the running time now becomes $O(k^{O(t)} nd$, and hence is linear in $n$ for fixed $k,t$. They also show that the resulting set of centers has approximation ratio of about 50. The analysis has the same structure as that of [10], though analyzing $t$-swap is much more tricky, and requires new insights. ", "questions": "1. Please compare with [9] in terms of theoretical and experimental results. \n2. It may be worth explaining what new ideas are needed in the analysis as compared to [10]. ", "limitations": "No negative impact. ", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "1. Usually $t$-swap heuristcs are expensive to implement, this paper gives a non-trivial implementation which is efficient in practice. \n2. The analysis of $t$-swap requires new ideas. ", "weaknesses": "1. The improvement in experiments is only marginal (in most cases less than 5%). \n\n2. The usual local search heurstic (with single swap heuristic), ie.., reference [9] in table 1, has running time $O(nkd \\lof \\Delta)$ and has approximation ratio of 25. This running time almost matches the running time of the result in this paper (upto log \\Delta factor), and has better approximation ratio. So it is not clear if the theoretical contribution is significant. I also don't see a comparison with [9] (with single swap heuristic) in the implementation. ", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the positive rating and the thoughtful comments. In the following we address the concerns.\n\t\n**Question 1: Please compare with [9] in terms of theoretical and experimental results**\n\t\nResponse: We thank the reviewer for raising this question. We are sorry for the confusion caused by the typo on the time complexity of the algorithm in [9]. The term of $O(n^{t}k^{t}d\\log\\Delta)$ in Table 1 should be $O(n^{t+1}k^{t}d\\log\\Delta)$.\n\nOur proposed MLS algorithm can achieve a $(50(1+\\frac{1}{t})+\\epsilon)$-approximation in time $(ndk^{2t+1}\\log(\\epsilon^{-1}\\log k))$. It can be seen that even with a single swap strategy ($t=1$), the algorithm in [9] (denoted as LS algorithm for short) has at least quadratic running time in the data size while the running time of our proposed MLS algorithm has a linear dependence on the data size. For the LS algorithm, it is not surprising to get a better ratio with quadratic running time. \n\nFor experimental performance, as shown in [8] (the FLS algorithm), even with single-swap strategy, the LS algorithm can hardly handle datasets with size over 10,000 due to its quadratic time complexity. Thus, we did not include the results of LS algorithm for comparison in our experiments. In the following table, we give additional results (see table 1) to compare our MLS and MLSP algorithms with the LS algorithm on different datasets used in our experiments with size smaller than 10,000. To ensure a fair comparison, we maintain a consistent number of 400 local search iterations, as employed in our previous experiments. It can be seen that our proposed MLSP algorithms outperform LS algorithm in terms of both clustering cost and running time. On average, the average clustering cost is reduced by 1.15\\% compared with LS using our MLSP algorithm while the running time of our MLSP algorithm is more than 1000 times faster than LS algorithm.\n\n**Table 1: Comparison results of LS, MLS and MLSG on datasets with size smaller than 10,000**\n|Iris(150*4)|Best Cost|Average with Std|Time(s)|Abs_FL(720*21)|Best Cost|Average with Std|Time(s)|\n|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|\n|LS|29.79|30.001$\\pm$0.23|300.61|LS|**1.0786E+06**|1.0873E+06$\\pm$7.8E+03|1513.83|\n|MLSP|**29.74**|**29.761$\\pm$0.03**|0.43|MLSP|**1.0786E+06**|**1.0786E+06$\\pm$0**|0.91|\n|MLS|29.93|30.214$\\pm$0.19|**0.17**|MLS|**1.0789E+06**|1.0970E+06$\\pm$1.1E+04|**0.47**|\n|**SEEDS(210*7)**||||**TR(980*10)**||||\n|LS|214.951|218.9401$\\pm$5.79|436.54|LS|**762.16**|767.88$\\pm$4.1|2049|\n|MLSP|**214.523**|**215.079$\\pm$0.49**|0.63|MLSP|**762.16**|**764.71$\\pm$1.9**|0.94|\n|MLS|214.954|219.294$\\pm$2.82|**0.46**|MLS|785.01|805.89$\\pm$11.5|0.56|\n|**GLASS(214*9)**||||**SGC(1000*21)**||||\n|LS|**251.859**|252.9835$\\pm$0.91|424.28|LS|1.1741E+08|1.2026E+08$\\pm$3.5E+06|2146.91|\n|MLSP|**251.859**|**251.994$\\pm$0.54**|0.52|MLSP|**1.1734E+08**|**1.1752E+08$\\pm$5.1E+04**|9.39|\n|MLS|253.294|254.001$\\pm$0.55|**0.46**|MLS|1.1735E+08|1.1846E+08$\\pm$9.6E+04|**0.62**|\n|**BM(249*6)**||||**HEMI(1955*7)**||||\n|LS|**375974**|378169$\\pm$2739.33|528.74|LS|2.7073E+06|2.7346E+06$\\pm$5.7E+05|3713.27|\n|MLSP|**375974**|**376276$\\pm$265.69**|0.52|MLSP|**2.7070E+06**|**2.7123E+06$\\pm$6.7E+03**|1.08|\n|MLS|378649|384982$\\pm$4131.91|**0.37**|MLS|2.7292E+06|2.7886E+06$\\pm$6.9E+04|**0.59**|\n|**UK(258*5)**||||**pr2392(2392*2)**||||\n|LS|**29.268**|29.602$\\pm$0.34|538.55|LS|**5.3578E+09**|5.4745E+09$\\pm$9.0E+08|4958.77|\n|MLSP|**29.268**|**29.286$\\pm$0.01**|0.61|MLSP|**5.3578E+09**|**5.3668E+09$\\pm$2.3E+07**|1.16|\n|MLS|29.411|30.004$\\pm$0.27|**0.39**|MLS|5.3629E+09|5.4203E+09$\\pm$6.2E+07|**0.58**|\n|**HF(299*12)**||||**TRR(5456*24)**||||\n|LS|**6.9604E+10**|7.2964E+10$\\pm$4.1E+09|634.03|LS|**1.3796E+05**|1.3834E+05$\\pm$566.8|12646|\n|MLSP|**6.9604E+10**|**6.9604E+10$\\pm$0**|0.71|MLSP|**1.3796E+05**|**1.3829E+05$\\pm$253.1**|2.99|\n|MLS|**6.9604E+10**|7.0148E+10$\\pm$4.3E+08|**0.4**|MLS|1.4591E+05|1.4939E+05$\\pm$2.0E+03|**0.64**|\n|**WHO(440*8)**||||**AC(7195*22)**||||\n|LS|**3.3631E+10**|3.3764E+10$\\pm$1.7E+08|913.23|LS|**1163.7**|1169.19$\\pm$8.42|18036.1|\n|MLSP|**3.3631E+10**|**3.3676E+10$\\pm$4.6E+07**|0.89|MLSP|**1163.7**|**1168.2$\\pm$9.01**|1.51|\n|MLS|3.4324E+10|3.5212E+10$\\pm$5.6E+08|**0.4**|MLS|1165.5|1181.7$\\pm$12.8|**0.64**|\n|**HCV(572*12)**||||**rds_cnt(10000*4)**||||\n|LS|1.1312E+06|1.1458E+06$\\pm$4.4E+03|1311.96|LS|1.6104E+06|1.6364E+06$\\pm$2.3E+04|20544.7|\n|MLSP|**1.1311E+06**|**1.1410E+06$\\pm$2.2E+04**|0.9|MLSP|**1.6099E+06**|1.6105E+06$\\pm$7.2E+02|0.88|\n|MLS|1.1505E+06|1.2135E+06$\\pm$4.2E+04|**0.49**|MLS|1.6146E+06|1.6520E+06$\\pm$2.5E+05|**0.42**|\n\n**Question 2: It may be worth explaining what new ideas are needed in the analysis as compared to [10]**\n\t\nResponse: Thanks for the great suggestion. For the LS++ method with single swap strategy in [10], its theoretical guarantee relies heavily on the one-to-one matched swap pairs constructed. When directly applying the ideas in [10] to multi-swap analysis, it has the following challenges: (1) the theoretical bounds on the clustering cost after multi-swaps can not be established; (2) for multi-swap local search, a successful swap involves that data points close to each center of a subset of optimal clustering centers should be simultaneously sampled. \n\nHowever, the LS++ method can only guarantee that a data point close to a single optimal clustering center can be sampled with high probability. To establish the bounds of clustering cost after multi-swaps, we extend the notions of swap pairs to swap sets and propose a new consecutive sampling method to construct candidate centers for swap. To analyze the success probability in each local search step, we propose new structures that divide optimal clusters into different groups for establishing a lower bound of sampling success probability. In summary, new definitions for swap pairs and new structures for analyzing the success probability lower bounds are needed in the analysis as compared to [10].\n\n"}, {"review_id": "739uxVTrfy", "reviewer": "Reviewer_hjYC", "summary": "The authors study the well-known k-means problem: Given a set of points in the Euclidean space, compute k centers such that the sum of squared distances between points and their closest center is minimized. A constant-factor approximation to this problem is achieved by local search: Start with arbitrary k centers and improve the solution by choosing a point which is currently not a center and replace a center in our current solution by this point. If we exchange at most t centers instead of one in the local search method Kanungo et al. showed that the algorithm computes a (3+2/t)^2-approximation.\n\nLattanzi and Sohler proposed a combination of k-means++ sampling and local search (LS++): Start with the k centers computed by the k-means++ solution and improve it via local search. In every step of local search the new center is sampled proportionally to its current cost in the solution. Lattanzi and Sohler showed that their algorithm computes a 509-approximation after O(k*log(log(k))) local search steps. \n\nThe authors extend the algorithm of Lattanzi and Sohler (MLS): Instead of sampling only one point, sample t points S simultaneously and search for a subset T of S such that replacing |T| points in the current solution by T reduces the cost significantly. The main result is the following: After O(k^O(t)*log(epsilon*log(k))) local search steps the algorithm returns a (50(1+1/t)+epsilon)-approximation in expectation. Each local search step takes time O(n*d*k^O(t)), where n is the number of points and d is the dimension of the euclidean space. \n\nIn the experimental part of the paper the authors implemented the MLS-algorithm for t=2 and proposed a second algorithm (MLSP): 1) Use k-means++ for an initial solution 2) For R rounds do the following in each round: Apply the MLS-algorithm, check whether the current solution C is the best observed so far. If yes, again apply MLS to improve it, if not exchange some (or all) centers from C according to a procedure developed by the authors and apply MLS with new center set.\nSo MLSP restarts MLS with a new set of centers if the local optimum which MLS computes is not good. This restart is performed R times. The authors compare the MLS and MLSP algorithm to three algorithms LS++, Lloyd and FLS on several data-sets (50 000- 11 000 000 points) for k=10. In most cases the MLSP algorithm computes the best solution but has a significantly higher run-time than MLS, Lloyd and LS++. ", "questions": "L 69: Should this be n^t instead of n*k^t? If not, please explain. (appears also in L 76, Table 1)  \nL 88: Missing O(...) in the run-time\nL 123: I feel that you should state Lemma 2 as a theorem\nL 136-138: You should mention that this only holds when the cost of the current solution is larger than ((50+1/t)+epsilon)*OPT.\nL 138: we can get an approximation -> we get an approximation\nL 145: In the lower bound on the cost of C there is missing the cost of an optimal solution.\nL 167: Maybe just introduce a map \\phi which assigns every center in C* the nearest center in C? Then you do not need the notation s_c and you can denote the set of optimal centers captured by a center c in C by \\phi^{-1}(c). I feel that the notation should be simplified at a lot of places throughout the paper.\nL  172: find a set \\sigma_h of unused lonely centers -> let \\sigma_h be an arbitrary set of unused lonely centers\nL 173: in the definition of A: don't we have c_h^* in \\phi(s_c_h^*) already by definition as c_h^* is captured by s_c_h^*?\nL 176-177: Please specify how you construct type 2 matched swap pairs. At first it seems that every unused lonely center from C should be in at most one type 2 matched swap pair, but this is not possible. I think what you later need is that every unused lonely center and every c_h^* with |\\phi(s_c_h^*)|>t should form a type 2 matched swap pair? \nL 179-182: For the union of clusters with centers in V you use X(V) while for the union of clusters with centers in Q you use Z(J(Q)), this is very confusing. I suggest to simplify the notation (maybe just directly write the union of clusters instead of introducing a new notation) and also don't change between center sets, clusters sets and union of cluster sets that much in the whole paper when it is avoidable.\nL 212-222: maybe it is better to state this as a lemma\nL 219-220: missing two \" ) \" in this chain of inequalities, also please put this into an align environment as it is hard to read if you include it directly in the text\nL 227:  cost of those optimal clusters -> cost of optimal clusters\nL 230: |Q|>1 : Isn't this automatically the case for sets in H_2? \nL 235: c_h^* can find clustering centers close to it -> c_h^* is close to a center from C\nL242: you use H_1,H_2  etc. are sets of centers, now H_2^* is some set of clusters, this is a little bit confusing\nL267: Shouldn‘t the probability Omega(1/k) also depend on t because of the definition in Line 231?\nL268-269: I belive that this should be independency of events, not the union bound\nL281: comma at beginning of line\nL281: Could you explain why we add Q_S‘‘ and not Q_T to H_G?\nAlgorithm 3 Line 11: Could you explain how this swapping exactly works? Also why don‘t use simply Lloyd here?\n\nAppendix:\nL428: In the following sequence of equations: first inequality remember the reader why s_{o_p} is not in V\nL446: rename one of the (c_h^*,c_h^m) in the argmin equation\nL448: \\mu is previously used as centroid of a set, maybe find an other letter\nL449: In the following sequence of equations, first inequality: should \\mu(L) be \\mu(M_3)? If not please define \\mu(L) \nL452: Please explain this upper bound\nL460: please state the bound for the relaxed triangle inequality as lemma somewhere\nL463: In the following sequence of equations: Z(Q_T) → Z(J(Q_T)), Z(Q‘_S)-Z(J(Q‘_S)), Z(Q)->Z(J(Q))\nL474-476: It is not clear if you mean the failure probability of not sampling at least one point correctly or all points correctly\nL484: In the following inequality: shouldn‘t „t“ be „|v(W)|“ in the exponent?\nL488 In the following inequalities, last inequality: \\Delta(P,C)→ 300\\Delta(P,C)\nL496: „unique lonely center“ misleading as one could think that l(c_h^*) is specific to P_h^*\nL515: \\lambda t^{-k} → \\lambda k^{-t}\nL520-528: this paragraph equals L512-520 → remove", "limitations": "n.a.", "rating": 6, "confidence": 5, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "The extension of LS++ to an algorithm which samples more than one center is natural and the resulting improvement in the approximation factor is significant. Even though some parts of the proof are similar to the approach of Lattanzi, Sohler and Kanungo et al. the proof involves some non-trivial steps and is original in my opinion.\n\nThe experiments show that the algorithm MLSP outperforms other algorithms with respect to the quality of the solution even in the setting where all algorithms run for a fixed amount of time. This may be due to the restart of the local search procedure when it runs into a bad local minimum, so it seems to be a good heuristic to implement in any local search procedure. ", "weaknesses": "Especially the technical part of the paper is very hard to read, which could have been prevented since it is mostly not due to the difficulty of the proofs but because of the complex notation and sometimes imprecise statements. While the result in the main theorem seems to be plausible I am not completely sure about the correctness of the proof.\n\nIt would be nice to have experiments for k>10, since the choice of k<=10 seems not to fit the considered data sizes of 50 000 to 10 000 000 points.  The run-time of MLS and MLSP rises faster with respect to k than the run-time of LS++ so it would be nice to see how their algorithms perform for larger values of k.   \n\nFurthermore the authors outline that the MLS algorithm is relatively fast, especially when compared to its easier version LS++ where only one center is exchanged in every step, which was confusing at first. However the speed-up in the run-time is achieved by an approximate computation of the clustering cost of a solution which could have also been implemented for LS++ for a fair comparison to MLS in my opinion.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the positive rating and the thoughtful comments. In the following we address the concerns.\n\t\n**Q1: Regarding $n^t$ instead of $nk^t$ in L69**\n\t\nResponse: The swap size $O(nk^t)$ used in L69 is actually $O((nk)^t)$, which is a typo. Sorry for the confusion.\n\t\n**Q2: Typos and misleading statements**\n\t\nResponse: We thank the reviewer for pointing these out and we apologize for being unable to respond to each typo error and misleading statement individually due to the character limitations in the rebuttal. In our revised version, we we will carefully correct all the typos and misleading statements in the paper.\n\t\n\t\n**Q3: Maybe just introduce a map $\\phi$ which assigns every center in $C^{*}$ the nearest center in $C$?**\n\t\nResponse: Thanks for the great  suggestion. To make the notations simpler and easier to follow, in our revised version, we will introduce a mapping function $\\phi$ to map every center in $C^*$ to the nearest center in $C$, and use $\\phi^{-1}(c)$ to denote the set of centers captured by $c$.\n\t\n**Q4: Regarding the definition of $A$ in L173**\n\t\nResponse: Thanks for pointing this out. In the construction of $A$, $c_h^*$ already belongs to $\\Psi(s_{c_h^*})$ and $\\Psi(s_{c_h^*}) \\cup \\{c_h^*\\}$ is indeed $\\Psi(s_{c_h^*})$.\n\t\n**Q5: L176-177, please specify how you construct type 2 matched swap pairs**\n\t\nResponse: We are sorry for the confusion caused by the statements on construction of type 2 matched swap pairs. In L176, every unused lonely center should participate in the construction of type 2 swap pairs with every $s_{c_h^*}$ satisfying the condition of $|\\Psi(s_{c_h^*})|>t$. In our revised version, we will make the statement clearer\n\t\n**Q6: Regarding the notations in L179-182**\n\t\nResponse: We thank the reviewer for the kind suggestions. In the revised version, we will make the notation for the union of clusters with centers in $V$ and $Q$ consistent.\n\t\n**Q7: L230, $|Q|>1$. Isn't this automatically the case for sets in $H_2$**\n\t\nResponse: $|Q|>1$ is implicit in the definition of $H_2$. We will remove it in the revised version.\n\t\n**Q8: Regarding the notations in L242**\n\t\nResponse: We are sorry for the confusion caused by the notations. In our revised version, we will make the notations consistent.\n\t\n**Q9: Regarding the probability bound in L267**\n\t\nResponse: We thank the reviewer for pointing this out. Since the swap size $t$ is usually a constant and could be much smaller than $k$, we omit the dependence on $t$ in the probability lower bound. In the revised version, we will add $t$ to the lower bound.\n\t\n**Q10: L281: Could you explain why we add $Q_S''$ and not $Q_T$ to $H_G$**\n\t\nResponse: Thanks for raising this question. $H_G$ is denoted as the collection of optimal clusters whose clustering cost is relatively small with respect to $C$. Recall that $Q_T = Q_S'' \\cup Q_L$. For each $P_h^* \\in Q_L$, with good probability, data points close to $c_h^*$ can be sampled in the $m(P_h^*)$-iteration of the independent $t$ sampling iterations. Thus, $Q_S''$ is added to $H_G$ while $Q_L$ is not added to $H_G$.\n\t\n**Q11: Regarding L11 in Algorithm3**\n\t\nResponse: We thank the reviewer for raising this important question. In L11 of Algorithm3, the swapping works as follows. For each center $c_h \\in C$, the algorithm finds the nearest 50 data points in $P$ to it (denoted as the set $N(c_h)$). For each point $x \\in N(c_h)$, a swap pair $(c_h, x)$ is constructed. If one of the swap pairs constructed can induce a clustering cost reduction, the algorithm will execute the swap. In L11 of Algorithm3, the nearest neighbor searching and the swap pairs construction processes will be repeated until the current clustering reaches a convergence.\n\t\nThe reason that Lloyd's method is not used in L11 is to prevent the algorithm from falling into a poor local optimum too early. Once Lloyd's method is used to adjust the centers, a local optimum is obtained, which reduces the possibility to find potentially better solutions.\n\t\n**Q12: Regarding the notations in L449**\n\t\nResponse: The $\\mu(L)$ term appears in the first inequality actually refers to $L$, which is a typo. Sorry for the confusion.\n\t\n**Q13: L452: Please explain this upper bound**\n\t\nResponse: We thank the reviewer for raising this question. In the following, we give a brief explanation on how to obtain the upper bound in L452, which will appear in the revised version. Let $L1 = \\cup_{c_h^* \\in L}s_{c_h^*}$. By the definition of $\\mu(M_3)$ and $L$, for each $s_{c_h^*} \\in L1$, we can find a set $z(s_{c_h^*}) \\subseteq \\mu(M_3)$ with size $|\\Psi(s_{c_h^*})| - 1$ such that $z(a) \\cap z(b) = \\emptyset$ for any $a, b \\in L1$. For each $s_{c_h^*} \\in L1$, since $|\\Psi(s_{c_h^*})| > t$, $|\\Psi(s_{c_h^*})|/|z(s_{c_h^*})| \\le 1 + 1/t$. Taking a summation for centers in $L$, we have $|L| = \\sum_{c_h^* \\in L}1 = \\sum_{s_{c_h^*} \\in L1} |\\Psi(s_{c_h^*})| \\le \\sum_{s_{c_h^*} \\in L1} |z(s_{c_h^*})|(1+1/t) \\le |\\mu(M_3)|(1+1/t)$.\n\t\n**Q14: Regarding the relaxed triangle inequality in L460**\n\t\nResponse: We thank the reviewer for the nice suggestion. In our revised version, we will include a lemma to state the relaxed triangle inequality in the preliminary section.\n\t\n**Q15: Regarding the failure probability in L474-476**\n\t\nResponse: We thank the reviewer for pointing this out. To make the statement clearer, we will rewrite it as \"The probability of failing to sample a point close to each of the optimal clustering centers\".\n\t\n**Q16: L484 about $|\\nu(W)|$ in the exponent**\n\t\nResponse: We thank the reviewer for pointing this out. The reason why $t$ is in the exponent instead of $|\\nu(W)|$ is based on the following relation between $t$ and $|\\nu(W)|$.\n\nBy the definition of $\\nu(W)$, we have $\\nu(W) = \\lbrace{P_h^*:P_h^* \\in Q''_S\\rbrace}$ and $Q''_S$ is a subset of $Q$. Note that $Q \\in H_2$. By the definition of $H_2$, we have $|Q| \\le t$. Thus, $|\\nu(W)| \\le t$ and $t$ can be used to replace $|\\nu(W)|$ in the exponent of the probability bound."}, {"review_id": "097DIM9ZRf", "reviewer": "Reviewer_bF7j", "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 constant factor approximation and (2) run in time linear in the dataset. In the past literature, one can distinguish essentially 2 types of local-search algorithms for this problem. Single-swap local search which have a simple swap procedure but can only guarantee approximation ratios around 500 at least. Multi-swap local search which use a richer swap structure to give better guarantees (the current best is 9+\\epsilon – approx) but are often impractical since finding a good swap often implies to enumerate over all subsets of some candidate set, which makes the algorithm slower than linear (in the dataset size) time. \n\nThis paper tries to reconcile these two directions by giving a new Multi-swap local search algorithm which gives and approximation factor of 50+\\epsilon and runs in essentially linear time. One of the main ideas in this work is to use sampling techniques in the spirit of the famous k-means++ algorithm and only after this try to swap current centers with the sampled candidates. Because the set of sampled candidates is much smaller than the actual dataset, the running time is significantly improved. They also give experiments that show that their algorithm performs well in practice.\n", "questions": "The authors actually propose 2 algorithms, first MLS for which they prove theoretical guarantees. And a second MLPS which they say is “more practical” but I do not see anywhere in the paper that the same approximation guarantee of 50 holds for MLPS as well. In the experiment sections, MLS seems to perform well, but it is unclear to me if this is the same MLS as the one that was analyzed in section 3 (line 336 the authors say they use the sampling method of MLPS for the MLS algorithm). Does this change the theoretical guarantee of MLS? It would be great of the authors could clarify these points.\n", "limitations": "yes", "rating": 7, "confidence": 3, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "In my opinion, the strengths of the paper are as follows.\n1)\tInteresting result and techniques. Not being an expert, it seems that it is the first time that the sampling technique is used in the context of multi-swap local search.\n2)\tThe algorithm seems to perform well in practice. Experimental section is quite thorough.\n", "weaknesses": "Some minor weakness in my opinion is that there is still a significant gap between the lower bound of 9 for local search and the proven upper bounds.\n\nTypos : \nline 88: benefit\nline 105: much *more*\nline 145: I guess the OPT value is missing in the inequality\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank the reviewer for the positive rating and the thoughtful comments. In the following we address the concerns.\n\t\n**Question 1: The authors actually propose 2 algorithms, first MLS for which they prove theoretical guarantees. And a second MLSP which they say is \"more practical\" but I do not see anywhere in the paper that the same approximation guarantee of 50 holds for MLSP as well. In the experiment sections, MLS seems to perform well, but it is unclear to me if this is the same MLS as the one that was analyzed in section 3 (line 336 the authors say they use the sampling method of MLSP for the MLS algorithm). Does this change the theoretical guarantee of MLS? It would be great of the authors could clarify these points**\n\t\nResponse: We thank the reviewer for raising this important question. The proposed MLSP algorithm is actually a heuristic algorithm for achieving better practical performance in experiments. The main heuristic strategies for MLSP are: 1) random sampling for accelerating the clustering cost updating process; 2) recombination for finding potentially better solutions. As a result, the theoretical bound $50(1+\\frac{1}{t}) + \\epsilon$ might not be applicable to MLSP.\n\nThe proposed theoretical guarantee for MLS algorithm is based on the standard worst-case analysis. However, in practical scenarios, heuristic strategies can help to improve both the clustering quality and runtime of the proposed algorithm. For the MLS algorithm used in experiments, as stated in line 336-337, the sampling strategy is also used to accelerate the clustering cost updating process. Thus, the bound of $50(1+\\frac{1}{t}) + \\epsilon$ might not hold in this case.\n\nIn experiments, it can be seen from table 2 that the running time of MLS and MLSP algorithm can be accelerated by using heurstic strategies such that they can achieve better performance compared with other local search methods, and MLSP performs much better than MLS. \n\nHow to give a theoretical bound on the approximation guarantee for the MLSP algorithm is an interesting problem that deserves further study."}]

Author rebuttal:
We thank the reviewers for their positive feedback and valuable suggestions and we highly appreciate the effort paid by the reviewers to provide in-depth reviews that helped us to improve our work. In the following, we will address the reviewers' comments in detail with separate responses.

Meta-review: This paper proposes a linear-time, local-search based, algorithm for k-means. All reviewers have a positive view of the paper.
I thus recommend it for acceptance.
Decision: Accept (poster)