Fast Local Search Algorithms for Clustering with Adaptive Sampling and Bandit Strategies
Authors: Junyu Huang
Zhen Zhang
Beirong Cui
Jianxin Wang
Qilong Feng
Venue: NeurIPS 2025 poster ()
Keywords: Local Search
Clustering Approximation
Primary area: theory

Abstract: Local search is a powerful clustering technique that provides high-quality solutions with theoretical guarantees. With distance-based sampling strategies, local search methods can achieve constant approximations for clustering with linear running time in data size. Despite their effectiveness, existing algorithms still face scalability issues as they require scanning the entire dataset for iterative center swaps. This typically leads to an O(ndk) running time, where n is the data size, d is the dimension, k is the number of clusters. To further improve the efficiency of local search algorithms, we propose new methods based on adaptive sampling and bandit strategies. Specifically,  adaptive sampling can well approximate the distance-based sampling distribution without maintaining pairwise distances between data points and the centers, enabling fast and accurate sampling in sublinear time after an $\tilde{O}(nd)$ time preprocessing step. The bandit strategy models the best swap pair selection as a bandit problem, where a grouping strategy is proposed for fast identification of the optimal swap pair. With these techniques, our proposed algorithm can achieve constant approximation in expected running time $\tilde{O}(nd + k^4)$ under mild assumptions on optimal clusters and swap pair distributions. Our approach also extends naturally to the k-median objective, achieving constant approximation in expected running time $\tilde{O}(nd + \sqrt{n}k^3)$ without distributional assumptions. Empirical results demonstrate that our algorithm achieves up to 1000× speedup over existing local search methods on datasets with 100 million points, while delivering comparable clustering quality. Compared to coreset-based approaches, it provides up to around 80× speedup and consistently yields better clustering results.

Fast Local Search Algorithms for Clustering with Adaptive Sampling and Bandit Strategies
Junyu Huang', Zhen Zhang?.3, Beirong Cui', Jianxin Wang" Qilong Feng!
School of Computer Science and Engineering, Central South University
Schooo of Advanced Interdisciplinary Studies, Hunan University of Technology and Business
3 *Xiangjiang Laboratory, Changsha, China
The Hunan Provincial Key Lab of Bioinformatics, Central South University, Changsha 410083, China
j junyuuuaggcsu@ffxmmaal com, esuzz@foxmal com
2347110980csu edu. en, jxwangdmail.c csu. edu. cn, csuf ssufeng@mail.. csu. www.cuu.uuc.eeeeeooo cn
csu.
Abstract
Local search is a powerful clustering technique that provides high-quality solutions with theoretical guarantees. With cestance-based sampling strategies, local search methods can achieve constant approximations for clustering with linear running time in data size. Despite their effectiveness, existing algorithms still face scala- bility issues as they require scanning the entire dataset for iterative center swaps. This typically leads an Ondkk) running time, where is the data size, d is the dimension, & is the number of clusters. To further improve the efficiency of local search algorithms, we propose new methods based adaptive sampling and bandit strategies. Specifically, adaptive sampling can well approximate the distance- based sampling distribution without maintaining pairwise distances between data points and the centers, enabling fast and accurate sampling in sublinear time after an ((dd) time preprocessing step. The bandit strategy models the best swap pair selection as a bandit problem, where grouping strategy is proposed for fast identi- fication of the optimal swap pair. With these techniques, our proposed algorithm can achieve constant approximation in expected running time Ö(nd + k*) under mild assumptions on optimal clusters and swap pair distributions. Our approach also extends naturally to the m--mddaan objective, achieving constant approximation in expected running time Ö(nd + Vnk*²) without distributional assumptions. Em- pirical results demonstrate that our algorithm achieves up to 1000x speedup over existing local search methods on datasets with 100 million points, while delivering comparable clustering quality. Compared to coreset-based approaches, it provides up to around 80x speedup and consistently yields better clustering results.
Introduction
Clustering fundamental unsupervised learning problem with wide applications in machine learning. It aims to partition a dataset into clusters such that points within the same cluster can share high similarity. Among various formulations, k-means and k-median are widely studied for their practical effectiveness. This paper focuses on these two problems in d-dimensional Euclidean space.
Over the past few decades, numerous heuristic and approximation algorithms have been developed for the k-means and k-median problems. Among them, the Lloyd's algorithm [27] remains the most widely used in practice. However, a potential issue that Lloyd's algorithm may require an
"orresponding Author
39th Conference on Neural Information Processing Systems (NcurIPS 2025).

---

unbounded number of iterations to reach a convergence M. To address this, Arthur and Vassilvitskii M proposed the k-means++ seeding method, which achieves an O k) approxximato for both k-means and k-median in O(ndk) time, where n is the data size, d is the dimension, and k is the number of clusters. Intuitively, k-means++ selects k centers sequentially, choosing each new center with probability proportional to the (squared) distances from data points to their nearest centers (also known as the D²-Sampling (or D-Sampling) strategy (25)). To further enhance the clustering quality, local search methods are commonly employed to achieve constant approximations 23 25 I 21 77. Starting from a random initial solution, local search aims to find iterative
improvements clustering quality through center swaps between data points and the centers selected. Kanungo et al. (23] proposed the first local search algorithm for k-means, achieving a polynomial- time (9 + <) aapproximation. Arya et al. (2 gave the first local search algorithm for k-median, achieving a polynomial-time (3+ + -aapproximaaionn Recently, Cohen- Addad et al. [13] improved this result for k-median polynomial-time (2.836 -approximatioon by introducing distance- based potential functions, which is the current best approximation result for k-median with local search.
As contemporary datasets scale to hundreds of millions of entries, even algorithms with modest polynomial complexity can become impractical. Hence, in d-dimensional Euclidean space, there has been growing interest in developing local search algorithms with linear runtime in the data size. In this setting, Lattanzi and Sohler (25) proposed the LS++ algorithm. Instead of exhaustive enumerations, LS++ selects exactly one point for swapping in via Da-apliin strategy. They proved that LS++ can achieve an (1))-pprrxiimation for both k-means and k- k-median objectives in time (ndkk² log log k). using O(k log log k) local search steps. Choo et al. [IT further reduced the number of local search steps to (ck)) for LS++, achieving an "aap/""aapooiiimatiio Fan et al. [9 narrows the swap candidates to exactly 2 using greedy strategy, achieving ((1)-pporxi in Ondk log log k) time under mild assumptions on average cluster sizes. While single-swap local search outperforms the seeding methods, there is still gap on clustering quality guarantees compared to multi-swap strategies. To bridge this gap, recent studies have focused on sampling-base -baed multi- swap methods. Huang al. (2 proposed (50(1 + 1/t) + >)approxximatinn algorithm for k-means with 0(ndk2*++1 log(c" log k)) time, where t is the swap size. Independently. Beretta et al. 77] gave a 1..44-aapprximmttio algorithm for k-means in (Idpolyyk:) time, which is the current best result for k-means with linear runtime in data size. These local search methods can naturally be extended to the general metric space to achieve ( (pp-xxopoonnnaoionn and linear running time in the data size.
For clustering in general metric spaces, it has been proved that achieving any constant- feattot ap- proximation requires Ω(nk) time [29]. Therefore, natural question arises:
further improve this bound in d-dimensional Euclidean space? Along this line of research, several algorithms have been proposed. One approach is incorporate the Euclidean spanners (20) with a near-linear time graph-based clustering algorithm 33), where an -(c))xopoximmaion can be achieved in O(n +1/s) time for k-median. However, its extension to the k-means objective remains unclear. Recently, la Tour and Saulpic [24] proposed an ( apprrxxmxiiiinon framework with Ö + n 1 running time for k-clustering, using geometrie hashing and ball-covering strategies for center selections. Alternatively, another direction focuses on coresets constructions [77 16 [4] 12], which aims to
reduce the data sizes to Ō(k/e) with preserved (1 + >)aaproxximation guarantees. While traditional coresets methods require ((ddk) time, Draganov et al. [18] broke this barrier and gave an ((nd) time coreset method. Combined with LS++, this yields an O app)-xpproinnattoon in Ö(nd + k. d) time.
Although some existing algorithms can achieve faster than ((ddk) running times for k-median and k-mmans problems, they appear to be less practical and may still encounter difficulties for handling real-world large-scale datasets. For the almost linear-time algorithms 24 32., while they can offer compelling guarantees on clustering quality and running time, the practical implementations often involve complex data structures. Consequently, it remains unclear whether these methods can effectively scale to datasets with sizes over 100 million. As for coresets methods, even with sampling and spread compression strategies 18) our experiments (see Appendix E.2) show that these methods still face scalability issues. It can take over 2 hours of computations for datasets with 100 million points (while our algorithm can produce higher-quality solutions within just minutes).
In this paper, we aim to bypass the coresets construction methods and develop practical clustering algorithms with theoretical guarantees for large-scale datasets. Given that local search is widely used in large-scale clustering applications, our primary focus is still on local search algorithms.
"The Ô notation is used to hide polylogarithmic factors in and d.

---

Specifically, our goal is to design local search algorithms with running time O(nd + n nnpooy(k)) for a fixed constant Ş E [o, 1). This improves upon the running time of previous local search algorithms in the common practical settings where k is relatively small compared to data size n. In the following, we briefly remark on the challenges encountered in this setting. Firstly, it a non-trivial task to break the ((ndk)) running time for sampling- based swap pair constructions. Existing linear-time local search algorithms typically rely on D2 -a-mmilin [25] 19 7 to avoid exhaustive enumerations. However, this approach requires maintaining distances between data points and their nearest centers, leading to an (nd) update time for each local search step. Since (2()) steps are usually required for theoretical convergence, the total complexity for the sampling process becomes (Iddk). To accelerate the sampling process, several approximation methods have been proposed. including approximate Sampling 3 4. approximate nearest neighbor search [5], and projection based sampling 9 However, these methods are primarily designed for finding initial seedings, where centers sampled sequentially. It remains unclear whether they can effectively handle the setting for local search, in which centers may be frequently inserted and deleted. Secondly, identifying the best swap pair at each local search step is computationally expensive, evaluating the clustering cost change for even a constant number of candidate swaps requires scanning the entire dataset, leading to an (ndkk) complexity. To accelerate the best swap pair identification, recent approaches proposed to formulate swap pair selection as a bandit problem 34 33 [ Nevertheless, these methods still require a search space that scales linearly with the data size n, resulting in a runtime complexity of (ndk)). Moreover, their theoretical guarantees rely on restrictive assumptions, requiring the cost changes induced by each swap pair to follow Gaussian distribution N/7,, 1) with a constant mean y.
1.1 Our Contribution
The main objective in this paper is to develop faster local search algorithms for clustering while bypassing coresets construction methods. For k-means clustering, we first propose a sampling-based algorithm achieving constant approximation in expected running time Ö(nd / vnk**) assuming uub-Gaussian prior on cost changes of swap pairs. To accelerate the sampling process, we propose an adaptive sampling method that dynamically maintains tree structure to approximate sampling distributions. The tree structure can be initialized in ((nd) time and updated within sublinear time. Unlike previous tree embedding methods (i.e., HSTs [15] 24 ), our approach dynamically maintains weight distribution for each internal node, enabling an approximate D Sampling process to be executed on the tree. For the best swap pair identification, we model the task as a bandit problem, regarding swap pairs arms. Existing bandit methods usually require strict Gaussian assumptions for cost changes of each swap pair [6
with constant mean and deviations. To relax this assumption, we propose combine sampling with swap pair grouping strategy, which can eliminate unnecessary swap pairs assuming sub- Gabbaaassinn prior on cost changes of swap pairs with bounded deviations. With theses techniques, the best swap pair can be determined in subliner time Ö(k*)).
34 33
To further enhance the scalability of local search algorithms, we propose a more practical local Hastings strategy for approximating the D2 sampling process. Metropolis- strategy generates candidate samples from uniform proposal distribution and employs rejection sampling to determine transitions between samples. We show that our algorithm can effectively approximate the D2-Sampling distribution under mild assumptions regarding the sizes and tail behavior of optimal clusters, which enables each sampling step to be executed in ((/kk) time. By combining with bandit methods, a constant approximate solution can be obtained in expected time Ö(nd + k:4). Experiments show that our algorithm can achieve clustering quality comparable to existing local search methods, with up to 1000x speedup on datasets containing 100 million points. Compared to coresets, it provides up to around 80x speedup with better qualities.
search algorithm that adapts Metropolis- Metropolis-Hastingg
Finally, extend our adaptive sampling approach to the k-median objective by introducing weighted sampling strategy. This strategy can estimate the k-median costs for swap pairs using weighted samples drawn from the dynamic tree structure maintained. We show that clustering costs can be approximated within (1 + e) factor using Ö(Vr/c) samples. By incorporating swap pair grouping strategy, sub-optimal swap pairs can be effectively filtered in Ó(vnk?/ /c) time for each local search step. Thus, our approach can give a constant approximation for the k-median objective in expected Ö(nd + vūke²) time, without requiring any distribution assumptions on data points.
Putting all these together, we can get the following results.
3

---

Theorem 1.
With sub- Gaussian prior on cost changes of swap pairs, there exists an algorithm for k-means that can achieve constant approximation within expected running time O(nd + vrk**)
Theorem 2. With sub- Gaussian prior on cost changes of swap pairs and mild assumptions on optimal clusters, there is constant approximation algorithm for k-means with expected runtime Önd + k*).
Theorem For the k-median problem, there exists an algorithm that can output a constant approximate solution in expected O(nd + vrkk²) running time without any data assumptions.
Table 1: Comparisons with existing local search methods (n: data size, d: dimension, t: swap size).
Objective
Ref
Appeoximation
Assumption
Running Time
9+
dpoly(is, A))
25
O(1)
O Olddk log k)
0(1)
O(rul(2)
19
0(1)
Average Cluster Sizes
Onddk log log k)
k-means
50(1 + 1/4) +
O Qadkat++ logie log k))
10:46
((ndpoly(8))
O(1)
Suh Gaussian Prior
Ö(vd v
This Paper
O(1)
Sub Gaussian Prior, Tailed
Ö(ud+ + k4)
Behavior and Sizes of Optimal Clusters
3++
- dpolyln, R) Olndk log log k) O(udk"))
O(1)
11]
((1)
m-median
((1)
Average Cluster Sizes
Oladk log logg k))
M1
O(1) (<c 50(1 + 1/8) +
O (ndk 2++1 logle" -1 log R))
O(1) (<. 10.46)
((ddpoly(8)
Paper
O(1)
O(nd rik-1)
Table presents detailed comparisons with state-of-the-art results, For the k-means objective, our algorithm achieves constant approximation while reducing the joint dependence on and k, improving the efficiency when k < n. For the k-median objective, our method also ensures constant approximation with faster running time without relying on any distributional assumptions.
2 Preliminary
In this paper, we use P C Rdd and k to denote the set of the given dataset and the number of clusters to be opened, respectively. For any two data point p. Q Rd let 5(p, q) and 5² (p. q) denote their Euclidean and squared Euclidean distance, respectively. Given a set C' C Rd of data points, denote 8(P,C)= = Epep minlee 8(p, c) and (P, = Epep mincec 82 (p e) as the clustering cost of P with respect C' for the k-meddia and k-means objectives, respectively. We use OPT to denote the optimal clustering cost of the given clustering instance. The goal of the k-median (k-means) clustering is to find a set C C Rd with size while minimizing 8(P, C²) (or (P, C') for k-means).
The dimension of given clustering instance can be reduced from d to Ollog n) via the standard Transform, which incurs only constant factor approximation loss 22 Without loss of generality, we can assume that d = Ollo n) through an Ó(nd) time preprocessing step for dimensionality reduction.
3 The BanditLS Algorithm
In this section, we present the BanditLS algorithm, which integrates adaptive sampling and bandit- based strategies to accelerate local search for k-means. Our main objective is to design a local search algorithm with running time Ö(nd + /npoly(k)), thereby reducing the combined impact of parameters n and i on the running time. This can improve the efficiency of previous local search methods, especially when & is relatively small compared to the data size n. While coreset-based approaches can achieve Ö(nd + k-dd) time, they often fail to scale well in practice due to time overheads. For datasets with sizes over 100 million, the constructions for coresets can require several hours of computation, making it less practical for large-scale scenarios (see Appendix E.2 for details). The key challenges here lie in bypassing the coresets construction methods while accelerating both the sampling process and identification of the best swap pairs during the local search swaps.

---

The formal description for the proposed BanditLS algorithm is given in Algorithm We first outline the main idea behind. To eliminate the ((ndk) overhead caused by explicitly maintaining pairwise distances for D Samplinnn we propose an adaptive sampling method. Specifically, we first embed the dataset into tree structure in near-linear preprocessing time, aiming to discretize the candidate distances between points and centers from I: down to O (nd)). Instead of computing pairwise distances, we then dynamically maintain weight distributions at internal nodes, capturing the distance distribution between data points and their centers. Based on traversals of the maintained weight distributions from root to the leaf, the D'-Sampling process can be approximated within sublinear time. To efficiently select the best swap pair without exhaustive scans during the swaps, we formulate the task as bandit problem. Unlike previous methods relying on strict Gaussian assumptions for all swap pairs, we propose a sampling-basse based grouping strategy that only requires assumptions for high-impact swap pairs. This approach can filter out ineffective swaps using only Ö(k²) samples, reducing the swap selection complexity from O(ndk) to be independent of the data size
Starting from a random solution, BanditLS proceeds in three main stages: (1) Tree construction (steps 1-5) where the data is embedded into a counded-height tree T' and centers are opened or closed by marking the corresponding leaf nodes; (2) Swap pair construction (steps 7-9) where an adaptive sampling strategy is used to approximate the D2-Sampling distribution via dynamic weight distributions maintained, allowing efficient sampling without explicitly tracking nearest-center distances; (3) Swap pair identification and tree updates (steps 10-17), where the best swap pair is selected via bandit model. To avoid restrictive assumptions, a grouping strategy is proposed to include virtual swap pairs with zero cost changes, enabling sublinear time identification of promising swap pairs under sub- Gaussin prior. These components together yield the improved runtime.
Algorithm 1 BanditLS( k, C,n,o)
Input: A k-means clustering instance (P, k; d). set C' C Rd of random initial clustering centers, parameter 7) E (0, 1). a parameter o >> representing the prior knowledge on deviations. Output: set C RRd of clustering centers.
1: T = crre-Construction P, k, d). & Coee-Construttin is detailed in Algorithm 4 in Appendix A Poee--onversion((7
2: T'
Tree-Conversion is detailed in Algorithm 6 in Appendix A
3: Call FIRE-INII(T") to initialize the dynamic data structure.
> Algorithm in Appendix
4: for E C do
5: Call the TREE-OPEN(T algorithm to mark the leaf node in ''' associated with ess active
and update the tree structure T TTREE-OPEN is detailed in Algorithm in Appendix A
6: for = 1 to Ö(vnk) do
7:
æ
C). Adaptive-Samplinn is in Algorithm 10 in Appendix A
8:
Set 8 arg minnge 0(x, c) and randomly sample center o' E ()(2)))
9:
S {(x, 8x), (82, 82 82)). S2 = {(x, d"), (4',
01 BanditCS
d, S, σ). 02 BandiiCS(PP k. C,
> Algorithm
11: if Joı1 > 1 and |02| > 1 then
Continue.
13:
else
Let O for 1,i E (1, and randomly choose an from C- fu), where (u, v) E
15:
Call the TREE-CLES CHEE-CLOSE(T' v) algorithm deactivate the leaf node T associated with
and update the tree structure T'. TTREE-COSE is detailed in Algorithm 9 in Appendix A
17:
Call the PREE-OPEN(T', u) algorithm to mark the leaf node T associated with u as active and update the tree structure T'.
18:
return C.
3.1 Tree Construction and Adaptive Sampling
In this subsection, we give the analysis for the proposed tree construction and adaptive sampling methods, which are the key components for accelerating the D² process. The core idea is to discretize the distances between data points and their centers by embedding them into a dynamically maintained tree structure. This reduces the number of distinct distance categories from k to Oloog (nd)), corresponding directly the tree's 's height. Thus, by traversing from root to the
5

---

leaves, the D2 Sampling process can be approximated within (log (nd)) time. Due to space limit, we present the main ideas here and leave the detailed algorithmic descriptions in Appendix A.
Our approach begins by embedding the clustering instance into a tree structure T (see Algorithm 4 in Appendix A. for details) with height H (llog(d))) where A is the aspect ratio of the clustering instance Each node E T is assigned level (e()), where leaf nodes have level 1, and the levels increase by one each node when moving upward toward the root (see step 7 of Algorithm For any two points P; q E P. let Up.o denote their lowest common ancestor in T. The tree distance between p and q is defined as (p, q) = 2vd (2/o(vp.a) - 2). Given set C' of centers, define (p, C) = min.ec ôT (p, c) as the tree distance from data point p to the center in C that shares the lowest common ancestor with p in T. If each internal node in the tree T has at most two children, let I(v) and r(u) denote its left and right children, respectively. If v has only one child, the child of is treated its left child for the ease of analysis. Given node of the tree T, let B, denote the set of points corresponding to the leaf nodes that are descendants of node v. Denote ôt (B., C) = En'gs, ô7(v', C). The following lemma shows that the tree structure T can well approximate the pairwise distances.
Lemma 1. (Cohen-Addad et al. [B ( s (p.9) and El8T(p. 4)) S O(d) 8(p.9).
Based on the tree structure, we will maintain weight distributions at each node of the tree T. enabling a simulated D2 D'-Sampling process by traversing from root to the leaves. However, since the internal nodes may have (nn) children, sampling could require traversing ((n) branches, resulting in Qn log(d)) running time. To address this issue, a tree conversion algorithm (Algorithm 6 in Appendix A.2) is proposed to convert the tree structure T into new tree structure T' with bounded branches by creating auxiliary nodes. The following lemma shows that this process only adds an Olog(dA)) factor to the tree height while preserving the bounds for the expected pairwise distances.
Lemma 2. Let T' be the modified tree structure of the output for Algorithm Then, each vertex in T' has at most 2 children. and the height of the tree T' can be bounded by Olog"(dA)) Additionally, it holds that 8(p, g) S (p, q) and Elōy- 4) < ((d) 5(p, g).
Next, we describe how to maintain a data structure 'D to approximate the D2- -amppnggn distributions. Let C' be the set of centers dynamically maintained by the local search process. The data structure D on T' should satisfy the following properties.
Each non-auxiliary node v in the tree T' is assigned a level le(v), with leaf nodes assigned with level 1. Along any path from leaf to the root, the level increases by at each non-auxiliary node, while auxiliary nodes retain the same level as their parents.
Each node e is associated with three values n²,- nt, and n(v). For non-leaf nodes, ni, and m. denote the number of leaf nodes in the left and right subtrees of v, respectively, n(v) represents the number of active leaf nodes (or opened centers) that are descendants of v. Each node u E T is also associated with a weight w(e) where wwv) = 8T (B., C*).
Each lonfleaf node u E ' is associated with a probability list p(v), which represents the proportion of the total weights contributed by each of its children. If has only one child, then p(v) [1, 0]. 6y (B)co).C)
If v has two children, p(v) = [e], where Lg
(B-ri)))) and 4 + To = 1. (B,.C)
87, (B..C)
The data structure D can be initialized via a bottom-up traversal from the leaf nodes to the root (see Appendix A.3 for details). However, center swaps during the local search process may influence the properties of the maintained data structure. To address this issue, we propose two operations: TREE-OPEN and TREE-CLOSE (as detailed in Appendix A.3). Specifically, TREE-OPEN marks a leaf node as active and updates D along the path from the leaf to the root, while TREE-CLOSE deactivates a marked leaf node and also updates the structure along the path from the leaf to the root. The following lemma shows that both operations can achieve sublinear update time.
in time (dA)) to update the tree structure. Let C' be the set of centers after a swap. Then, each non-leaf node u E T holds a probability list p(u) = Ul..rol where ly
Lemma 3. The TREE-OPEN and CREE-CLOSE operations
57 C") and re (B.(B. C')
by (B.,C")
(B,)
A
max p. 8(pig) ' "...o
and can be compressed to poly(n, d) in time Ö(dd) (Draganov et al. [8]

---

Using the TREE-OPEN and TREE-CLOSE operations, we can always maintain tree structure T' satisfying the properties of data structure D. The adaative-Sampling algorithm (Algorithm I0 Appendix A.4) then performs a sampling process on 'I' to approximate the D2 Di-Sampling distribution: initializes an empty set S and repeatedly performs sampling on the tree TT' until point accepted. Each sampling step starts at the root and proceeds down the tree, choosing children based on the probability lists. Upon reaching a leaf, the corresponding point is selected using a rejection sampling rule. The following lemma shows that D'-Sampling can be approximated within sublinear time.
Lemma 4. The Adaptive-Sampiln Sampling Algorithm takes expected ((k) time to sample a data point x: with probability 91 ((,,C))
(P.)
3.2 Bandit-Based Center Swap
In this subsection, we present the bandit-based center swap method, which aims to accelerate the best swap pair selection process. Specifically, we reformulate the computation for clustering costs as best-arm identification problem, where each swap pair corresponds to an arm in a multi-armed bandit framework. In this setting, we propose grouping strategy to estimate the clustering costs. Due to space limits, we leave the intuitive ideas here and leave the detailed analysis Appendix B.
The proposed algorithm is given in Algorithm which dynamically maintains confidence intervals for the clustering costs of swap pairs via successive sampling (steps 7-8). As the sample sizes increase, suboptimal swap pairs can gradually be eliminated. To achieve sublinear sample sizes, prior bandit methods require strict Gaussian distribution assumptions on each swap pair 6 34 [ with constant mean and deviations (see the proofs in Appendix I of (6] for details). To relax this, we propose to include virtual swap pairs with zero cost changes and group them with other swap pairs. This enables the elimination of suboptimal swap pairs after sampling Ó(k2) points, under a weaker assumption that the swap pairs follow a sub-Gaussian distribution with bounded deviations.
Algorithm 2 BanditCS k. C,S,
Input: k-means clustering instance (P, ki d), a set C Ied of clustering centers, swap group S, probability parameter 0 < 7 < 1, and an >> 0 representing the prior knowledge on deviations. Output: group Sawap of swap pairs.
1: Swwap = S, Spe P. niet = 0, B =
2: Set jı, 0 for
(u, v) E S.
3: while nher < Ō(k2) and [Sswap > 1 do
4: Sample a data point 6 randomly and uniformly from P, and set B = (b).
5: For each s (u, v) Sswap. let Gb) 8 (b, (V( {u)) (b,
6: for s E Swwap do
p/g Aad HG.
7:
Peet +1
8: F. F.to
2 2m(1/n)
9: Sewap = {s E Sawap Po - F. S minyesma (By + Fu)}, Dyef + Thref + 1.
10:
return Swap.
Let C' be the set of the centers opened. Given data point b randomly and uniformly sampled from the dataset P and swap pair 8 = (u' v²) (u' E P, o' E C), the clustering cost changes induced by tu',"" on b is defined C, = 82 (b, C' U U {u')) (b, C). Under the following assump- tions on swap pairs with sub- Gaussian distribution and bounded deviations, with high probability, Algorithm can return a swap pair with minimum clustering cost using ((k2) samples. Putting all these together, Theorem can be proved (the proofs is given in Appendix B).
Assumption 1. Let C be the set of centers before executing each step 10 of Algorithm ! For each swap pair s E St U S2, it is assumed that G.(b) is o² sub-Gaussian,fer for a randomly sampled point b E P. with a known bound o² O((62 (P, ())) and o² 2 o².
Lemma 5. Let S be set of swap pairs as input for Algorithm 2 which contains at least one swap pair with (1 1/100k) fraction of cost reduction. For n O(n-t) Algorithm 2 can remove non-promising swap pairs in S with high probability using Ō(k2) samples.
7

---

In practice (also in our experiments), since is usually unknown, we estimate for each swap pair as the standard deviation xx STD (TDyes(G,(y)) which can be computed from the of sampled points in B accordingly.
4 A More Practical Algorithm for k-means Clustering
In this section, we propose more practical local search algorithm for large-scale clustering scenarios. The proposed algorithm is presented in Algorithm B where the intuitive idea behind is to combine our bandit-band baadit-based strategy with Metropolis-Hastie Due space limit, present the high-level idea here and leave the analysis in Appendix C.
technique to accelerate the sampling process.
Algorithm 3 BandifausLSS(Pkk,d..C
Input: A k-means clustering instance (P, k, d), set C' C E of random initial clustering centers, parameter 7) E (0, 1), a parameter a representing the prior knowledge on deviations. Output: set C Rd of clustering centers.
1: for 1 to ((k) do
2: Randomly sample a data point æ E P.
3: i to (m//k²) do Randomly sample data point y E P. with probability min {I, 8 (..)) (a..C)
4:
5:
Set ret +
6: $. = arg mingec 0(),z).
7:
Randomly sample a point o' from C'{sz),
8: S = (z. (8,, 8z) 1, Sg = [(z d)), (d'
9:
0n BanditCS k d, C,),,,,,), 02 BanditCS k. C, 1,))
10: if |oi > 1 and [02| > then
11:
Continue.
12:
else
Let O = (o: 10 L, [1,2}), and randomly choose swap pair o' from O. C= C Uu, where (u, v) E o'.
14:
return C.
Starting with a random initialization, the BanditFastLS algorithm consists of the following two stages: (1) approximate sampling (steps 2-6); (2) bandit-based swap pair identification (steps 7-14). In the approximate sampling stage, a Hatropolis-Hasting strategy is used to approximate the D²-Sampling distribution, where a Markov Chain of bounded length (n/kk²) is constructed. This reduces the rumning time for D -a-mmmlin by eliminating the need to maintain distances between data points and their closest centers. In the bandit-based swap pair identification stage, the bandit- method proposed in Section 3.2 is adapted to accelerate the swap pair selection process.
Under mild assumptions on optimal cluster sizes and tailed behavior of optimal clusters, we show that the D Sampling distribution can be well approximated within (nnkk² rounds of sample transitions, where an improved rumning time can be achieved for swap pair construction. Putting all these together, Theorem can be proved (the proofs are given Appendix C).
Assumption 2. Assume that the given dataset P is average where each optimal cluster has size (2kk2). For each optimal cluster PN, let eî, be the optimal center for Pž. We assume that each Pž follows distribution F over Rdd with exponential tails, i.e., 3c, such that (, (4) > a] S ce holds for I P* where c, f are constants and μ is the mean of the distribution.
Lemma 6. After Ö(n/k²) Metropolis- Hastings sampling steps (steps 2-5 of Algorithm 3). if the given clustering instance satisfies the properties Assumption 2, we can sample a data point æ e P with probability at least 0.582 ())
5 Extension to the k-median Objective
Our adaptive sampling and bandit-based methods naturally extend to the k-median objective. To climinate the assumptions and sub-Gaussian prior from the bandit process, propose weighted

---

sampling method that can estimate the k-median cost of each swap pair using a small number of samples taken from the dynamic tree structure (see Algorithm 12 Appendix for details). This method achieves a (1 = <) estimation error with only Ö(vn/c) samples. By adapting key ideas from bandit methods, we show that a constant-factor approximation can be obtained in Ö(nd + vnk*) time for k-median, without relying on any distributional assumptions. Due to space limit, we leave the detailed algorithmic description and analysis in Appendix D. Putting all these together, Theorem 3 can be proved (the proof's are given in Appendix D).
Experiments
In this section, we compare our proposed BanditFastLS algorithm (the proposed practical algorithm in Section 4) with other state-of-the-ar local search methods and coresets methods. For hardware, all the experiments are conducted on machine with 100 Intel Xeon Gold 6348 CPUs and ITB memory.
Datasets. We evaluate our algorithm on both small and large-scale datasets from prior k-means studies [21 25). Large-scale datasets include SYN (IM x 2). USC_ CCCC__9900__0 (2,45M x 68). SUSY (5M x 17) and HIGGS (11M x 27) from UCI Machine Learning Repository 4 We also include a dataset SIFT (100M x 128)
30 Small datasets are with sizes ranging from 150 to 50,000 used in [2]
Algorithms. We mainly compare our algorithm with the following: LSDS++ 19 (the fastest single-swap local search with constant approximation), LS++ [25], MLS (21) (a fast version of multi-swap local search), and BanditPAM++ (33] (a near-linear time bandit-based algorithm). For our algorithm (BanditFastL.S. we set 7) = 1/(2n") ) to satisfy theoretical guarantees. The number of Metropolis- Metropolis-Hastinnss steps is fixed at 20, and we limit the maximum bandit samples to 50,000.
Experimental Setup. We evaluate our BanditFastLS algorithm against other local search algorithms, following the setup in [25] 19. [1]. For fair comparison, the initial centers are randomly selected instead of using a seeding method. For datasets with fewer than 50,000 points, we set k ranging from 3 to 10. For larger datasets, we set k ranging from 10 to 100. Each algorithm is executed for 10 times, and we report the average results with deviations and average runtime. To ensure faireess, all algorithms perform 400 local search steps. Following 225). we apply 10 iterations of faster Lloyd's algorithm (mini-batch Lloyd's 's algorithm (3I) after each local search algorithm to finalize the centers. For large-scale datasets (such as HIGGS and SIFT), since LS++'s 's nearest-neighhor index is cemmory-rooibiiiive we implement it via direct pairwise distance computations.
Results. Table 2 and Table 3 (Appendix E) present comparison results large-scale datasets with more than 1 million points. On the SIFT dataset with 100 million points, even with moderate number of clusters to be opened (i.e., k 10), only MLS, LSDS++, and BanditFastl.S S algorithms can return feasible solutions within 24 hours. As the parameter l increases, it takes over 24-48 hours for MLS, LSDS++, and LS++ algorithms return a feasible solution. In contrast, our BanditFastLS algorithm scales efficiently, producing comparable (or even better) clustering results within minutes. Specifically, on dataset SIFT with A: 10, BanditFastLS is over 200x faster than competing methods. With h = 100, it achieves up to 1000x speedup. By calculating the average results across all scale datasets, BanditFastLS is approximately 210x faster than other local search methods, with only a 2% increase in clustering costs. When comparing with the baseline algorithm mini- bati- k-means++, it also achieves better clustering quality and runtime on dataset SIFT for different choices of parameter k: Table 4 (Appendix E) presents results on smaller datasets (datasets with fewer than 50,000 points), where all local search methods perform well, and MI. and LS++ slightly outperform others in clustering quality. However, as data sizes and k increase, BanditFastl. consistently delivers faster runtimes while maintaining comparable clustering performance.
We also evaluate the impact of varying parameters and compare our method with coresets in Appendix E. The input parameters for the BanditFastLSS algorithm mainly include the probability parameter and the number of Metropolis- Hastings steps used to approximate the D2-Sampling distribution. Table 5 presents the results for datasets SYN, USC 1990, and SUSY. with a fixed number of clusters k = 50 and a maximum sample size of 50,000, while the Metropolis- Hatstigs steps vary from 20 to 100. The results indicate that increasing the number of Setropolis-Hastingg steps can improve the overall clustering quality, with only slight increase in the running time. However, the overall clustering performance remains consistent, demonstrating that the algorithm is robust to variations in
9

---

Table 2: Results on the HIGGS (11 million) and SIFT (100 million) datasets with varying k, where BanditPAM excluded as it failed to return solution within 24 hours in all cases.
Method
Dataset
k
Cust
Time(s)
Dataset
h
Cost
Time(s)
k-means++
1.6619E488
1.3E+06
6.51
1.1160=41322..4--100
181.03
MLSS
1,5651E+48
5.6F+06
451.91
1. 1.07077+1332.20--1
24286.19
LSDS++
HIGGS
10
1.5618E+08 2 2,9E405
1162.77
SIFT
10
1.07006433:1.00--11
19994.26
LS++
1.5704E+08= -
6902.73
2-44
Ours
5.73
1.6411E+11::3,8:::110
99.81
k-means++
331438+08=1.3E444
9.25
1.094B+13311.76--00
190.74
MLS
1.4273E+08- 1.2E+06
756.35
24h
LSDS++
HIGGS 20
1,4247E+08-- -7.2E+00
1432.91
SIFT
20
>8h
LS++
1.4321E+08_-77
-.72440
18845.10
244
Ours
1.4335F+08- 5.21+005
5.73
12
101.33
1.4296E+008
1.1E+06
13.18
9.6697E=+1::152E5-09
276.36
MLS
1136.02
4-8h
LSDS++
HIGGS 30
1.350GE+0S=11.B44
1883.43 43
SIFT 30
2-4
LS++
1.5501E+0=11.0---
35030.13
48h
Ours
1.3527E+08 4.7E+005
6.92
9.2386E+12 1,47++
105.85
1.3119E+08-1
20.58
9.1771E412.1.122-099
381.74
MLS
1.2658F=008 1.++00
1658.98
4-8h
LSDS++
HIGGS 50
1.2638E+08-2.9EN
2227.33
SIFT
I.S++
>-24h
3-8h
Ours
1.2636E00-333888 3..3E+00
7.88
8.7442E+12-+11060
125.89
k-means++ MLS
1998F 00004088-666.000 95
40,98
1.47745412-5.44442
677.93
1. 141444+00- 1.0E+05
3181 .7
2-8h
LSDS++
HIGGS 100
1.1438E+08
3676.71
SIFT
100
LS++
2-48h
> --88
1.1434HH08
1.3++05
the number of Metropolis Hatropolis-Hastingss steps. This indicates that fixing the number of transition steps 20 sufficient to achieve high- -uali-u clustering results with fast running time.
Tables 7 and 8 present comparison between coreset construction time and the final clustering time for our algorithm on large-scale datasets. The results show that our algorithm outperforms coreset methods, achieving 20x speedup average compared to state-of-the-att coreset construction techniques. In general, the experimental results show that our algorithm is robust across different parameter settings and consistently outperforms coreset-based methods.
7 Conclusion, Broader Impact Discussion and Limitations
In this paper, we propose fast and practical clustering algorithms leveraging adaptive sampling and bandit methods. Under mild statistical assumptions data distributions, the proposed algorithms can achieve better scalability while maintaining comparable clustering quality. Experiments show that, compared with the state-of-the-art local search methods, our proposed algorithm can handle datasets with sizes over 100 million while achieving up to 1000x speedup.
The main purpose of our work is to provide algorithmic insights to accelerate large-scale clustering tasks, with no foresceable negative socictal impacts, One potential limitation is the statistical assumptions. Such assumptions unknown. Hence, designing fast local search algorithms without any data distribution assumptions is an interesting future direction that deserves further studies. Another potential limitation is the trade-off between clustering quality and efficiency. While our proposed method reduces the combined influence of parameter A: and data size n the runtime complexities, may not scale well for sufficiently large k (i.e., (.))).. This motivates future work on local search algorithms design with strong approximation guarantees and running time independent of the number of clusters kt
difficult to validate in practice, as the optimal clusters usually
Acknowledgments and Disclosure of Funding
This work was supported by National Natural Science Foundation of China (62502545, 62432016), the Science and Technology Innovation Program Hunan Province (2025RC3207), Open Project of Xiangjiang Laboratory (25XJ03009), and Central South University Research Program of Advanced Interdisciplinary Studies (2023QYJC023). This work was also carried out in part using computing resources at the High Performance Computing Center of Central South University.
10

---

References
[1] David Arthur and Sergei Vassilvitskii. k-means++: The advantages of careful seeding. In Proceedings of the 18th Annual ACM-SIAM Symposium Discrete Algorithms, pages 1027- 1035, 2007.
[2] Vijay Arya, Naveen Garg, Rohit Khandekar, Adam Meyerson, Kamesh Munagala, and Vinayaka Pandit. Local search heuristic for k-median and facility location problems. In Proceedings of the 33rd Annual ACM Symposium on Theory of Computing, pages 21-29, 2001.
[3] Olivier Bachem, Mario Lucic, S Hamed Hassani, and Andreas Krause. Approximate k-means+ in sublinear time. In Proceedings of the 30th AAAI Conference on Artificial Intelligence, pages 1459-1467, 2016.
[4] Olivier Bachem, Mario Lucie, S Hamed Hassani, and Andreas Krause. Fast and provably good seedings for k-means. Proceedings of the 30th International Conference on Neural Information Processing Systems, pages 55-63, 2016.
[5] Olivier Bachem, Mario Lucic, and Andreas Krause. Scalable k-means clustering via lightweight coresets. In Proceedings of the 24th ACM SIGKDD International Conference on Knowledge Discovery & Data Mining, pages 1119-1127, 2018.
[6] Vivek Bagaria, Govinda Kamath, Vasilis Ntranos, Martin Zhang, and David Tse. Medoids in almost-linear time via multi-armed bandits. In Proceedings of the 21st International Conference on Artificial Intelligence and Statistics, pages 500-509, 2018.
[7] Lorenzo Beretta, Vincent Cohen-Addas Silvio Lattanzi, and Nikos Parotsidis. Multi-swap k-meanstt. In Proceedings of the 37th International Conference on Neural Information Processing Systems, pages 166669926099 2023.
[8] Lorenzo Beretta and Jakub Tětek. Better sum estimation via weighted sampling- In Proceedings of the 33rd Annual ACM-SIAM Symposium Discrete Algorithms, pages 2303-2338, 2022.
[9] Moses Charikar, Monika Henzinger. Lunjia Hu, Maximilian Vötsch, and Erik Waingarten. Simple, scalable and effective clustering via -ie-dimeniional projections. In Proceedings of the 37th International Conference on Neural Information Processing Systems, pages 64618-64649, 2023.
[10] Siddhartha Chib and Edward Greenberg. Understanding the metroppolishhastings algorithm. The American Statistician, 49(4):327-335, 1995.
[11] 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, pages 1909-1917, 2020.
[12] Vincent Coden-Addad, Kasper Green Larsen, David Saulpic, Chris Schwiegelshohn, and Omar Ali Sheikh-Omar. Improved coresets for euclidean k-means. In Proceedingy of the 36th Annual Conference on Neural Information Processing Systems, pages 2679-2694, 2022.
[13] Vincent Addad, Anupam Gupta, Lunjia Hu, Hoon Oh, and David Saulpic. An improved local search algorithm for k-median. In Proceedings of the 33rd Annual ACM-SIAM Symposium on Discrete Algorithms, pages 1516-1612, 2022.
[14] Vincent Addan, Kasper Green Larsen, David Saulpic, and Chris Schwiegelshohn. To- wards optimal lower bounds for k-meddian and k-means coresets. In Proceedings of the 54th Anmual ACM SIGACT Symposium on Theory of Computing, pages 1038-1051, 2022.
[15] Vincent Addad, Silvio Lattanzi, Ashkan Norouzi-Fard, Christian Sohler, and Ola Svens- son. Fast and accurate k-means++ via rejection sampling. In Proceedings of the 34th Interna- tional Conference on Neural Information Processing Systems, pages 16235-16245, 2020.
[16] Vincent Cohen-Addad, David Saulpic, and Chris Schwiegelshohn. A new coreset framework for clustering. In Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory of Computing, pages 169-182, 2021.

---

[17] Yichuan Deng, Zhao Song. Yitan Wang, and Yuanyuan Yang. A nearly optimal size coreset algorithm with nearly linear time. arXiv preprint arXiv:2210.08361,, 2022.
[18] Andrew Draganov, David Saulpic, and Chris Schwiegelshohn. Settling time vs. accuracy tradeoff's for clustering big data. Proceedings of the ACM on Management of Data, 2(3):1-25, 2024.
[19] Chenglin Fan, Ping Li, and Xiaoyun Li. Lsds++: Dual sampling for accelerated k-means++. In Proceedigns of the 40th International Conference on Machine Learning, pages 9640-9649, 2023.
[20] Sariel Har-Peled, Piotr Indyk, and Anastasios Sidiropoulos. Euclidean spanners in high dimen- sions. In Proceedings of the 24th Annual ACM-SIAM Symposium on Discrete Algorithms, pages 804-809, 2013.
[21] Junyu Huang, Qilong Feng. Ziyun Huang, Jinhui Xu, and Jianxin Wang. Linear time algorithms for k-means with multi-swap local search. In Proceedings of the 37th International Conference on Neural Information Processing Systems, pages 45651 45651-45680, 45680, 2023.
[22] William B Johnson, Joram Lindenstrauss, et al. Extensions of lipschitz mappings into hilbert space. Contemporary Mathematics, 26(189-206):1, 1984.
[23] Tapas Kanungo. David M Mount, Nathan S Netanyahu, Christine D Piatko, Ruth Silverman, and Angela Wu. A local search approximation algorithm for k-means clustering, In Proceedings of the 18th Annual Symposium on Computational Geometry, pages 10-18, 2002.
[24] Max Dupré la Tour and David Saulpic. Almost-linear time approximation algorithm to euclidean k-median and k-means. arXiv preprint arXiv:2407.11217, 2024.
[25] Silvio Lattanzi and Christian Sohler. A better k-means++ algorithm via local search. In Proceedings of the 36th International Conference on Machine Learning, pages 3662-3671, 2019.
[26] HW Lewis. Medians and means in probabilistic risk assessments. Engineering, 91(2):220-222, 1985.
Nuclear Science and
[27] Stuart Lloyd. Least squares quantization in PCM. IEEE Transactions on Information Theory, 28(2):129 137 1982.
(28) Konstantin Makarychev, Yury Makarychev, and Ilya Razenshteyn. Performance of johnson- lindenstrauss transform for k-means and k-medians clustering. Proceedings of the 51st Annual ACM SIGACT Symposium on Theory of Computing, pages 1027-1038, 2019.
[29] Ramgopal R Mettu and C Greg Plaxton. Optimal time bounds for approximate clustering. Machine Learning, 56(1):35 (1:3:)35366)) 66, 2004.
[30] Jiayang Ren, Kaixun Hua, and Yankai Cao. Global optimal k-medoids clustering of one million samples. In Proceedings of the 36th Annual Conference on Neural Information Processing Systems, pages 982-944, 2022.
[31] David Sculley. Web-scale k-means clustering. In Proceedings of the 19th International Conference on World Wide Web, pages 1177-1178, 2010.
[32] Mikkel Thorup. Quick k-median, k-center, and facility location for sparse graphs. SIAM Journal on Computing, (4(2))405-433) 2005.
[33] Mo Tiwari, Ryan Kang, Donghyun Lee, Sebastian Thrun, Chris Piech, Ilan Shomorony, and Martin Jinye Zhang, Banditpam++: Faster k-medoids clustering. In Proceedings of the 37th International Conference on Neural Information Processing Systems, pages 73371-73382, 2023.
[34] Mo Tiwari, Martin Jinye Zhang, James Mayclin, Sebastian Thrun, Chris Piech, and Ilan Shomorony. Banditpam: Almost linear time k-medoids clustering via multi-armed bandits. In Proceedings of the 34th International Conference on Neural Information Processing Systems, pages 10211-10222 2020.
12

---

Technical Appendices and Supplementary Material
A Tree Construction and Rejection Sampling
A.1 Tree Embedding
The formal description for the tree embedding method presented Algorithm f Given a k-means clustering instance (P, k; d), the algorithm first scales the data points such that the minimum pairwise distance data points in P is (step of Algorithm 4). Then, it calculates an upper bound dmax for the maximum pairwise distance between data points in P using greedy strategy (step 3 of Algorithm 4). which takes ((d)) time. Then, a random shift S S dmax (step Algorithm 4 applied to each coordinate of the dataset P. The root node of the tree represents axis-aligned hypercube with side length 2dmmax that contains the entire dataset (step Algorithm/4 Then, the root node is iteratively decomposed using tree decomposition algorithm (step of Algorithm 4 to form a tree structure. First, the root node r decomposed into 2d cris-aligned subcubes (step 3 of Algorithm 5). where the obtained subcubes assigned to the subcube that fully contains its coordinates. For every non-empty subcube, we create node and add it child of the root (step 5 of Algorithm 5) This process is recursively repeated for each node until the side length of the subcube is smaller than Va (step I of Algorithm 5). The resulting tree structure T has leaves the height, with maximum tree height of H O(lggdAA) (A is the aspect ratio of the given clustering instance), and at most n. nodes in each layer.
associated with side length d'nax and each data point
Algorithm 4 Prre-Construction(P, k; d)
Input: A k-means instance (P, k; d).
Output:
embedded tree structure 1: Rescale the data points in such that the minimum pairwise distance between points P is 1.
2: Initialize empty tree structure TT
3: Randomly pick a data point p from P, and set d'mux = m 6(g. p).
4: Add a random shift
d'max to each coordinate of each data point in dataset P. 5: Create a hypercube r centered at with side length 2d and set r as the root node for T.
6: Call the Tree- Tree-Decomposition" , 2d max) algorithm to recursively decompose the root node.
7: Assign level le(v) to each node v E T, where leaf nodes have level I with levels increasing progressively toward the root.
8: return T.
Algorithm 5 Tree- Prre-Decompositioon77 (T, BB)
Input: tree structure a tree node B with side length (g. Output: A decomposed tree structure T.
1: lg < 1/Vd then
2: STOP.
3: Decompose B into 2d axis-aligned subcubes with side length 4g/2
4: Assign data points to the subeubes that fully contain them.
5: Set the non-empty subcubes as children of B.
6: for each non-empty child node E' of 1 do
7: Execute Tree-Decomposition(? (1/2))...
A.2 Binary Tree Construction
Our objective is to maintain a weight distribution for each node in the tree structure T to simulate D²-Sampling process by traversing from the root node to the leaves. However, since non-leaf nodes can have up to (n) children for the tree structure T. the sampling process on the tree structure may involve (n)) branches, leading to running time of 22 log(dA) (d)dd))). To address this issue, we propose an algorithm to transform the tree structure into binary tree T'' while preserving the pairwise distance approximation on the tree metric.
13

---

The tree conversion algorithm is presented in Algorithm 6 For each node in the tree structure T with more than 2 children, the tree conversion algorithm creates auxiliary nodes to connect the children and their parent to construct the binary tree structure (steps 3-5 of Algorithm 6 To maintain the expected approximation of pairwise distances, auxiliary nodes are assigned the same level as their parent nodes (step 7 of Algorithm 6). The following lemma shows that this process only adds an (log factor to the tree height while preserving the expected pairwise distances.
Lemma 2. Let T' be the modified tree structure of the output for Algorithm 6 Then, each vertex in T has at most 2 children, and the height of the tree T can be bounded by Olog 2 (dA)). Additionally, it holds that 8(p q) S 87 (p, q) and EleT (p. q) S &(dd) 8(p,)) q).
Proof. Consider an arbitrary non-leaf node e E T. Let Ue, Vey ... Uc. be the children of v. In Algorithm 6 it adds s auxiliary vertices '/2, Vey...., 1'c,- According to the connection rule stated in steps 4-5 of Algorithm each auxiliary vertex V/2 for some [2, s] can have at most 2 children e'y,, and U Therefore. each vertex in T' has at most two children. Since each node has at most n children T (where is the data size), adding n auxiliary nodes is enough to connect the children with the parent and auxiliary nodes. This adds a tree depth of at most (logn) for each node. Therefore, the overall height of the tree T' can be bounded by Ollog"(1)) (n)).
has only one child V'2 By renaming Vegs..., Ve, as Ve,11 Ve, 12 Vez.'
For the pairwise distances approximation, we consider any two point p, E P. Let Up.g be the node representing the lowest common ancestor p and g T. Denote v/,, as the lowest common ancestor of p and q T'. Since auxiliary nodes maintain the same tree level as their parent nodes, it follows directly that le(U2.) = le(o).g) Consequently, the properties of the approximate pairwise distance stated in Lemma I also holds for the transformed tree structure T',
Algorithm 6 Tree- Tree-Conversion(7)
Input: tree structure T.
Output: A new tree structure T'
1: Add all nodes of T into T', and keep their original levels.
2: for each v E T do
3: Consider each non-leaf node T with more than 2 children, where its children denoted as Go Vegs .... Ue,, Create & auxiliary nodes Ve Vep Ve.
4: For each E [s]. set the parent of Ve, in T as "(444//11""
5: For each i E [2, s). set the parent of vic, in T as v/ "/aj""""
6: Set the parent node of v',, as v.
7: Assign each auxiliary node level that is the same as its parent node.
8: return" T
A.3 Maintaining a Dynamic Data Structure D
Once the tree structure constructed, clustering centers are opened by marking the corresponding leaf nodes as active. However, during the local search process, center swaps can significantly influence the tree structure constructed. Hence, to overcome this issue, we propose TREE-INIT, TREE-OPEN and TREE-CI.OSE operations to maintain a dynamic tree structure during the local search process while maintaining the pairwise distance approximation. In general, the TREE-INIT operation sets up the necessary components for the dynamic structure D by traversing from the leaf nodes to the root. TREE-OPEN marks a leaf node as an active center and updates 2 along the path from the leaf to the root, while TREE-CI -LEE-C deactivates a marked leaf node and also updates the structure along the path from the leaf to the root. In the following, for better consistency, we also list the properties for the data structure as presented in the main context. Let C be the set of centers dynamically adjusted by the local search process. The data structure D on T'' should satisfy the following properties.
Each non-auxiliary node e in the tree T' is assigned level le(u), with leaf nodes assigned with level 1. Along any path from leaf to the root, the level increases by at each non-auxiliary node, while auxiliary nodes retain the same level as their parents.
14

---

Each node v is associated with three values ni,, ny and n(v). For non-leaf nodes, ni, and n denote the number of leaf nodes in the left and right subtrees of u, respectively, n(v) represents the number of active leaf nodes (or opened centers) that are descendants of w. Fach node u E T is
also associated with weight w(e) where w(v) = 8T (B., C).
Each non-leaf node v E T is associated with a probability list p(v),, which represents the proportion of the total weights contributed by each of its children. If has only one child, then p(e) = [1.
(B(,),)) If & has two children, p(v) = [L., rn), where L, (B. ,v (B. C) and l, + ru 1.
We first show that the above data structure can be initialized in time (nd log (dA)) by traversing the tree T' in bottom up manner, which call TREE-INIT operation.
Algorithm 7 FRREE-INI(P?
Input: A tree structure T'.
Output: A tree structure T' initialized with data structure D.
1: For each leaf node e E T', initialize ni, = 1, n O, n(()) = 0. w(v) d'axe where d' d'max is
the estimation for the maximum pairwise distance during Tree-Contto (Algorithm 4
2: Enumerate in bottom-up manner 'T' such that the children of a node always visited earlier
than the parent node.
3: for V E T' do
4: if v is non-leaf node then
5:
e has only one child o' then
6:
ni, = n²y + n.g, n. 0, n(v) 0, w(v) w(t'). p(v) = [1.0].
7:
else
8:
Let 0i and v/2 be the left and right child node for v.
9: ny = n + mii ny = m/g mig n(v) = 0, w(v)= = w(v() + w(r)).
10:
n n5 p(v) tng m
Specifically, for each v E T', the level le(v) has been set when performing the Cree-Construction and Tree-Conversion algorithms. By traversing from the leaf nodes to the root, nt nt, can also be initialized accordingly (step 6 and steps 8-9 of Algorithm 7). For each node v E TT. nn(v) is set as 0 since there are no active centers (centers that are opened) during the data structure initialization. For each node u with single child, the probability list p(v) is initialized as p(v) [1, 0] (i.e, v has only a left child node). On the other hand, for each node e with two children, the list p(v) is initialized as
n!. n p(u) ni +n o! + J. For each leaf node u E T', its weights w(v) is set as w(v) d' d'axx The weights of the internal nodes are updated by taking a summation of the weights of its children.
Next, we introduce the TREE-OPEN and TREE-OPEN operation marks the leaf node Up associated with as active and then traverses along the path from up to the root node for updating the probability list and weights for the nodes. The formal description for the TREE-OPEN process is described in Algorithm g
CLOSE operations. Given data point p P. the
The TREE-OPEN algorithm updates the tree structure T when a clustering center is opened. It begins by marking the leaf node ve corresponding to the opened center as active and setting its weight to zero (step I of Algorithm 8) The algorithm then traverses up the tree to the root, updating weights at each node based on whether its child nodes contain opened centers. Depending on whether the child node is left or right child node, subtree weights are adjusted accordingly (steps 4-15 of Algorithm 8). The total weight of cach node is updated as the weight sum of its left and right subtree,
and the probability list is updated accordingly (step 15 of Algorithm 8). This process dynamically adapts the tree structure while maintaining the data structure D, ensuring that the D²-Sampling can be approximated without requiring full tree reconstruction.
tree structure T', recall that we use I(v) and r(v) to denote its left child node and right child node, respectively. For the ease of analysis, if v has only one child node v', then we regard v' as the left child for v. Let B, denote the set of points corresponding to the leaf nodes that are descendants of node In the following, we use C' to denote the set of centers that have been opened before executing a single local search step. Then, data point p E P is ready to be inserted to C' newly selected center. For each center g C U {p), denote Vg as the leaf
Given a node
15

---

Algorithm 8 TREE-OPEN T POOOPPEENTT
Input: A tree structure T', a center p to be opened.
Output: A tree structure T with marked leaves, updated weight, and probability distribution.
1: Let Up be the leaf node associated with p. mark Up as active and set w((pp) = 0, n((pp) = 1.
2: Traverse the tree towards the root forming a path uo, u1...., 2e, where Up Up and 2/2 = Ur'
3: for each i E [t] do
4: if 1 is the left child of o; then
5: (1;))
then
" = 2vd (2ic(v)) - 2) ni w/: = w(()-1).
6:
7:
else
8:
" = w( w = w(v--1).
9:
else
10:
if
nnu;) = 0 then
11:
wyy = 2v (2le(vs) -2) mby wh: = w(m
12:
else
13:
" v(((v:)) w,, = =w(vi-1)
14:
n(ns) (v:) + 1, w(v)) wh, + w
15:
Let p(ey) = [,,,, re.] be the probability list associated with node uy and set L, = w(lle:))
w([(es] [[ (r(ee)) and Tuy (((ce)))(((te)))
w((((.))
16: return T
node in T' that is associated with q. Define L(4) = (vo, Vi,-..,; wr} as the set of the tree nodes representing the path from the leaf node Vy to the root node ur, where Up = Uy and v/ = ey. Recall that ÖT (B,.. C*) = Ew'ts (T' C?) is the sum of the tree distances (defined on T') between the leaf nodes in B, to the centers in C. The following lemma shows that the probability list maintained for each tree node can well represent the D2- D2-Sampling distribution.
Lemma A.1. For each q E C' CCU {p} and a tree node v, E L(g), holds that Lu, (Bre CU(p))
by,
5y (Brco)) and w(vz)= = (B,, (p))
Vi;
(B,
Proof. The proof is based on an induction.
For each center q e C' U {p). denote Vi the active leaf node associated with q. Consider the base case for C = {s} where is the first center added to C. We will show that the properties stated in Lemma A.I holds for Ug. Note that through INIE-INIT, each non-leaf node E T' is assigned
P probability list p(v) = rel. where 4 -
Each leaf node is also assigned
n +
weight w(v) = d'max- For ease of analysis, without loss of generality, we can assume that 22 E C(s) is the left child of 21+1 E (ss.. The proof for the base case also uses an induction. Firstly, consider the node o1 E Cs)). Since thp = vs is marked as active, the weight for Up is set as ve(vo)
(step of Algorithm (). Then, for the node in we have Bi(,,) = } and wh set as wi w(vo) = (step 6 of Algorithm g which corresponds to the sum of the distances of data points in Bi(..) to the opened center s. Then, since n(v1) is set as 0 through TREE-ING -NIT, we have wo, 2v (22 22) -n (step 6 of Algorithm 8). Since all the paths from leaf nodes to the root node have the same length, it holds that 8 (a, 8) 4vd for each a E Br(x.) and 8y: (B)())))) 8) = 4vd- n. = wh By updating w(v)) as wh, + "h we have w((1) 8 (B.,
Now, assume for each i E [j 11 and p(v;)) = Moes rn.. it holds that 4. = (B
To. (1.t0)))))))))) and w(rz)= = (B,, 8). Then, we consider the node uy. Since n(v;)) is set as zero byd (B., .s)
(B,,)))))))
through TREE-INIT, INIT, we have wh, w(U--1) (B, s) and wt, 2vd (21e(1,) - 2) nt, (step 6 of Algorithm (. Since all the leaf nodes are at the same height and the auxiliary nodes do not increase the node level, it holds that 8T (a,s) 2vd. (2/e(s)) - 2) for each E Br(e,). Hence, we can get that 87. 5) = 2vd (21e(19) - 2) m,y Therefore, we have w'y, = 8T 8) and 8π 8T (Br(e,).)) (Br(e)) s). Then, it holds that w(ey) wh, +wi, = ôT (Bug, s) (step 14 of Algorithm & wh, =
16

---

which proves the case for vy. Thus, we can conclude that C' = { where s is the first center added to C,, holds for Lemma A.I.
Next, given a set C" of opened centers, we assume that for each q E C" and a tree node oy E L(1), it
holds that Lr = (0118110)00""" Te, = ST ",r"B,, (') and w(w;)) = 87. (B,,, C*). Let g be a center that (B,, FT.
will be added to C", we will show the case holds for C' U'U(g (g)
The TREE- OPEN operation traverses from the leaf node vg to the root node v, to update the weights and probability distributions. For ease of analysis, without loss of generality, we can also assume that
V E L(g) is the left child of vi+1 E [g)- The proof for C" U {g} is also based on an induction. We
first consider the node u1. Since to f/g is marked as active, the weight for to is set w(vv)
(step 1 of Algorithm g). We have B(()n) = { and w', is set as 0 (step
step 8 of Algorithm 8
which corresponds to (B((,)) C' U {g]). Then, for analyzing the value of w there two subcases that may happen: (1) n((1) > O; (2) n(u1) = O. If subcase (1) happens, en must have been visited before g is added to C" by the TREE- OPEN operation, and e1 belongs to L(q') for some o' E C. According to the induction assumption, we have wo, w(r(vi)) = Ô (Br(e)) C*) (step 8 of Algorithm 8) Note that since g does not belong to (v,)), adding to C" will not induce distance change (no lower common ancestors are added) for the nodes in By(r))- Therefore, wt, = (Br()) C*) = 8T.(Br C' {g)). If subcase (2) happens, according to step 6 of Algorithm 8 we also have w", =r,)(B(()))))))) Note that since n(v1) = 0, it holds trivially that wh, 87 (Br(n)): {g}) (Br())) CU {g}). For both cases, by updating w(vr) as wh, + w we have w(v) 87 (Bu, C" (g}).
Assume for i E [i - 1] and p(u) = (lr,, it holds that l,, (B BT.(B,,
To:
(Brt C(U(H)) and w(t)) (T.(B., C" U (9)). We then consider the node vy- We have (B.
(B. (()) according to the induction assumption (step 6 step 8 of Al- gorithm 8). Then, there are also two subcases that may happen: (1) n(u;) 0; (2) n(v;) = 0. If sub- case (1) happens. Uj must have been visited by the TREE-OPEN operation and by belongs to L(y') for some q' C². According the induction assumption, we have wi, = w(((v))) = (B, C"). Note that since 9 does not belong to Br(v,), adding g to C' will not induce distance change for the leaf nodes in Br(ng) and hence wi, 8y. (Br(sy) C*) = (Br())))))) C' U {g}). If subcase (2) happens, we also have wi, = (T((Br())))))) Note that since n(v;) = O, it holds trivially
wh, =w( =
that wt, = 8T (Br (B))))) 8T C' {g]). By updating w(v,) as w/, + wh,- have w(t;)) C" {g)).
Finally, consider the tree nodes in L(q') for some q' E C" that do not belong to C(g). For each such node v, a key observation is that g does not belong to B,,. Hence, for each node o' E B,, the tree distance between v' to its closest opened center on the tree remains the same, where we have (B(o), and r, eBr (Bt,,, and w(v) (B,,, C (g))).
(B,
Putting all these together and using the induction, Lemma A.1 can be proved.
Lemma A.2. Assume that the current set of active centers (or opened centers) the tree T' is C. Given a data point p E P. TREE-OPER operation runs in time Oglog" (dA)) update the data structure D such that each non-leaf node v E T holds a probability list p(u) = [o,Tul. where (.(Buu,,CUPp) and r. (B,,,, (B.
1. (B,,),)
Proof. According to Lemma A. the properties stated in Lemma A.2 hold for each node v E L(q) where q C {p]. Then, we only need to consider the tree nodes that do not belong L(g) for any E (p). For each such node e, observe that is never visited during the TREE-OPEN operations, and hence B, does not contain any leaf node that are associated with the centers in C. For each such node v, since no leaf nodes that are descendants of v have been marked as active nodes, the tree distances between any leaf node v' E B, to the opened centers are the same. According to TREE- (,((B(( and (B.,
INIT, we have L,
R.. the and r, mo
ni.
-. Then, it holds that 7, ni. +
8 b (B,, (U(p)
(u
Tv
17

---

As for the rumning time, the TREE-OPEN process traverses from the leaf node associated with to the root node. Since each node along the path has at most 2 children and the tree structure T has height (dA)) the running time can be bounded by O (dA)).
Then, we can present the TREE- CLEEE operation. Based on set C of centers that have been opened, the TREE-CLOSE operation aims to delete center p E C. Specifically, the TREE-CLOST operation deactivates the leaf node Up associated with p and then traverses along the path from up to the root node to update the probability list and weights for the nodes. The formal description for the TREE- CLOSE process described Algorithm 9 The TREE-CLOSE algorithm updates the tree structure T' when clustering center is removed, ensuring the data structure D remains consistent. It begins by deactivating the leaf node Up associated with the p (step 2 Algorithm 9 Then, identifies the first ancestor yf that still has least two remaining active leaf node in its descendants (step 3 of Algorithm 9.. If o' does not exist, it assigns a maximum weight d' d'max for Up (step 5 of Algorithm 9). Otherwise, it updates the weight based on the tree level of 1 (step of Algorithm 99. The algorithm then traverses from Up to the root, adjusting weights and probability distributions at each node to reflect the center closure. During the updates, the weights structure (steps 9-25 of Algorithm 99. The probability list is updated accordingly to ensure that the tree still well approximate the D22-Sampling process. By dynamically modifying the weight distributions and probability values, TREE-CLOSE ensures that center closures are properly handled without requiring a full tree reconstruction. The following lemma shows that the probability list maintained for each tree node can still well represent the D²-Sampling distribution based on tree distances between data points. Since the proof for Lemma A.3 is similar to those for Lemma A.2 and Lemma A.1, we omit the proof here.
recalculated to maintain the hierarchical
Algorithm 9 TREE-CLOSE(T') .p)
Input: tree structure T'
Output: A tree structure 'I' with updated leaves, weights and probability distributions.
1: Let up, U1,..., be path along the leaf node Up associated with p the root node Ur, where Vo = Up and 1/2 = Uy.
2: Deactivate Up. and set n(u,)) n(up) 1.
3: Traverse along the nodes Up to find the first node v' such that n(v') > 1.
4: if o' does not exist then
5: Set w(e))) = d'axx.
6: else
7: w(v)) = =2vd- (2(e(v")) -2). 2).
8: for each i [t] do
9:
if 0;-1 is the left child of then if (n(u)) 1= 0 then if u' does not exist then
11:
w = d'nax nh,,, w d'nax n
13:
else
w 2/d (2le(0') -2) nhy w/, 2/d (2te(v)) 2) n/y..
15:
else
16:
1. w(((v))), w = w(vi-1).
17:
else
18:
n(.) - 1- then
19:
if u does not exist then
20:
w = -'nax ny w = d'nax n/..-
21:
else
w., = 2vd (2le(o') - nhy, w/y = 2vd (2le(1') - -2) nl...
22:
23:
else
24:
w/ w(((v,)), w., = w(w;-1).
25:
re(1)) (vx) 1, w(v)) wh, + :
26:
Let p(v:) = Fey Foy] be the probability list associated with node v. and set L, = w(r(e:))
=(le:))+=u((e:)) and ru. v((c)))+((((())))
w(l(e:))
27: return T'.
18

---

Lemma A.3. Assume that the current set of centers opened is C, where C also corresponds to C, PREE-CLOSE operation runs time Olog" (dA)) to update the tree structure such that each non-leaf node T holds a probability list p(u) = Tol, where 1, (, (Mpp) and r. 6 (V) (B((B.
the set of active leaf nodes T'. Given a center
Putting all these together, Lemma can be proved.
A.4 The Adaptive Sampling Strategy
Based on the TREE-OPEN and TREE-CLOR operations, during the whole local search process, always maintain tree T' such that each non-leaf node E T' holds a probability list (Brcs)))
p(u) = [4,roo, rol. where 4 = (Br(Bu)) and r. = 5 (B However, since the tree structure (B,(B_,
only guarantees the approximation loss in expectation (see Lemma 1 and Lemma for details), directly performing a sampling on the tree structure may lead to unbounded success probability. Partly inspired by previous work [15], we propose an adaptive sampling algorithm to overcome this issue, which coverts the expected distance-basee sampling distribution into the true distance- bassd sampling distribution. The formal adaptive sampling algorithm is presented in Algorithm T0
Algorithm 10 Adapiive-Samplnn((7 CC)
Input: tree structure T', set C' clustering centers. Output: A sampled point æ.
1: S = 0.
2: while S = do
3: Let v be the root node of the tree T',
4: while v is not leaf node do
5: Let p(v)
[4, ry] be the probability list for 6: Update v to be its left child with probability 4 and to be its right child with probability Te-
p(v) is detailed in Section 3.1
7: Let be the point associated with the leaf node v.
8: Add æ to S with probability min { (,. where 8 C) is the tree distance between
6(x.C)
point to C' in the tree structure T" return S.
D. Tree distance is detailed in Section 3.1
9:
The Adaptive-Sampling algorithm samples data point by leveraging the tree structure T' to align with the D²-Sampling distribution. The process begins by initializing an empty set S and repeatedly drawing samples until valid point is accepted (steps 4-8 of Algorithm 10p. The sampling process starts from the root node and traverses the tree downward, selecting the left right child based on the probability list associated with each node (steps 5-6 of Algorithm TT). Once leaf node is reached, it identifies the corresponding data point x and applies rejection sampling step (step 8 of Algorithm 10 where x is accepted with probability min 1, 6(xrCC }. The process repeats until a
(x,C),
point is successfully accepted, where the sampled point then returned by the algorithm.
The following lemmas show that the rejection sampling process can well approximate the D2- Sampling distribution with bounded expected rumning time.
Lemma A.4. The probability of inserting a point æ to the set in Adaptive-Sampling algorithm is 8(a.C) Σpen 8(p.05
8y-(2.C) 1. Based on Lemma A.2 and Lemma A.3, during the local search process, the tree T always maintain a probability list p(u) [4,, for each non-leaf node T' where 4 (8,
Proof. According to Lemma 2 we have (x, C) 2 8(z, C). and thus 5(a,C)
645
Hence, the probability of sampling data point & P becomes (B.(B..C) and r, 5 (B..
(x, 8(x,C)
Σperdr:
(T. C)
(x,
19

---

which is independent of the term (2, Since the term Σpep (p, C) is the same for each E P and does not depend on 2, the probability of adding any point to S is 5(a.C) according Eper 6(2,C)
to the rejection sampling rule.
Lemma A.5. The expected number of iterations for Adaptive-Sampling to return a sampled point P can be bounded by Ö(d).
Proof. Consider a single step for sampling data point E P using Adaptive-Samplng According 8(a.C)
Lemma A.4, the probability of adding a data point E to is Ever (p.C) Denote the random variable N the number of steps executed when at least one data point x E P is added to
S. Let y the probability that a data point x E added to within a single sampling step (i.e.,
steps 3-8 of Algorithm T0 Then, we have
Exep 8(2, C) Σp 8T (p,
=
This indicates that data points are added to S in expectation within single sampling step. We can get that
+00
E[N] = Σ - (1-))=-1 10(-1))-0 i=1
q
= 1/4.
Ege (p.C))
By combining with Lemma we have E|N]
for a fixed tree T'. By taking the Exce Exen ((x,C) < Ö(d), which
Zeer 6(x,C)
expectation over the randomness of tree-embedding, we have E proves the lemma.
Lemma A.6. The (2,(C))
takes expected time Ö(k) to sample a data point E wiith
probability at least
V1S-(P,C)
Proof. According to Lemma A.4, get that the probability of sampling data point x E P is 8(x.C)
Epen Slp.Cy. By using Caushy Inequality, we have
2
Σ PEP
Σ 12
Σ PEP
Σ5(..
2
PEP
which implies that
((P.C) Vn. (8)(P.C)
Then, can get that each data point E P can be sampled with probability at least 8(x,C) 8(P.C)
2
6(a.C) J
(x where the last step follows from (xxx,)) (P.C)
< 1.
(P(.
(P(P.C)
As for the rumning time, TREE-CLOSE and TREE-OPEN operations take time (log (dA)). Each sampling step traversing from the root node of the tree to the leaf node of the tree takes time Ologg" (dA)) (steps 3-7 of Algorithm 10). Querying the true distance between the sampled data point I to its closest center in C' takes O(kd) time (step 8 of Algorithm
10). According to Lemma A.5, the expected number of iterations for sampling a data point is (dd). Finally, with ((nd) preprocessing time before tree construction, the aspect ratio A and the dimension d can be reduced poly(n, d) and Oloog n.), respectively. The overall expected running time for sampling data point E P using Adaptive-Sumpling is Ö(k)
20

---

B Bandit-Based Center Swap
In this subsection, we present the bandit-based center swap method. For the existing local search methods, even with sampling-based swap pair construction, a clustering cost estimation process is required to pick the best swap pair among the candidate swap pairs. However, in the worst case, this process requires scanning the entire dataset for distance and data point reassignment calculations, which leads to a running time of (((d)) for each center swap step.
To further accelerate the best swap pair selection process, our proposed bandit-based method (Al- gorithm 2) can identify the swap pair inducing the best clustering cost reduction with sub-linear runtime assuming sub-Gaussian prior on swap pairs. The high-level idea behind is to reformulate the deterministic clustering cost estimation process as a statistical bandit problem. By regarding each swap pair as an arm in the bandit framework, selecting the swap pair that yields the best clustering cost reduction can be modeled as a best-arm identification problem. To give the detailed analysis, we first introduce the following folklores for statistics.
Lemma B.1. a < X S holds with probability 1 for a random variable X, then X (5 Suu--aussiart where ( is the variance proxy.
Lemma B.2. (Hoeffding's Inequality) Let Xn be i.i.d. sub-Gaussian random variables with variance proxy o². Then, for any 0, it holds that P. (1 Ei= X; F[X] 2 r) Se 202
We begin the formal analysis of the proposed bandit-based center swap method by first presenting a reduction of the best swap pair selection problem to a bandit problem. Given group S; (i = 1, 2) containing several swap pairs as the input for Algorithm 2, let Swap = S; and See = P be the sets of the reference swap pairs and reference points, respectively. In the modeled bandit problem, each swap pair s E S is regarded as an arm and the goal is to identify the swap pair (u, v) E Sswap with the maximum clustering cost reduction, i.e., (u, v) = arg min (') (P, CU (u)A)(())
}). Given a data point randomly sampled from the dataset P. for each swap pair s = (u', (u','') E Swapp the clustering cost change induced by (u', e²) on b Seer is defined as (b) (b, ( (v())
(2,b, C²). Given set H. of data points, let G2 (HL)= = dEH (b) Define μ.. G.(P) Sperl The best swap pair also has the smallest value of G2 (Pe/)Seet.
The bandit-based center swap algorithm is outlined in Algorithm 2 (denoted as the BanditCS algo- rithm). BanditCS algorithm iteratively samples data points from the referenced points to estimate the average changes in clustering cost for the swap pairs, where confidence interval is maintained for each candidate swap pair s E S In each iteration, the algorithm updates both the mean estimate م. of the average cost reduction and the interval length F, (step 7 of Algorithm 2) for each 8 E Sswap based on the sampled points. From these updates, confidence interval [is - Fs, fis + FF.] is computed (step of Algorithm 2). The confidence interval captures the uncertainty around the estimated average clustering cost reduction based on the current samples. As more samples are collected, the intervals are narrowed to reflect the true average clustering reduction for the swap pairs, The algorithm then climinates swap pairs with less promising estimates (step 9 of Algorithm 2). This process continues until only one swap pair remains or predetined number of sample sizes reached. We show that under the Sab-Gausiin prior on swap pairs with significant clustering cost reduction, the bandit- based algorithm (Algorithm 2) can always return the best swap pair with high probability.
Assumption 1. Let C be the set of centers before executing each step 10 of Algorithmm [ For each swap pair s Si U Uz. it is assumed that G, (b) is o² ---ub--auussann for a randomly sampled point b E P, with a known bound o² = O((5² (P, C) ())}) and o² 2 d²
Lemma 5. Let S be a set of swap pairs as input for Algorithm 2 which contains at least one swap pair with a (1 - 1/100k) fraction of cost reduction. For n = On +)), Algorithm 2 can remove non-promising swap pairs in S with high probability using ( samples.
Proof. We first prove that the best swap pair will never be climinated during the whole bandit process with high probability. The proof begins by showing that with probability at least 1/n?, all confidence intervals maintained by Algorithm 2 can well approximate the real mean of the clustering cost change for each swap pair 8 E Syway. Let B denote the set of the sampled points obtained after
'A random variable X is o² --u--uu--uuue Gau-Gaus if P.|x > 1] 2e 62
21

---

executing step 4 of Algorithm 2. Then, we consider fixed swap pair E Sswap. Let Pa be the average of the [B] i..dd samples in B drawn from Stt where Ser P. Since each (b) has a bounded value, random variable X, G, (B) can be regarded as o--sub- suu--ausssinn variable for some o3 > 0. By using the Hoeffding's Inequality (Lemma B.2), we can get that
PA(lp Pol F.)
2m3
According to the definition of F, (step 8 of Algorithm 2), we have F, = a.
ll((1/m)
where is
Rett+
the input parameter for Algorithm to control the success probability, nhar the number of samples that have been taken, is provided as input to Algorithm 2 as prior knowledge (serving as upper bound on each o2). Then, we have
18 2l((//))
2e
207
2e
207 (1)
2e In(1/n),
where the last step follows from or + 1 1 (B] and 0² > o according to Assumption 1.
Let 17 2n We have Pr(lu, A. > F.) S 1/n². Then, Ns [A, F.. i + Fs] holds for each 8 E Swap with probability at least 1 + by taking union bound over the success probability of all the swap pairs in Sywap- Since there are most [P] sampling steps, can get that with probability at least 1 ,. the intervals maintained can always reflect the true mean of the clustering cost reduction, ie., Ms E [Ps- Fs, + F.. holds for each s E Sswap during the whole bandit process.
Let s* = arg min (P) be the optimal target swap pair. Since A.- F.- S Mar, we have A..- FF,- S ps S. fig' + Fa' holds for any s' E SSwapp such that s' * s" This indicates that the optimal swap pair with maximum clustering cost reduction will never be removed.
According to the grouping strategy, there are two cases that may happen (see Lemma B.3 for details): (1) either S₁ S2 contains swap pair (u, v) such that 5² (P, C* fo) U fu}) I (1 - 1 )82 (P, C); (1 TooE )62 (P.C).
(2) (2) and S. do not contain any swap pair (u, v) such that 82 CP (w))
If case (2) happens, since each swap group contains a virtual swap pair with zero clustering cost change, the optimal swap pair will not be removed. Hence, if the of swap pair returned by Algorithm 2 has size 1, the swap pair will not induce clustering cost increase. Then, we consider that
case (1) happens. For swap pair Swaap= let
Ms 142 For nhet
200g In n. have
G
2(F, + F.)) 40 / 2n((n)) 1)
S So
1b 14
Then, we can get that
No F. 2 us - 2FF
u4. + Co- 2F,
M + 2F,
As F..
This implies that after randomly and independently sampling 200 In points from Siet the swap pair Sywap can be removed from Sawap because its confidence interval is entirely outside the confidence interval of s"
Furthermore, in case (1), the given swap pair set Sywap contains a swap pair s* = (u, v) satisfying 52(P. C (u](fe]) s (1 100k )82 C) and virtual swap pair = (u, u) with zero clustering (PP.C) Since the virtual swap pair
cost change. Then, it holds that s* = (u, v) and Ma* S Sutt 100k
can only induce zero clustering cost change, we have C.- 2 [8a] Then, observe that o² S (( C)/\S,al)) according Assumption 1, can get that
52(P.C)
200g² Inn < Ollog)) o ((81P.1)) (100k)Swrl2
Nhet
2
S
(6Z(P, = O(k log n),
22

---

which implies that by randomly and independently sampling Ö(k2) points from reference point set Stets the best swap pair in S can be returned by Algorithm 2.
During the swap pair construction process, the BanditLS algorithm (Algorithm 1) constructs swap pairs by finding the center 8z E C closest to and a randomly selected center o' E C \ {sy) (steps 8-9 of Algorithm I), where A = {xx 8x), (z, q'}} serves as the candidate swap pairs. The following lemma shows that if is sampled with probability 8 (x, ()/62(P, C), with constant probability, there exists at least one swap pair that can induce certain fraction of clustering cost reduction.
Lemma B.3. (Fan et al. (99) Let x be data point sampled with probability (xx (), P. If (P, > PT,0OPT, then either with constant probability, there exists at least one swap pair (u, v) E A such that (P, C [u]{d)))= S (1 10ok
C) from
1 )8² (P, C). or with probability ((2/k)),
there exists least one swap pair (u, u) E A such that 8 CU (u)A(u) S (1 T0ox )6 (P, C). where X denotes the number of lonely centers (1 S X S k 1).
According to Lemma 5, for the swap groups S₁ and S2 provided input to Algorithm 2, the algorithm can return a set of for Si and a set of for Sg. Algorithm guarantees that if |o = 1. then of contains the best swap pair of S; for each E [1, By taking union bound success probability over the entire local search process, Algorithm 2 can always return the best swap pair with constant probability for each input swap pair group S, According to Lemma B.3, with certain probability, there exists one swap pair that can induce at least a (1 - 1/100k) fraction of reduction on clustering cost in each local search step. Then, with high probability, Algorithm can return set of of swap pairs with |oo1 = 1. Since the final swap pair is randomly picked from of with |oo| = 1, we can get that with certain probability, the clustering cost can be reduced for each local search step.
Lemma B.4. After (ynk (n)) rounds of local search steps, Algorithm I can return a constant
approximate solution in expectation.
Proof. According to Lemma 5 and Lemma B.3, if sampling with distribution (1 (P,C) (where z E P), then there are two cases that may happen: (1) with constant probability. the clustering cost can be reduced by factor of at least (1 100k for a single local search step; (2) with probability
((A/k), the clustering cost can be reduced by a factor of (1 robx). Based on Lemma A.6, each data
point can be sampled with probability at least (a,
(n.8((P.C) (where x E P) using
strategy. Hence, either the clustering cost can be reduced by a factor of at least (1 100K with probability 22(
(nn or the clustering cost can be reduced by a factor of at least (1 - Toox with probability Q k.7n Let R 2ark log(n where is a sufficiently large constant. We begin by examining case (1) of the local search process. Following the formulation in 225, we define a random process X that starts from an initial clustering cost SP(P, C*), where C' is an )( >aapprxximation. During R/2 iterations, this process reduces S²(P, C²) by multiplicative factor of (1 100k
1 ) with
probability 1 B while adding an additive term of 20000PT after R/2 iterations for some constant, 9. From this, we can derive that
R/2
E[X] 20000PT + (P,.C)) Σ
R/2
(1-))) (1-1) R/2-0 (1
1
100k
=((PC') (1
1J
R/2
20
100k
8PP.C") < n + 20000 PT.
Next, we consider case (2) and define a new random process X'. Over R/2 iterations of sampling and swaps, the process decreases 82 (P, C') by a multiplicative factor of (1 1 with probability
n Fynkk and increases the final cost (P,, C*) by an additive term of 20000 PT for some constant 3. In this case, we can similarly derive that
23

---

1/2
R/2
- n/) (1 100X
B[X*] 2000OPT (P,CC))
Σ
(1))
t=1
= - (1
r/'
R/2
2000OPT
100A
8P(P,C" 20000PP. nd
For R = 2nkyn (n4) iterations, there are at least /k//n log(nA) iterations such that either case (1) case (2) happens. It is easy to see that the random process X X' can dominate the local search process and hence [P, C)] < E[X'] and (P, C)] < E hold, where C' is the set of the final clustering centers returned by the local search algorithm. Then, can get that
(P,0)) & (| (18( S E(X')
S (E((PCC)))))
P((( (P((.C") + 20000PP)
nd
(P,
+ 2000OPT.
nA
Since the C' is an O(nA)-ana apprrxxmmae initial solution, we have E16 0(1)PP by setting a as a large enough constant.
Putting all these together, Theorem can be proved.
Theorem 1. With sub- Sab-Gaassian prior on cost changes of swap pairs, there exists an algorithm for k-means that can achieve constant approximation within expected running time Ó + /nkk*)
Proof. The approximation guarantees clustering quality follows from Lemma B.4. Then, analyze the running time. The tree initialization (steps 1-3 of Algorithm 1) takes time O(nd log (nt)), )). Note that A and the dimension can be compressed A poly(n, d) and Olog respectively, using the methods proposed in [18 28 Hence, the running time for initialization becomes Ö(n)). During the local search steps, there are most ((vrk) iterations. The Adaptive-Sampling takes time Ō(k) for each local search step according Lemma A.6. For the bandit-based center swap process, finding the closest center in C' for the sampled data point æ takes time Ō(k) in step 6 of Algorithm For each swap pair constructed, the bandit-based clustering cost estimation takes Ō(k²) samples. The algorithm needs to calculate the true distances between the samples to their nearest centers, which takes ((k3) time. Hence, the bandit process takes overall running time of Ō(k*) (steps 9-13 of Algorithm 1). Finally, updating the tree structure through TREE-OPEN and TREE- CLESECLOE operations take time ((1) for each local search step. Since there are most Ö(v local search steps, the total running time Ö(nd vnk²1).
C A More Practical Algorithm with Metropolis Hastings
The Metropolis- Haaropolis-Hastingss strategy is widely used in statistics for sampling from complex probability distributions, where several initialization schemes with sub-linear (or linear) running time were proposed 4. In the following, we show that under mild assumptions, the Metropolis-Hastings strategy can be used to approximate the D2-Sampling distribution within Ó(n rounds of sample transitions during the local search process. The following lemma is folklore for Metropolis-Hastingss.
Lemma C.I. (Chib and Greenberg (T) Let æ denote the target distribution defined on P that wish to approximate using Metropolis Hastinggs process. Given parameter 3 > 0, if the
24

---

proposal distribution q satisfies that g(x) 2 [T(x) for all E P, then after T = O(1/3) steps, the distribution of the Metropolis Hasting chain dominates 7/2.
The following lemma shows that under a tailed behavior of distributions on data points, the ratio between the maximum distances to the optimal clustering center and the average clustering cost can be bounded by a function that is sublinear in the data size
Lemma C.2. (Bachem al. (3) Let F be probability distribution over IR with finite variance that has at most exponential tails, i.e. 3c, t such that P. P, p) > aj I ce where e, t are constants. Let X be set of points independently sampled from F. Then, with high probability, for sufficiently large n, a maxxex x 6 (x)
Ex'ex ex (P(r' = Ollgg* n).
Based on Lemma C.I, Lemma C.2 and Assumption 2, Lemma 6 can be proved.
Assumption Assume that the given dataset P average where each optimal cluster has size S2(k?). For each optimal cluster Ph: let cj be the optimal center Ph. We assume that each Př follows distribution F over Re with exponential tails, i.e., 3c, f such that P,(8 p) a] S ce holds for x E P% where C, f are constants and μ is the mean of the distribution.
sampling steps (steps 2-5 of Algorithm 3). if the given clustering instance satisfies the properties in Assumption 2, we can sample data point E P with probability least 0.56² CP/82
Lemma 6. After Ö(n/k k22) Metropolis-Hasse
Proof. For any data point p P and given set of clustering centers opened, let 8p arg min ò(p,c') be its closest center in C. We first show that if the given clustering instance satisfies the property of average sizes (i.e., each optimal cluster has size (2(k))) the distance between p and 8p can be bounded. Without loss of generality, we can assume that lies in some optimal cluster Pj where of is its optimal clustering center. Let OPT, be the clustering cost of Pj with respect to oj i.e., OPT, = (P) (G))) Denote Sez as the nearest data point in C to oj- Then, we
can get that 82 (c(),, S ming (c(c), (q) S EyEc; -))
20PT; +2
|c31
1C,1
O [OPT,+6(03
where the second step follows from the triangle inequality, and the last step follows from the assumption that |Cj S(k2) (Assumption 2) holds for each E [k]. Denote dmax= maxzery 52 cg). Then, have
82 (p, 8p) S (1, sx;))
< 28 ct) + +c)
$ 2d + (8' C)/K² + +OPT,/13).
We now prove that p(x) 2 (2)/ where p(x) is the probability of sampling using Ö(n/k3) steps of Metropolis Hastings (with a uniform proposal distribution 1/n) and T(x) = (2 (P,88 is the target D²-Sampling distribution. Let a for some large enough constant т. For each optimal cluster Pj. according to Assumption 2, we assume that Pj follows distribution F 10 Ollog n).. This with exponential tails. Then, according to Lemma C.2, we have a, (3))) implies that
2a2
((C),C)+OPI) 12
8 (p. =0
-log² n.
25

---

For each 2: E P. let ((x) be the index of the optimal cluster that x belongs to, i..,, E P d(x)- Then, it holds that
hz
(x,
Bm(x) S
Tn log n
(P(C.
h²
8 (x)
s Tm log" (P(PCC)
8' (CCaan
llog
s 12
12
n
log
n
S -, n
( holds for each Cj. since Cj C P and (P, C) 2 OPT. The proof can be concluded using Lemma C.I. A
where the last step follows from the fact that O (ICC
Putting all these together, Theorem 2 can be proved.
Theorem 2. With sub- Gaussian prior on cost changes of swap pairs and mild assumptions on optimal clusters, there is constant approximation algorithm for k-means with expected runtime Ond + k*).
Proof. We first prove the approximation guarantees clustering quality. According to Lemma 6, with O(n/ki Metropolis-Hasinee Hastings steps, we can sample data point x E P with probability at least 0.55" (x, (P. C') for each data point P. Then, similar the analysis for Lemma B.4, after Ō(k) local search steps, we can obtain a constant approximation in expectation.
As for the running time, each approximate D²-Sampling step based on Metropolis-Hastine strategy takes (n//k² transitions, where each transition needs to calculate the distance between the data point sampled to its nearest center in C. Hence, the overall rumning time for the approximate sampling process in single local search step Ó(nd/ k).. According Theorem 1, each bandit step takes time Ö(kad). d). Finally, since (k)) rounds of local search steps are required to return a constant approximate solution, the total running time is Ö(nd + k:4) using dimension reduction techniques.
D Extension to the k-median Objective
In this section, we show how to extend our proposed adaptive sampling method the k-median objective without relying on any data distribution assumptions. The formal algorithm is presented in Algorithm IT] where the high level idea behind replace the bandit process with a weighted sampling strategy (Algorithm) 12) to remove the assumptions on uub-Gaussian prior for swap pair.
Similar the idea for the k-means objective, the algorithm begins by embedding the given k- median instance into tree structure T', where the pairwise distances can be well approximated on the tree metric (steps 1-3 of Algorithm [T]. Then, during each local search step, an adaptive sampling process is used to approximate the D-Sampling (which refers to the sampling distribution of 6(p. C)/8(P, C') for given set of centers opened) process. Instead of modeling the best swap pair identification task as a bandit problem, the algorithm adapts a weighted sampling strategy to estimate the clustering cost induced by each swap pair (Algorithm [2). For each swap pair to be algorithm first performs the swaps on the tree structure T' based
estimated, the on TREE-OPEN and TREE- CLOSEE operations (steps 1-3 of Algorithm 12). Then, independently takes ((vī/c) samples using the adaptive sampling method (step 8 of Algorithm 12). The intuitive idea is to leverage the properties by adaptive sampling method such that each data point E P can be sampled with probability linearly related to the distance of to its nearest centers opened. This would give good approximation to the inverse value of the clustering cost for k-median objective using weighted summation of the samples drawn. With this technique, the time complexity for clustering cost estimation can be reduced from O(nk) to be sub-linear in the data size
26

---

Algorithm 11 L.S-KMedian(PPkkk.k
Input: A k-median clustering instance (P, ki d), set C' C Rdd of random clustering centers.
Output: A set C' C. Rdd of clustering centers.
1: T Trre-Coontruuctionf/ P, k, d). Tree-Construction is detailed in Algorithm 4 in Appendix A
2: T' (rree-oonversion(]]
Convernnvrsinn detailed in Algorithm 6 in Appendix
3: Call TREE-INIT(T") to initialize the dynamic data structure. > Algorithm in Appendix A
4: for c C do
5: Call the TREE-OPEN(T', ,c) algorithm to mark the leaf node in T' associated with active
and update the tree structure T'. TREE-OPEN detailed in Algorithm in Appendix A
6: for = 1 to (k)) do
7: z = Adaptive- Addative-Samplingg Sampling(T C) Adaptive-Sampling is in Algorithm 10 in Appendix A
8: Set 82 = arg mincec 6(2, c) and randomly sample a center o' E (x2(
9:
8x), 82), (z, q'), !')).
10:
for o; E o do
11:
F.ST (o:) = -eeighted-Sampling Sampling P.k. d. d. 04, o4, 1/400k). = [(r, 82). (8x, }. 02 = ((x,g)) (d'
12:
13:
Im argginggg migmingger E.S Imz arg gin EST
14:
for g1 do
if ES) 400k-1 EST( then
100k+T
16:
remove g from 91- for E g2 do if (g)
17:
18:
400k-1
FST (gmg) then
1006+1
19:
remove g from 92.
20:
Let = (g: gg = 1, E {1, 2],, and randomly choose a swap pair of from O.
21:
C= U where (u. v) o'.
22:
Call the TREE-CLOS CHREECLOSE(T', v) algorithm to deactivate the leaf node in T' associated with v
and update the tree structure T'. TREE-CLOSE is detailed in Algorithm 9 in Appendix A
23:
Call the TREE-OPEN(T'. u) algorithm to mark the leaf node in T" associated with u as active
and update the tree structure T
24:
return C.
In the following, give the detailed analysis for the proposed algorithm. The following lemma can be directly obtained via Lemma A.6, which is an extension from the k-means objective to the k-median objective.
Lemma D.1. The Adaptive-Samplinn takes expected time Ö(k) to sample a data point E witth 6(a.C) probability J(P.C).
According Lemma D.I, each data point is sampled with probability proportional to the distance to its nearest centers opened. In the following, we will show that by modeling the clustering cost estimation sum estimation problem, the clustering cost be well approximated within factor of + E) using Ö(Vn/e) samples.
For the clustering problem, given a swap pair (u, v), v), the clustering cost of 8(P, C. U (u}) can be decomposed into 8(P, C* (u) U (ut)= = 5(p, CV) U {u)), U which can further be modeled as estimation task. Our main objective is estimate the sum W Eper 8(p, C' (v) U (u))
using as small number of samples possible. The following lemma shows that if proportional sampling is allowed (i.e., the D S-mmplingg,.
can estimate the sum using ((Vrfe) samples.
Lemma D.2. (Beretta and Tëtek (8) Given a U where each item E U has a weight w(a), let
W Eaev w(aa) be the summation of the weights in U. Let S = {ai,..., am be set of elements sampled mn O( O(V/U//e) times from such that each a; E is sampled with probabilities proportional to their weights. For each E S, define c, to be the number of times s is sampled. Then,
W' = (():-) (ZOES
(
-1
gives an estimation for W that induces multiplicative error of (1 ee
w[s]
with probability at least 2/3.
27

---

Algorithm 12 Brighted-Samplin((P.. ,,, ,)
Input: A k-median clustering instance (P. k. d), a tree structure T', a swap pair s, and a parameter
0 < c < 1.
Output: The estimation of the clustering cost induced by the swap pair s.
1: Let s = (u, v) be the swap pair to be estimated.
2: Call the SHEE-CLOSE(T' v) algorithm to deactivate the leaf node associated with v and update
> 'TREECLOSE is detailed in Algorithm 9 in Appendix A 3: Call the TREE- PREE-OPEN(TT) u) algorithm to mark the leaf node associated with u as active and update the tree structure T',
the tree structure T'.
TREE-OPEN is detailed in Algorithm 8 in Appendix A
4: Initialize = 0, C² = Cv(v) (u).
5: for m = 1 to Olog n) do
6: Initialize an empty set .M= 0, and set mp = O.
7: for 1 to O( (//n/e) do
8:
8 Addapive-Sampllng
C"), mn = mnn +
9:
if s E .M then
10:
C - else
11:
12:
C = -1 (s}.
13:
W (mr.) EEEM MG..C))
w wu(
14:
Sort the values in W with increasing order and let W be the median value. return W.
15:
Lemma D.2 indicates that through D-Sampling (sampling with probabilities proportional to the weights), the weighted sum can closely estimated using O (v/n/c) samples. Since Lemma D.1 guarantees that the Adaptive-Sampllig process can sample each data point E P with probability 8(a, CVto U (u))//(P, CV(v) U (u)), the following lemma can be directly obtained.
Lemma D.3. Given a swap pair (u, v) and a set C of centers opened as the input for Algorithm 72 let C' = Cu) Uv be the set of centers obtained after performing the swap. For the m-th iteration of Algorithm 12 (steps 6-13), an estimation W' can be obtained in step 13 of Algorithm 12 such that W' E [[1 - () C²), (1 E))(P, C")] with probability at least 2/3.
However, since there are Ō(k) local search steps, need boost the success probability to guarantee a union bound of success probability over all the local search steps. According to the analysis in weighted sampling process (8) 1 1/W" is unbiased estimation for 1/W and the success probability is established by using Chebyshev inequality. Hence, based on a standard Median of Means trick (26) (step and step 14 of Algorithm (2). the success probability for finding an accurate estimation can be boosted to 1 - 1/n² by constructing Ollog n) groups and taking the median value. Hence, we have the following result.
Lemma D.4. Given a swap pair (u, v) and a set C of centers opened as the input for Algorithm 12 let C' = CV (v) be the set of centers obtained after performing the swap. An estimation W can be returned by Algorithm|] 12 such that W E [(1 - () C"), (P C²)) with high probability of at least 1/n².
Next, we will show that by carefully adjusting the parameter e and introducing a judging criteria, the clustering cost can be reduced by at least factor of (1 - 1 1000k) with certain probability for each local search step. The basic idea behind is to add a virtual swap pair (x) for some r E C' to serve baseline, where swap pairs that fail to induce clustering cost reductions can be filtered by the virtual swap pair. The following lemma shows that using the virtual swap pair as a baseline, the optimal swap pair among the swap pair set that induces the maximum reduction on clustering cost will never be removed (steps 12-19 of Algorithm IT].
Lemma D.5. Let gi (i E (1,23) be the swap groups constructed in step 12 of Algorithm 11 For each swap group gi. I = arg mi,ea, (P, (v) U (u)) be the swap pair in gi with minimum clustering cost. Then, t; will never be removed during the subsequent steps 14-19 of Algorithm 11
28

---

Proof. Given a swap pair E o (where o is the set of swap pairs constructed in step 9 of Algorithml| IT. we use cost(a) 8(P. CMf U {u}) to denote the true clustering cost induced by the swap pair a. where a = (u, v). According to Lemma D.4, by taking union bound over the success probability, we can assume without loss of generality that the clustering cost estimation (.) for each swap pair 0g E o should satisfy (o;) E [(1 1 /4000))cos((o)) ), (1 + + 1/400kk)
by setting = 1/400k. Consider an arbitrary swap group 9i for some i E 2. Let gm, be the swap pair in 9: with minimum clustering cost estimations (step 13 of Algorithm I. According to the judging criteria (step 15 and step 18 of Algorithm IT swap pair g E gi is removed only if its cost (g.()..) For the best swap pair t,, it hods that 400k I -
estimation (g) satisfies (ST
400k-1
1004++1
I (ST((m.)) ). where the first and third inequalities follow from Lemma D.4, the second inequality follows from that t; is the swap pair with minimum clustering cost. This indicates that the optimal swap pair will never be removed from the set 9; of candidate swap pairs.
Et TOUR+II I cost(t,) (1+1/40000 1/00/)) 100k+I
400k-1
Lemma D.5 implies that the optimal swap pair in each swap group gi E (1,23) will never be removed. Then, we will demonstrate that the clustering cost can be reduced significantly with certain probability for each local search step. According to Lemma B.3, with certain probability, for single local search step. there exists at least one swap pair (u, v) E o such that the clustering cost can be reduced by a factor of at least (1 - 1/100k). The following lemma shows that if such swap pair exists, the virtual swap pair and swap pairs with clustering cost larger than the virtual swap pair must be removed.
Lemma D.6. If there exists at least one swap pair (u, v) E 9i such that o(P, C.P.CMHT U (uf) S (1-1/100k)0(P) C), then swap pairs in gi with clustering cost larger than 8(P, C) will be removed.
Proof. Given a swap pair 1; = (u, v) such that 5(P, (PP,C\tt)U s 1/100k) 66 C), it holds that
ES S (1+ + 10004000)))s((((
<1 + 1/400k) (1 -1- (111/100))8(PPC)
S (1 3/ 400k))
where the first inequality follows from Lemma D.4, and the second inequality follows from the
definition of 4. For any swap pair g y: with clustering cost cost(g) 2 8(P, C),, it holds that
400k: 1
400k 1
EST(g)-
Z cost(g) (1 - 1/4004)) 400k: + 1
400k +
8(P.C) (400k- 1)²
400k(400k + 1)
>
101-3340000 3/400k).
where the second inequality follows from the definition for g. This indicates that y will be removed from dr according to the judging criteria.
By combining Lemma D.5 and Lemma D.6, the clustering cost will not increase during each local search step. Then, based on Lemma B.3, in each local search step, with certain probability, there exists at least one swap pair (u, v) E o such that the clustering cost can be reduced by a factor of at least (1 - 1/100k). By grouping such (u, v) with the virtual swap pair, the virtual swap pair will be removed from 9; for each i E {1. 23. Since the final swap pair is randomly picked from gi with [g.1 = 1, we can get that with certain probability, the clustering cost can be reduced by factor of at least (1 - 10o% for each local search step.
Corollary D.4. For the LS-KMedian algorithm, either with constant probability, there exists at least one swap pair such that the clustering cost can be reduced by a factor of (1 -1- 11/100),, or with probability
there exists at least one swap pair such that the clustering cost can be reduced by a factor of (1 - 1/100X) where X denotes the number of lonely centers (1 < < k -
Putting all these together, Theorem 3 can be proved.
Theorem 3. For the k-median problem, there exists an algorithm that can output constant approximate solution in expected Ö(nd + vnk²) running time without any data assumptions.
29

---

Proof. We first prove the approximation guarantees on Clustering quality. According to Corollary D.4, with certain probability, the clustering cost of P with respect to C' can be reduced by a certain fraction for each local search step. Then, similar to the analysis for Lemma B.4, after O(k)) local search steps, we can obtain a constant approximation in expectation.
As for the running time, since A can be compressed to A = poly(n, d) using the method proposed in [8]. the running time for initialization becomes ((nd). Each weighted sampling step requires to take Ö(v samples for = 1/400k. During this process, the Adaptive-Sampling algorithm used to take samples, where each sample can be picked within time Ō(k). The weighted sum calculations (step 13 of Algorithm 12) involves distance calculation between the sampled point to the centers after swaps, which takes O(k) time. Hence, the overall running time for weighted sampling during each local search step can be bounded by O( D(vnk)) since d can be compressed to (log with (mdd) preprocessing time. Finally, since Ō(k) rounds of local search steps are required to return constant approximate solution, the total running time is O(nd + +vkk²))
E Complementary Experiments
In this section, we present complementary experiments on both large-scale and small datasets. We also evaluate the robustness of the BanditFastl.S S algorithm by testing it with varying parameters.
E.1 Comparison results other large-scale datasets
Table 3 presents the comparison results with other local search algorithms (with mini-batch k-means++ as baseline) on datasets SYN, USSC_1900 and SUSY, with l varying from 10 to 100. The results demonstrate that our algorithm still outperforms other local search methods in terms of running time, while achieving much better clustering quality than the mini-batch k-means++ sceding method.
Table 3: Comparison results on datasets SYN, USS__1991 and SUSY with varying k where algorithms
fail to return a feasible solution within 24 hours not included
Method
Daniee
Cost
Tume()
Dataset
k
Cost
Time!s)
Dataeet
Cost
Timels)
3.250EE-0711.000 DB +05
3.85
20226 +08 26-06
2.30
6.071 442 11
28 404
0.59
MLS
7B-05
172.2
2.39
023-0802222200 A6-00
129.5
5,6302E+18
27.64
LSDS++
SUSY 10
440, 25
USC 1960
10
3111 411-08 3::300000000000
508
SYN
10
3,6000HH11 1:2 (0
50.19
LSa+
3.0820E-07H
$49.06
2.36508: 2.90
S19.16
5.62268441
-92
Outs
30147E40752.300
3,67
1.12
6.7609E405
11:70-06
MLS
2. 2.5095E-07:.5. 1B404
1899004-00 -08
1.7812E405:5
LSDS++
SUBY
20
0.5050E-07.2..00 0.7.222555-00
616.36
18CC 1990
20
1.949 1.1
372.19
SYN
20
1,041E+ (1:11111::111.
LSee
2. 2.4979
1308,59
18334E 1.22
552.91
1.5655K 13:3--44-0
Ours
2. 9.503 F-00F-12---00
2.68
1862 1F-088
+11F5406
1.24
[ 10233-405-70 9F-m
0.52
R-mecunes MLS 1 +538++
2,4013E-072...0 48404
8.06
1.66KE +08
4.03
3.22/2E404::2:30
2.05
3 3.22992-47:22 -7713777::33
435.32 32
1. 1.5849E +08
294.56
62.41
1.Sm
SUSY
2.2337E.407+2 78-06
2.2366F--07 0R-04
1410.85
719.61
USCC
8. 1.5399 408 1.38406
54658-08-
|1-00
558.81
450.36
7. 7.543 56515--444+ 38400
118.18
91 19.09
Ours
5E+405
3.22
1. 1.5797E-08:
3.35466
1.69
7.7674B444::: 1.540
0.69
Fomcanees
2.0866F-477-9
2F-01
9.67
1.31675--0--
16406
5.69
MLS
1.954F=407+1 77E-44
612.43
493.50
4.523PP444
1.5SS+
SUSY
50
1-9512E407H-3 0.--00
USCC 1950
50
1225554082
575.61
SYN
50
2.827084044:1
141.38
1.9474H077
2083.05
592.95
2.
195776600
2,10
31991154044122
L.749E407844 2B400
21.54
11.56
MLS
1251.1
1.01000:-08811 11330-066
L151 46
144580-00221.. 18+0
195.48
LSUS8+
SUNY
100
16633--07112 85-05
1595.44
LSSC
1900
100
9199 N5-477 1---05
915.23
SYN
100
7. 7.36556001857 01-5566------244-000
211.39
1.
16 7F-44
29
9.61828--47 +3.2--
921 009
$14.95
Ouzs
6.54
9.9665E
3.58
L.77:E::0::2.1E000
1.79
Table 4 presents the comparison results on small datasets with sizes smaller than 50,000, where the number of clusters ranges from 3 to 10. The results indicate that all local search algorithms achieve fast running times, with the MLS and LS++ algorithms delivering better clustering quality.
E.1.1 Experiments with Varying Parameters
In this subsection, we present experiments on the BanditFastLS algorithm with varying parameters. The input parameters for the BanditFastL.S algorithm include the probability parameter v) and the number of Hutropolis-Hasting steps used to approximate the D2-Sampling distribution. As outlined in the main text, we fix the number of Metropolis- Hastings steps at 20 and set vj = 1/(2n²) to meet the success probability guarantees specified in theoretical analysis. However, since controlling 77 to
30

---

Table 4: Comparison results small datasets
Timms
MLS5
707107 A (903)
460003 (UP41
01001
3300
22.577 9300
0L285 A DAETA
502045 à 14440
(2306 L 20001
B:126
212-33 928
02290 A 0J157
S0661 1998
IC19 1001
19605
2119
11 1996
100 000 194
R0B31 BARML
16:009
8.80
70:00:000 1770
LS05
2943379 1.48s
IJE33 D6L80
3802541
122241
01840
0004
2350977
44482
02125 à 00450
2943287 0.031
(3541 D0856
401 2407 0.500
04800
8 3233
346637 63763
0.890 à 0.134
MLA
291 13345
WIEES 0003)
412 22222711
016/9 3155
2130055 23031
UUISS 650
1Sen
B3253 8333
$4.366 183366-- 83208
0.102
3276
3324T18 14851
RI9-1 01351
On
626 1571 593861
0.300 60065
(55.600
467000
0:155
3556
544 0011 91966
AMM22 + 6A6AS
CARRE-10 $75.09
64887
355a60
6.667
23155-13 9.10m
6A215
0:129
LSUS
211100-00 335-25
5.81TE-00
345-08
12212
33995-12
0.170
8.15000110 185-25
1/236 03344
55118-00
155
132-42
DARSA
IJ833 RARE3
11993.500
B8631
8.991
3560197.550 a588 166222831
2000 1545
15503500000
1563 5191
0331
3015
2557000 5400 -02288851
02333 63372
1797780900
3822.900
120223765
DOBES
Backatiblle
1X123
& 0.2122
153601.266
06429
354220000
207223-772
6330
33555
8250
276750601
M.S
-1977220 2011
4972841
6381 396
(1241 110066
11710773
11327 0160N
LIDAN 1487 8357266
1536 140
3673
L-NES 9 UOONS
Da
102-200
6293 60037
6301--331 03157
HONSAN9 13519847
ESWOLLIIO 5535315
BORS7
B8745
71992784
6.10% 0.813
8,886 8.332
1.05
EA 177.
1175-02560 327160
5030 340
38-3301
68831
0.77
383.230
7AANT & 04895
174 5440 14.544
11 11.390
33852399 46411155
2403 340 2 212305
3144 PANL: (18-111
balance the efficiency of the algorithm is impractical, we instead manage the bandit process by fixing the maximum sample size at 50,000.
Table 5 presents the results for datasets SYN. USC_,990, and SUSY, with fixed number of clusters k = 50 and a maximum sample size of 50,000, while the Metropolis- Hastings steps vary from 20 to 100. The results indicate that increasing the number of Metropolis-Haann the overall clustering quality, with also a slight increase in the running time. However, the overall clustering performance remains consistent, demonstrating that the algorithm is robust variations in the number of Metropolis-Hastine steps. This indicates that fixing the number of transition steps at 20 sufficient to achieve high-quality clustering results with fast running time.
steps can improve
Table 5: Results with varying Eatropolis-Hastingg steps on datasets SYN, USC_1990, and SUSY,
where k is fixed at 50 and the maximum sample size is fixed at 50,000.
Method
Datasets
Steps (Transition)
Cost
Time(s)
20
3.0940F:-0442.1BB-000
1.42
40
3.14899=00422111-0
1.46
Ours
SYN
60
3.1231E+04+1.11443
1.71
80
3.1253E+0411.1B-03
1.73
100
3.0834E+0418.9E--2
1.88
20
1.2757E408+113344006
3.78
40
1.29996-08:33355--66
5.47
Ours
DSC__190
60
1.2771E+08+1.5E-06
4.72
80
1.2893E+08:334E5006
4.96
100
1.2740=+08+1.6E-006
4.35
20
1.9573E+07+2.2B--55
3.93
40
1.9526E+07+74B--0
3.75
Ours
SUSY
60
1.9615E++07+2285--5
3.95
80
1.9621E+07+2.4B400
3.91
100
1.9466E+07#1100--05
4.26
Table 6 presents the results for datasets SYN, USC_ 1990, and SUSY with a fixed number of clusters
k = 50 and 20 Metropolis-Hasting steps, while the maximum sample size varies from 50,000 to 250,000 during the bandit process. The results show that increasing the sample size improves the clustering quality, with also slight increase in the running time. However, the overall clustering
31

---

performance remains stable, indicating that the algorithm is robust to different choices of maximum sample sizes. Fixing the sample size at 50,000 is sufficient achieve high-qualig, clustering solutions with fast running time.
Table Results with varying maximum sample size on datasets SYN, USC_1990, and SUSY, where step is fixed at 20.
k is fixed at 50 and the
Method
Datasets
Maximum Samples
Cost
Time(s)
50,000
3. 1.404E+0411.EE-33
1.38
100,000
3.2266==0442245403
1.48
Ours
SYN
150,000
3.0020E+04:7.99E-22
1.79
200,000
3.0758E=04+1.9E503
2.18
250,000
3.0456E+0411.3B-003
2.17
50,000
L.2790=-88:32EE-06
4.55
100,000
1.2734F:40+
4.14
Ours
USC_1990
150,000
1.2756E--08+11.6----
3.98
200,000
1.2761E+008+114E--00
4.73
250,000
1.2859E+08+1.8E-06
4.23
50,000
1.9563E-08:11.66-055
3.57
100,000
1.9576E-0881.77--05
3.97
Ours
SUSY
150,000
1.9500E=-08:17444444
4.33
200,000
1.9713E-48422.66-05
3.61
250,000
1.9479E+0814.9E444
5.20
E.2
Comparison with Coresets Methods
In this subsection, give additional experiments by comparing proposed BanditFastLS algorithm with coresets methods.
Algorithms. We compare our BanditFastL.S algorithm against 2 different fast coresets construction methods, which are summarized as follows.
Lightweight Coreset: The lightweight coreset is constructed by sampling data points using mixture
[5]. The sensitivity is calculated by s(p) = P/ + ()) (P, μ), where P is the given dataset, fs is the 1-means solution for P. This method runs in O(nd/e c) time but provides a weaker guarantee, where there is an e μ)) term of additive error for some E (0, 1).
of sampling distributions
Fast-Coreset: This is the coreset construction method proposed in 18 where aspect ratio and rejection sampling methods are used to construct the coresets. The proposed algorithm can run in ((nd) time for constructing coreset with (1 + apappoxxinaiionnn guarantee.
Parameter Settings. In all datasets, a preprocessing step used remove the data points that are co-located at the same position. Following the settings in [18), we default the coreset size as 40k and we only run the dimension-reduction step on datasets with high dimensions (d 2 50).
Experimental Setup. For a fair comparison, we report both the time taken to construct the coreset and the overall time required obtain the final clustering solutions. After constructing the coresets, we use the LS++ algorithm (25), which provides significantly better clustering quality than other local search algorithms, to finalize the clustering centers on the coresets. Each algorithm is executed 10 times, and we report the average clustering costs along with deviations and the average running time.
Datasts. We mainly compare our algorithm with coresets methods on large-scale datasets, as local search methods have already demonstrated fast and accurate clustering performance smaller datasets. Large-scale datasets include SYN (IM 2), USC_1990 (2.45M M x 68), SUSY (5M x 17) and HIGGS (1IM 27) from UCI Machine Learning Repository We also include larger dataset
SIFT (100M x 128)
'https:/arcciie..c ics hwwp:h/wwiiiioosccceee..mennnnnp http://corpus texmexxsetetxxttttteeeemmeeeeeee .risa....iiiiffff
32

---

Tables 7 and 8 present a comparison between coreset construction time and the final clustering time for our algorithm on large-scale datasets. The results show that our algorithm outperforms coreset methods, achieving 20x speedup on average compared to state-of-the-art coreset construction techniques.
Results.
Table 7: Comparison of the coreset construction time and the final clustering time of our algorithm on the SYN, USC__S_C00 1990, and SUSY datasets.
Method
Dataset
k
Time(s)
Dataset
Time(s)
Dataset
Time(s)
Fast-Coreset
407.41
1052.54
918.94
LightWeight
SYN
10
2.09
USC_1990
10
96.15
SUSY
10
54.28
Ours
0.41
1.12
3.67
Fast-Coreset LightWeight
445,09
1186.93
954.93
SYN
20
2.11
USC_ 1990
20
95.25
SUSY
20
55.13
Ours
0.52
1.24
2.65
Fast-Coreset
441.29
1 1199.86
917.64
LightWeight
SYN
30
2.09
USC_ 1990
30
94.91
SUSY
30
54.72
Ours
0.69
1.69
3,42
Fast-Coreset
598.33
1034.49
941.23
LightWeight
SYN
50
2.09
USC_1990
50
95.49
SUSY
50
54.52
Ours
0.95
2.10
4.22
524.78
1037.78
903.0 00
LightWeight
SYN
100
2.12
USC_ 1990
100
94.49
SUSY
100
54.53
Ours
1.79
3.58
6.54
Table 8: Comparison of coreset construction time and final clustering time for our algorithm on the datasets HIGGS and SIFT.
Method
Dataset
Time(s)
Dataset
k
Time(s)
Fast-Coreset
2311.19
24h
LightWeight
HIGGS
10
185.56
SIFT
10
7687.46
Ours
5.73
99.81
Fast-Coreset
2479.21
24h
LightWeight Ours
HIGGS
20
184.01
SIFT
20
7.31.09
5.73
101.33
Fast-Coreset LightWeight Ours
2193.73
24h
HIGGS
30
183.76
SIFT
30
7007.9
6.92
105.85
Fast-Coreset
2126.66
24h 6964.60
LeigtWeight Ours
HIGGS
50
184.81
SIFT
50
7.88
125.89
Fast-Coreset
2195.06
24h
LightWeight
HIGGS
100
184.92
SIFT
100
6965.
Ours
9.87
163.82
Tables 9 and 10 present comparison between coreset-based local search methods and our Bandit- FastLS algorithm on large-scale datasets. The results show that our algorithm not only delivers better clustering quality but also provides at least 10x speedup the rumning time. On average, algorithm achieves a 5% improvement on clustering quality, while the running time is 50 times faster than coreset-bee sareset-based local search methods.
33

---

Table 9: Comparison of the coreset construction methods and our algorithm on SYN, U99, and SUSY datasets using LSS+ as the solver on the coresets.
Method
Dotaset
Cost
Time's)
Datzset
Cost
Timets)
Dutaset
Timers)
Corcee
8.00ENE++ 6F-04
413.76
2.598408--5 18406
1070.35
3.218407+ 7F:05
957.71
LightWeight
SYN
8.53
USC 1990
2.5984 8.3-83
11895
SUSY
3,24E407:+4 4E:05
94.28
Ours
6, 10--4
0,41
2,40E,085 4.15:06
1.2
301E=07 1233005
3,67
CarceCarcsett LightWeight
2,0E445+9 SE-OR
451.46
1.978408=2244 277440
1200-99
2.67E407=1 1E:05
997.4)
1.72E+061 18-40-03
8.72
1990
1.92600::22 2.5846
12306
SUSY
2.6/ 188-05
18-005
103.83
Our
L62E448 +4.9
0,52
1.80F 1,3F-06
1.34
2.50E +2.2000---22--0
2,65
Fists Corcset
9.1SE404 4E-08
448.03
1.668408=222 7.00::5555055555
1226.02
2.37E407455 SE404
965.04
LightWeigtt
8.00F+04+2 4F-03
8,94
1990
1611+00 114406
127.45
SUSY
2,36E4077: 22-00
101.54
Qurs
7.76F+04+ 1.54403
0.69
1.57F -3,2
1,09
2.23R:07 1,5E405
3.42
Fast. Coreset
636.26
1.75 L4
1065 X7
2,07E40711 1E+05
996.34
LighWWiiht
3. 3.10F-04++9.5 SF-02
9.56
3114000 535406
138.99
SUSY
2,05F407++ IF-00
110.29
Qurs
0.95
1,28F 1.6::066
2,80
1.95:407=2224446
4.22
Past Coreset
S.00E-0337.44E
534.09
1092. 75
1.73E00788 9E-04
98129
LightWeigh
8.77E401+0 IF-02
11.12
1038408- 3.58406
161666 d
SUSY
1.71 5F404
128.02
Qurs
1,17E404+2 1E-03
1.79
1.98E447=1155--000
3.58
L63E-07=29E40
6.54
Table 10: Comparison of the coreset construction methods and our algorithm on HIGGS and SIFT datasets using ISS+ as the solver on the coresets.
Method
Dataset
k
Cost
Time(s)
Dutaset
k
Cost
Time(s)
Fast-Coreset
2422.77
24h
LightWeight
HIGGS
10
1.6548E
9E+05
297.14
10
1.1444
8853,91
Ours
L.5601E
5.73
99,81
Fast-Coreset LightWeight Ours
14994E+088555
2587.18
24h
HIGGS
20
+. 2F+05
282.25
SIFT
20
1.0161B++3+1144++10
8470.58
1.4335E:008 +5,2
5.73
9.6759E:+1-+16677-111
101.33
Fast-Coreset
2310.14
>: 24h
LightWeight
HIGGS
30
1.41177--0:22.7E--000
291.58
SIFT
30
9:6
8586.20
Ours
1:7
105.85
Fast-Coreset LightWeight Ours
E.313:E 1E+05
2263.74
24h
HIGGS
50
1.3 +.31335==0::55.11 1E+05
316.89
SIFT
50
9. 1200E+1247.1E+000
9267.05
1.2636E=48:3333EEM
7.88
8.
125.89
Fast-Coreset LightWeight
1.1935E+08 08:-:30088-3...0 2E+05
2398.43
>. 24h
HIGGS
100
1
373.28
SIFT
100
8.4956E+ 8E
10891.71
Ours
1.1434E:08 :1.3:
9.87
8.1365T::12:::4999
163.82
NeurIPS Paper Checklist
1. Claims
Question: Do the main claims made in the abstract and introduction accurately reflect the paper's contributions and scope?
Answer: [Yes]
Justification: The main claims made in the abstract and introduction can reflect the paper's contributions and scope.
Guidelines:
The answer NA means that the abstract and introduction do not include the claims made the paper.
The abstract and/or introduction should clearly state the claims made, including the contributions made the paper and important assumptions and limitations. A No or
NA answer to this question will not be perceived well by the reviewers.
The claims made should match theoretical and experimental results, and reflect how much the results can be expected to generalize to other settings.
It is fine to include aspirational goals as motivation as long as it is clear that these goals are not attained by the paper.
2. Limitations
Question: Does the paper discuss the limitations of the work performed by the authors?
Answer: [Yes]
Justification: The limitations of the work have been discussed in the Conclusion part.
34

---

Guidelines:
The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper.
The authors are encouraged to create a separate "Limitations" section their paper.
The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be.
The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with few runs. In general, empirical results often depend on implicit assumptions, which should be articulated.
The authors should reflect on the factors that influence the performance of the approach. For example, facial recognition algorithm may perform poorly when image resolution low or images are taken in low lighting. Or speech-to-tex system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon.
The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size.
If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairesss.
While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an impor- tant role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations.
3. Theory assumptions and proofs
Question: For each theoretical result, does the paper provide the full set of assumptions and complete (and correct) proof?
Answer: [Yes]
Justification:
The assumptions have been clearly stated in the main context.
Guidelines:
The answer NA means that the paper does not include theoretical results.
All the theorems, formulas, and proofs in the paper should be numbered and cross- referenced.
All assumptions should be clearly stated or referenced the statement of any theorems.
The proofs can either appear in the main paper or the supplemental material, but if they appear in the supplemental material, the authors are encouraged to provide short proof sketch to provide intuition.
Inversely, any informal proof provided in the core of the paper should be complemented by formal proofs provided in appendix or supplemental material.
Theorems and Lemmas that the proof relies upon should be properly referenced.
4. Experimental result reproducibility
Question: Does the paper fully disclose all the information needed to reproduce the main ex- perimental results of the paper to the extent that affects the main claims and/or conclusions of the paper (regardless of whether the code and data are provided or not)?
Answer: [Yes]
Justification: All implementation details, parameter settings, and datasets are fully disclosed ensure the experiments can be reliably reproduced.
Guidelines:
The answer NA means that the paper does not include experiments.
35

---

If the paper includes experiments, a No answer to this question will not be perceived well by the reviewers: Making the paper reproducible is important, regardless of whether the code and data are provided or not.
If the contribution is dataset and/or model, the authors should describe the steps taken to make their results reproducible or verifiable.
Depending on the contribution, reproducibility can be accomplished in various ways. For example, if the contribution is a novel architecture, describing the architecture fully might suffice, or if the contribution is specific model and empirical evaluation, it may be necessary to either make possible for others to replicate the model with the same dataset, or provide access to the model. In general. releasing code and data is often one good way to accomplish this, but reproducibility can also be provided via detailed instructions for how to replicate the results, access to a hosted model (e.g., in the case of a large language model), releasing of a model checkpoint, or other means that are appropriate to the research performed.
While NeurIPS does not require releasing code, the conference does require all submis- sions to provide some reasonable avenue for reproducibility, which may depend on the nature of the contribution. For example
(a) the contribution is primarily a new algorithm, the paper should make clear how reproduce that algorithm. (b) If the contribution is primarily new model architecture, the paper should describe the architecture clearly and fully.
(c) If the contribution is new model (e.g., a large language model), then there should either be way to access this model for reproducing the results or way to reproduce the model (e.g.. with an open-source dataset instructions for how to construct the dataset).
(d) We recognize that reproducibility may be tricky in some cases, in which case authors welcome describe the particular way they provide for reproducibilityy In the case of closed-source models, it may be that access to the model is limited in some way (e.g., registered users), but it should be possible for other researchers to have some path to reproducing or verifying the results.
5. Open access to data and code
Question: Does the paper provide open access to the data and code, with sufficient instruc- tions to faithfully reproduce the main experimental results, as described in supplemental material?
Answer: [Yes]
Justification: The paper provides open access to both data and code, along with detailed instructions in the supplemental material for faithful reproduction of the results.
Guidelines:
The answer NA means that paper does not include experiments requiring code.
Please see the NeurlPS code and data submission guidelines (https://nips. cc/
for more details.
While we encourage the release of code and data, we understand that this might not be possible, so "No" is an acceptable answer. Papers cannot be rejected simply for not including code, unless this is central to the contribution (c.g., for a new open-source benchmark).
The instructions should contain the exact command and environment needed to run to reproduce the results. See the NeurIPS code and data submission guidelines (https
//nips.cc/public/eee. www.ddddeeeooeuummmmccco ssi ssionoolicy) for more details.
The authors should provide instructions on data access and preparation, including how to access the raw data, preprocessed data, intermediate data, and generated data, etc.
The authors should provide scripts to reproduce all experimental results for the new proposed method and baselines. If only subset of experiments are reproducible, they should state which ones are omitted from the script and why.
At submission time, to preserve anonymity, the authors should release anonymized versions (if applicable).
36

---

Providing as much information as possible in supplemental material (appended to the paper) recommended, but including URLs to data and code permitted.
6. Experimental setting//etails
Question: Does the paper specify all the training and test details (e.g., data splits, hyper- parameters, how they were chosen, type of optimizer, etc.) necessary to understand the results?
Answer: [Yes]
Justification: This paper specifies all datasets details, including data preprocessing, hy- perparameters, selection procedures, and optimizer settings, ensuring the results are fully understandable and reproducible.
Guidelines:
The answer NA means that the paper does not include experiments.
The experimental setting should be presented in the core of the paper to level of detail that is necessary to appreciate the results and make sense of them.
The full details can be provided either with the code, in appendix, or as supplemental material.
7. Experiment statistical significance
Does the paper report error bars suitably and correctly defined or other appropriate information about the statistical significance of the experiments?
Question:
Answer: [Yes]
Justification: This paper reports standard deviations over multiple runs indicate the statistical significance of the results.
Guidelines:
The answer NA means that the paper does not include experiments.
The authors should answer "Yes" if the results are accompanied by error bars, confi- dence intervals, or statistical significance tests, least for the experiments that support the main claims of the paper.
The factors of variability that the error bars capturing should be clearly stated (for example, train/test split, initialization, random drawing of some parameter, or overall
with given experimental conditions).
The method for calculating the error bars should be explained (closed form formula, call to a library function, bootstrap, etc.)
The assumptions made should be given (e.g., Normally distributed errors).
It should be clear whether the error bar is the standard deviation the standard error of the mean.
It is OK to report I-sigma еrror bars, but one should state it. The authors should preferably report a 2-sigma error bar than state that they have 96% CI, if the hypothesis of Normality of errors is not verified.
For asymmetric distributions, the authors should be careful figures symmetric error bars that would yield results that are out of range (e.g. negative error rates).
show in tables or
If error bars are reported in tables or plots, The authors should explain the text how they were calculated and reference the corresponding figures tables in the text.
8. Experiments compute resources
For each experiment, does the paper provide sufficient information on the com- puter resources (type of compute workers, memory, time of execution) needed reproduce the experiments?
Question:
Answer: [Yes]
Justification: This paper provides sufficient information on computer resources needed to reproduce the experiments.
Guidelines:
37

---

The answer NA means that the paper does not include experiments.
The paper should indicate the type of compute workers CPU or GPU, internal cluster, or cloud provider, including relevant memory and storage.
The paper should provide the amount compute required for each of the individual experimental runs well estimate the total compute.
The paper should disclose whether the full research project required more compute than the experiments reported in the paper (e.g., preliminary or failed experiments that didn't make into the paper).
9. Code of ethics
Question: Does the research conducted in the paper conform, in every respect, with the NeurIPS Code of Ethics http::/ps///aarrg...ooooooo///aa
Answer: [Yes]
Justification: The research conducted in the paper conforms, in every respect, with the NeurIPS Code of Ethics.
Guidelines:
The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.
If the authors answer No, they should explain the special circumstances that require deviation from the Code of Ethics.
The authors should make sure to preserve anonymity (e-g., there a special consid- eration due to laws or regulations in their jurisdiction).
10. Broader impacts
Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?
Answer: [Yes]
Justification: The potential positive societal impacts and negative societal impacts have been discussed in the Conclusion part.
Guidelines:
The answer NA means that there no socictal impact of the work performed.
If the authors answer NA or No, they should explain why their work has no societal impact why the paper does not address societal impact.
Examples of negative societal impacts include potential malicious unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g.. deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.
The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster.
The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology being used intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology.
If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition attacks, mechanisms for monitoring misuse, mechanisms to monitor how system learns from feedback over time, improving the efficiency and accessibility of ML.)
11. Safeguards
Question: Does the paper describe safeguards that have been put in place for responsible release of data models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?
38

---

Answer: [NA]
Justification: The paper does not involve pretrained language models, image generators, or any datasets or models with potential misuse risks.
Guidelines:
The answer NA means that the paper poses such risks.
Released models that have high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere usage guidelines or restrictions to access the model or implementing safety filters.
Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images.
We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make best faith effort.
12. Licenses for existing assets
Question: Are the creators or original owners of assets (e.g.. code, data, models), used in the paper. properly credited and are the license and terms of use explicitly mentioned and properly respected?
Answer: [Yes]
Justification: All external assets, including code and datasets, are properly credited, and their licenses and terms of use explicitly acknowledged and respected.
Guidelines:
The answer NA means that the paper does not use existing assets.
The authors should cite the original paper that produced the code package dataset.
The authors should state which version of the asset used and, possible, include URL.
The name of the license (e.g., CC-BY 4.0) should be included for each asset.
For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided.
If assets released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, pwpppsaittcoooeooo thcode. com datasets has curated licenses for some datasets. Their licensing guide can help determine the license dataset.
For existing datasets that are re-packaged, both the original license and the license of derived asset (if has changed) should be provided. this information is not available online, the authors are encouraged to reach out to the asset's creators.
13. New assets
Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets?
Answer: [Yes]
Justification:
This paper does not release new assets.
Guidelines:
The answer NA means that the paper does release assets.
Researchers should communicate the details the dataset/code/moodll as part their submissions via structured templates. This includes details about training, license, limitations, etc.
The paper should discuss whether and how consent was obtained from people whose asset used.
At submission time, remember to anonymize your assets (if applicable). You can either create anonymized URL or include an anonymized zip file.
39

---

14. Crowdsourcing and research with human subjects
Question: For crowdsourcing experiments and research with human subjects, does the paper include the full text of instructions given participants and screenshots, if applicable, as well as details about compensation (if any)?
Answer: [NA]
Justification:
This paper does involve crowdsoureing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Including this information the supplemental material is fine, but if the main contribu- tion of the paper involves human subjects, then much detail as possible should be included the main paper.
According to the NeurIPS Code of Ethics, workers involved in data collection, curation, or other labor should be paid at least the minimum wage in the country of the data collector.
15. Institutional review board (IRB) approvals or equivalent for research with human subjects
Question: Does the paper describe potential risks incurred by study participants, whether such risks were disclosed to the subjects, and whether Institutional Review Board (IRB) approvals (or equivalent approval/teview based on the requirements of your country institution) were obtained?
Answer: [NA]
Justification:
The paper does not involve crowdsoureing nor research with human subjects.
Guidelines:
The answer NA means that the paper does not involve crowdsourcing nor research with human subjects.
Depending on the country in which research is conducted, IRB approval (or equivalent) may be required for any human subjects research. If you obtained IRB approval, you should clearly state this in the paper.
We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere the NeurIPS Code of Ethics and the guidelines for their institution.
For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review.
16. Declaration of LLM usage
Question: Does the paper describe the usage of LLMs if it is important, original, or non-standard component of the core methods in this research? Note that if the LLM is used only for writing, editing, formatting purposes and does not impact the core methodology, scientific rigorousness, or originality of the research, declaration is not required.
Answer: [NA]
Justification: The core method development in this research does involve LLMs as any important, original, or non-standard components.
Guidelines:
The answer NA means that the core method development in this research does not involve LLMs as any important, original, non-standard components.
Please refer to L.MM policy autee:://berrrppsccc///ofrreeerp for what should or should not be described.
40


Reviews:
[{"review_id": "rykpR33Erd", "reviewer": "Reviewer_jiTA", "summary": "The paper proposes faster local search algorithms for k-means and k-median clustering using adaptive sampling and bandit strategies, claiming significant runtime improvements while preserving constant-factor approximation guarantees. The authors introduce tree-based data structures to accelerate sampling and model swap selection as a multi-armed bandit problem.", "questions": "\"The authors mention that 'to ensure fairness, all algorithms perform 400 local search steps' (line 329). However, for large-scale datasets, such a high number of steps is often unnecessary in practice. Can the authors report the performance of each method with fewer steps, such as 50 or 100, to provide a more realistic comparison of practical efficiency?", "limitations": "The claimed theoretical guarantees depend heavily on assumptions such as sub-Gaussian priors on swap pairs and specific distributional behavior of optimal clusters (e.g., exponential tails). These assumptions are not standard or justified empirically", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 2, "strengths_and_weaknesses": "Strengths: The use of tree-based structures and bandit modeling is conceptually interesting.\nWeaknesses :\n    Over-reliance on strong assumptions:\n    The claimed theoretical guarantees depend heavily on assumptions such as sub-Gaussian priors on swap pairs and specific distributional behavior of optimal clusters (e.g., exponential tails). These assumptions are not standard or justified empirically, and no sensitivity analysis is provided to show how the algorithm performs when these assumptions are violated.\nPresentation:\n    The paper has typos and formatting issues.\n    The title itself has a typo: “Adpative” instead of “Adaptive, which is unusual.\nExperimental results: \nExperimental results show that the clustering quality is slightly weaker compared to MLS, which is expected since multi-swap strategies can naturally achieve better results than single-swap approaches.\nAll methods evaluated in the experiments show only marginal improvements in clustering quality over k-means++, which is unsurprising given that k-means++ already performs well when the clusters are well-separated. In such cases, local search via center swaps is not particularly effective for further improving the clustering quality", "was_revised": true, "final_justification": "The  rebuttal resolving circular reasoning concerns. However, some issue remain: (1) the sub-Gaussian assumption is theoretically elegant but impractical to verify in real applications; (2) no improved runtime bounds are provided despite empirical speedups;  (3) presentation issues persist. At the same time, the theoretical framework and synthetic validation provide some foundation for borderline paper.", "rebuttal": "We thank the reviewer for the thoughtful comments and constructive feedbacks. Below, we address the concerns.\n\n**Weakness 1: These assumptions are not standard or justified empirically, and no sensitivity analysis is provided to show how the algorithm performs when these assumptions are violated.**\n\nResponse: We thank the reviewer for raising this concern. For the $k$-median problem, we obtain theoretical guarantees without assumptions by applying weighted sampling on $\\ell_1$-norm objectives. However, this approach cannot be used for $k$-means due to the large approximation loss caused by squared distance definitions. To reduce the dependence on both $n$ and $k$ in the runtime, we apply bandit method for best swap pair selection. As shown in the literature (e.g., [6], [30], [31] in our ref list), to obtain theoretical guarantees in bandit settings, sub-Gaussianity is a standard assumption. \n\nPrevious assumptions for bandit methods are much stronger than ours. In [6], [30], and [31], it is assumed that the gain of each arm follows a fixed distribution with mean $\\gamma$ and variance 1, or the variance-to-maximum-gain ratio is bounded by a constant. In contrast, we only require the variance to be bounded by a function of the average cost, making our assumption much weaker.\n\nTo further improve runtime for $k$-means, we design a practical algorithm (Alg 3) based on the Metropolis-Hastings method. To achieve theoretical guarantees, both exponential tail and non-degeneracy assumptions are required in previous work (e.g., [3] in the ref list). In contrast, our analysis only requires an exponential tail and a cluster size condition, making it much weaker than previous ones.\n\nTo validate the assumptions, we conduct the following experiments. For the sub-Gaussianity assumption, we use the fact that a variable with range $2\\sigma$ is $\\sigma$-sub-Gaussian, and measure the ratio between gain range and average clustering cost $\\delta^2(P,C)/n$. The gain for a point $p$ under swap $(u, v)$  and center set $C$ is defined as $\\delta^2(p, C \\setminus u \\cup v) - \\delta^2(p, C)$. In Tables 1 and 2, we report the mean, variance, and maximum of the ratios calculated during the local search processes. The results show that this ratio is typically bounded by a constant (< 50) across most datasets, verifying the sub-Gaussianity assumption. Even on datasets like SYN, where the assumption may not hold (with large ratio), our algorithm still achieves over 20× speedup, demonstrating the empirical robustness.\n\nFor the exponential tail assumption, we evaluate it using approximate solutions from multi-swap local search, since the optimal clusters are unknown. Similar evaluation based on approximate solutions has also been adopted in prior work such as [1]. We rely on the fact that the log-CCDF (Complementary Cumulative Distribution Function) of an exponential distribution is upper bounded by a linear function of the variable. For each approximate cluster, we: (1) compute the log-CCDF of squared distances, and (2) apply linear regression using slope < 0 and p-value < 0.05 as criteria. As shown in Tables 3-4, 7 out of 12 datasets are likely to satisfy the assumption. All datasets with over 1M points meet the condition, while smaller ones (e.g., Glass) may not. Most datasets also satisfy the cluster size condition, though SYN ($k=100$) may not for its small clusters. However, our algorithm achieves over 20× speedup on datasets violating both sub-Gaussian and cluster size assumptions (e.g., SYN), demonstrating the empirical robustness.\n\n[1] Huang et al. The power of uniform sampling for k-median. Proc. ICML 2023, pp. 13933-13956.\n\n*Table 1: Ratio of gain range to average cost on large datasets*\n\n|**Dataset**|**$k$**|**Mean**|**Variance**|**Max**|\n|---|---|---|---|---|\n|SYN|10|9.75|11.89|77.92|\n||30|7.97|18.98|211.73|\n||50|6.54|19.50|185.31|\n||100|6.23|25.89|370.32|\n|USC1990|10|3.77|7.96|41.67|\n||30|1.71|2.81|48.99|\n||50|1.69|1.33|12.31|\n||100|1.62|1.17|11.00|\n|SUSY|10|4.05|4.83|34.08|\n||30|2.99|4.51|33.96|\n||50|2.87|3.74|28.71|\n||100|2.41|3.74|42.28|\n|HIGGS|10|1.66|1.14|12.82|\n||30|1.56|1.25|14.94|\n||50|1.41|1.08|22.56|\n||100|1.41|2.04|41.22|\n\n*Table 2: Ratio of gain range to average cost on small datasets*\n\n|**Dataset**|**$k$**|**Mean**|**Variance**|**Max**|\n|---|---|---|---|---|\n|iris|3|4.98|3.32|13.61|\n||5|3.54|3.16|15.42|\n||10|2.36|1.21|7.02|\n|seeds|3|3.84|2.39|11.40|\n||5|2.99|1.46|8.41|\n||10|2.66|1.45|9.81|\n|glass|3|4.31|3.15|16.21|\n||5|3.65|3.56|14.74|\n||10|3.01|3.12|15.88|\n|Who_FL|3|5.37|4.74|32.04|\n||5|6.19|7.74|37.71|\n||10|4.67|8.42|59.39|\n|HCV_L|3|7.21|13.31|62.19|\n||5|8.34|19.14|96.74|\n||10|3.68|8.22|78.61|\n|TRR_FL|3|1.31|0.34|2.53|\n||5|1.36|0.41|2.75|\n||10|1.24|0.38|2.73|\n|urbanGB_L|3|4.89|3.11|14.68|\n||5|5.47|4.22|22.74|\n||10|5.37|5.16|28.84|\n\n*Table 3: Cluster size and tail tests on large datasets*\n\n|**Dataset**|**k**|**MinClusterSize/$k^2$**|**Num of Exponential Tail Clusters (slope<0, pvalue<0.05)**|\n|---|---|---|---|\n|SYN|10|437.75|10|\n||30|6.51|30|\n||50|1.43|50|\n||100|0.14|100|\n|USC1990|10|889.5|10|\n||30|13.52|30|\n||50|2.05|50|\n||100|0.34|100|\n|SUSY|10|652.66|10|\n||30|39.61|30|\n||50|6.31|50|\n||100|0.78|100|\n|HIGGS|10|8907.19|10|\n||30|143.87|30|\n||50|11.49|50|\n||100|1.82|100|\n\n*Table 4: Cluster size and tail tests on small datasets*\n\n|**Dataset**|**k**|**MinClusterSize/$k^2$**|**Num of Exponential Tail Clusters(slope<0, pvalue<0.05)**|\n|---|---|---|---|\n|iris|3|4.33|3|\n||5|0.52|5|\n||10|0.04|9|\n|seeds|3|6.77|3|\n||5|1.24|5|\n||10|0.11|10|\n|glass|3|2.44|1|\n||5|0.08|4|\n||10|0.02|7|\n|Who_FL|3|5.66|3|\n||5|0.24|4|\n||10|0.02|8|\n|HCV_L|3|1.22|3|\n||5|0.28|5|\n||10|0.03|9|\n|TRR_FL|3|119.33|3|\n||5|36.96|5|\n||10|2.05|10|\n|urbanGB_L|3|2449|3|\n||5|894.52|5|\n||10|65.14|10|\n\n**Weakness 2: Regarding the typos**\n\nResponse: We apologize for the typos. In the revised version, we will carefully correct them to improve the presentation.\n\n**Weakness 3: The runtime comparison is also somewhat unfair, as MLS performs multi-swap operations while LS++ is based on single swaps. There is no clear reason why MLS should be faster than LS++**\n\nResponse: We thank the reviewer for raising this concern.\n\nAlthough MLS is a multi-swap method and LS++ is based on single-swap, as shown in prior work (e.g., [7], [21] in the ref list), both of them are compared together to evaluate the trade-offs in runtime and clustering quality.\n\nWhile single-swap methods are generally expected to be faster than multi-swap methods, MLS is faster than LS++ because it employs a heuristic that estimates swap gains using a small sample of size $O(k \\log k)$ instead of the full dataset. In contrast,  LS++ evaluates exact gains on the entire dataset, incurring a cost of $O(ndk)$ per local search step. This makes MLS more efficient despite its multi-swap nature.\n\nIt can be seen that the sampling step is the key factor behind MLS’s speed advantage. Additional experiments show that when MLS is modified to compute exact gains over the full dataset (denoted as MLSP in Tables 5-7), it becomes much slower than LS++.\n\nWe will clarify this in our revised version.\n\n**Question 1: Can the authors report the performance of each method with fewer steps, such as 50 or 100, to provide a more realistic comparison of practical efficiency?**\n\nResponse: We thank the reviewer for this insightful suggestion. In our experiments, 400 local search steps is set to ensure that all methods have sufficient iterations to converge. In the following (Tables 5-7), we conducted additional experiments on large datasets SUSY (5M), HIGGS (11M) and SIFT (100M) using 50 local search steps. The results show that our algorithm still outperforms other local search methods with 50x to 300x speedups and comparable clustering quality. \n\n*Table 5: Results on SUSY with 50 LS steps*\n\n|**Dataset**|**$k$**|**Method**|**Cost**|**Time(s)**|\n|---|---|---|---|---|\n|SUSY|10|k-means++|3.35E+07$\\pm$1.0E+05|3.83|\n|||LS++|3.03E+07$\\pm$4.7E+05|150.51|\n|||LSDS++|3.27E+07$\\pm$7.4E+05|135.49|\n|||MLS|3.08E+07$\\pm$4.2E+05|33.85|\n|||MLSP|**2.94E+07$\\pm$6.3E+05**|6286|\n|||Ours|3.05E+07$\\pm$2.5E+05|**1.98**|\n|SUSY|50|k-means++|2.40E+07$\\pm$7.4E+05|8.06|\n|||LS++|1.94E+07$\\pm$1.8E+05|328.23|\n|||LSDS++|2.03E+07$\\pm$2.1E+05|256.04|\n|||MLS|1.95E+07$\\pm$4.5E+04|106.09|\n|||MLSP|-|>6h|\n|||Ours|**1.94E+07$\\pm$7.4E+04**|**3.83**|\n|SUSY|100|k-means++|1.74E+07$\\pm$4.2E+03|21.54|\n|||LS++|**1.63E+07$\\pm$1.0E+04**|695.99|\n|||LSDS++|1.71E+07$\\pm$9.3E+04|387.79|\n|||MLS|1.64E+07$\\pm$2.3E+04|191.41|\n|||MLSP|-|>6h|\n|||Ours|1.64E+07$\\pm$3.0E+04|**9.12**|\n\n*Table 6: Results on HIGGS with 50 LS steps*\n\n|**Dataset**|**k**|**Method**|**Cost**|**Time(s)**|\n|---|---|---|---|---|\n|HIGGS|10|k-means++|1.66E+08$\\pm$1.3E+06|**6.51**|\n|||LS++|1.58E+08$\\pm$2.3E+06|507.53|\n|||LSDS++|1.62E+08$\\pm$2.8E+06|304.46|\n|||MLS|1.56E+08$\\pm$1.7E+06|65.31|\n|||MLSP|**1.53E+08$\\pm$7.8E+05**|10266.83|\n|||Ours|1.56E+08$\\pm$3.4E+05|7.39|\n|HIGGS|50|k-means++|1.31E+08$\\pm$1.1E+05|20.58|\n|||LS++|**1.26E+08$\\pm$1.2E+05**|13234.24|\n|||LSDS++|1.30E+08$\\pm$1.2E+06|1007.52|\n|||MLS|1.26E+08$\\pm$2.1E+05|209.14|\n|||MLSP|-|>6h|\n|||Ours|1.26E+08$\\pm$1.8E+05|**18.23**|\n|HIGGS|100|k-means++|1.20E+08$\\pm$6.9E+05|40.98|\n|||LS++|-|>6h|\n|||LSDS++|1.19E+08$\\pm$1.5E+06|2137.09|\n|||MLS|1.15E+08$\\pm$6.8E+04|637.82|\n|||MLSP|-|>6h|\n|||Ours|**1.14E+08$\\pm$1.5E+05**|**25.82**|\n\n*Table 7: Results on SIFT with 50 LS steps*\n\n|**Dataset**|**k**|**Method**|**Cost**|**Time(s)**|\n|---|---|---|---|---|\n|SIFT|10|k-means++|1.12E+13$\\pm$2.1E+10|181.03|\n|||LS++|1.05E+13$\\pm$2.8E+10|23560.27|\n|||LSDS++|1.65E+13$\\pm$3.0E+10|6265.61|\n|||MLS|1.05E+13$\\pm$2.5E+10|4691.80|\n|||MLSP|-|>6h|\n|||Ours|**1.05E+13$\\pm$2.1E+09**|**70.24**|\n|SIFT|50|k-means++|9.12E+12$\\pm$1.1E+09|381.74|\n|||LS++|-|>6h|\n|||LSDS++|-|>6h|\n|||MLS|-|>6h|\n|||MLSP|-|>6h|\n|||Ours|**8.75E+12$\\pm$1.1E+09**|**122.51**|\n|SIFT|100|k-means++|8.48E+12$\\pm$5.4E+09|677.93|\n|||LS++|-|>6h|\n|||LSDS++|-|>6h|\n|||MLS|-|>6h|\n|||MLSP|-|>6h|\n|||Ours|**8.13E+12$\\pm$8.1E+08**|**65.51**|"}, {"review_id": "Sm16Mws5Yb", "reviewer": "Reviewer_sbwN", "summary": "This paper proposes a new local search method to accelerate clustering on large datasets. Under a sub-Gaussian prior assumption on \"swap pairs\", the authors propose an $O(1)$-approximation algorithm with $\\tilde{O}(nd+\\sqrt{n}k^4d)$-running time for $k$-Means (and $\\tilde{O}(nd+\\sqrt{n}k^3d)$-time for $k$-Median without assumptions). Additionally, they achieve a better running time $\\tilde{O}(nd+k^4d)$ for $k$-Means under more assumptions on the optimal center set. This is faster than existing local search methods.\n\nIn the algorithm, they first convert the dataset by a tree structure, and use an adaptive sampling method to approximate $D^2$-Sampling distribution. Then they use a bandit method to accelerate the swap pair identification. Experimental results demonstrate that the proposed algorithm performs better than existing methods on large datasets, achieving a significant speedup.", "questions": "Could the authors provide a comparison between their new algorithm and a faster coreset-based approach over one million data points? For instance, a possible method is to apply uniform sampling to obtain a weak coreset and run $k$-Means++ on the coreset. If additional experimental results could demonstrate that the proposed algorithm still outperforms such coreset-based methods, I think it will increase the contribution.", "limitations": "Yes", "rating": 5, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 2, "clarity": 3, "significance": 2, "originality": 3, "strengths_and_weaknesses": "Strengths:\n* Compared with the prior work, the algorithm breaks the $O(ndk)$ running time using local search method, which is non-trivial.\n* The adaptive sampling method is novel and useful. Combined with a tree structure, it avoids maintaining distances to all data points and reduces the running time.\n* The experiments are thorough and convincing, conducted on multiple datasets to validate the algorithm's acceleration performance. The proposed algorithm is practical for large datasets.\n\nWeaknesses\n* The Sub-Gaussian Prior assumption seems not mild and difficult to verify in practice. Generally, the assumptions of an algorithm are expected to depend on the dataset and be easy to verify. However, this assumption concerns the “swap pairs” generated during the algorithm's execution, making it difficult to determine whether a given dataset satisfies the assumption in practice.\n* As mentioned in the introduction of this paper, there exist coreset-based algorithms achieving $\\tilde{O}(nd + k^3d)$ running time, which is better than the new algorithm theoretically. This weakens the contribution.\n* \"For datasets with sizes over 100 million, the constructions for coresets can require several hours of computation\". This conclusion may have some problems. By [Huang et al., 2023], a weak coreset can be obtained by uniform sampling under the balanceness assumption, which is very efficient.\n* Empirically, the algorithm improves only on very large datasets with over one million data points and has no significant advantages on small datasets.\n\n[Huang et al., 2023] Huang L, Jiang S H C, Lou J. The power of uniform sampling for k-median[C]//International Conference on Machine Learning. PMLR, 2023: 13933-13956.", "was_revised": true, "final_justification": "I appreciated the additional experiments in the author's response. It justifies the most concern on the reasonability of the assumption. Although the assumption is still not easy to be checked before running the algorithm, it provides evidences that they could hold in practice, largely increasing the applicability of the algorithm together with its provably guarantee.", "rebuttal": "We thank the reviewer for the thoughtful comments and constructive feedbacks. Below, we address the concerns.\n\n**Weakness 1: The Sub-Gaussian Prior assumption seems not mild and difficult to verify in practice. Generally, the assumptions of an algorithm are expected to depend on the dataset and be easy to verify. However, this assumption concerns the “swap pairs” generated during the algorithm's execution, making it difficult to determine whether a given dataset satisfies the assumption in practice.**\n\nResponse: We thank the reviewer for raising this important concern. Unlike the $k$-median problem with $\\ell_1$-norm objective, our weighted sampling strategy (which requires no assumptions) cannot be applied to $k$-means due to approximation loss for swap gain estimation under squared distances. To overcome this issue, we model the swap pair selection process as a bandit problem. In this setting, the sub-Gaussianity assumption is both necessary for theoretical guarantees and commonly adopted in prior work (e.g., [6], [30], [31] in our ref list).\n\nWe refer to our sub-Gaussian Prior assumption as mild because it is much weaker than those used in prior work. Existing bandit-based clustering methods (e.g., [6], [30], [31]) typically assume that the gain of each swap pair follows a fixed distribution with mean $\\gamma$ and variance $1$, or that the ratio between the variance and the maximum gain is bounded by a constant. In contrast, we only require that the variance is bounded by a function of the average cost, making it a milder assumption than previous ones.\n\nTo verify our assumptions empirically, we conduct the following experiments. We use the fact that any variable with range $2\\sigma$ is $\\sigma$-sub-Gaussian. Thus, the ratio between the gain range and the average clustering cost can be used to measure whether the sub-Gaussian assumption is satisfied. For each swap pair $(u,v)$ and center set $C$, we define the gain for a point $p$ as $\\delta^2(p, C \\setminus u \\cup v) - \\delta^2(p, C)$. Then the ratio can be calculated as the gain divide by $\\delta^2(P, C)/n$. For each dataset, we run our algorithm 5 times with 400 local search steps and record the gain-range-to-average-cost ratio for each swap pair. Tables 1-2 give the mean, variance, and maximum of these ratios. The results show that the maximum ratio is typically bounded by a constant (less than 50) across most datasets, supporting the validity of the sub-Gaussian assumption. Even on datasets like SYN, where the assumption may not strictly hold, our algorithm still achieves over 20× speedup with competitive clustering quality.\n\n*Table 1: Ratio of gain range to the average cost on large datasets*\n\n|**Dataset**|**$k$**|**Mean**|**Variance**|**Max**|\n|---|---|---|---|---|\n|SYN|10|9.75|11.89|77.92|\n||30|7.97|18.98|211.73|\n||50|6.54|19.50|185.31|\n||100|6.23|25.89|370.32|\n|USC1990|10|3.77|7.96|41.67|\n||30|1.71|2.81|48.99|\n||50|1.69|1.33|12.31|\n||100|1.62|1.17|11.00|\n|SUSY|10|4.05|4.83|34.08|\n||30|2.99|4.51|33.96|\n||50|2.87|3.74|28.71|\n||100|2.41|3.74|42.28|\n|HIGGS|10|1.66|1.14|12.82|\n||30|1.56|1.25|14.94|\n||50|1.41|1.08|22.56|\n||100|1.41|2.04|41.22|\n\n*Table 2: Ratio of gain range to the average cost on small datasets*\n\n|**Dataset**|**$k$**|**Mean**|**Variance**|**Max**|\n|---|---|---|---|---|\n|iris|3|4.98|3.32|13.61|\n||5|3.54|3.16|15.42|\n||10|2.36|1.21|7.02|\n|seeds|3|3.84|2.39|11.40|\n||5|2.99|1.46|8.41|\n||10|2.66|1.45|9.81|\n|glass|3|4.31|3.15|16.21|\n||5|3.65|3.56|14.74|\n||10|3.01|3.12|15.88|\n|Who_FL|3|5.37|4.74|32.04|\n||5|6.19|7.74|37.71|\n||10|4.67|8.42|59.39|\n|HCV_L|3|7.21|13.31|62.19|\n||5|8.34|19.14|96.74|\n||10|3.68|8.22|78.61|\n|TRR_FL|3|1.31|0.34|2.53|\n||5|1.36|0.41|2.75|\n||10|1.24|0.38|2.73|\n|urbanGB_L|3|4.89|3.11|14.68|\n||5|5.47|4.22|22.74|\n||10|5.37|5.16|28.84|\n\n**Weakness 2: As mentioned in the introduction of this paper, there exist coreset-based algorithms achieving running time of $\\tilde{O}(nd+k^3d)$, which is better than the new algorithm theoretically. This weakens the contribution.**\n\nResponse: We thank the reviewer for raising this concern. While it is true that some coreset-based algorithms for $k$-means achieve a theoretical runtime of $\\tilde{O}(nd + k^3d)$ under no distribution assumptions, this may overlook the practical overhead involved in coreset construction—particularly for large-scale datasets. Experimental results indicate that there is a sharp increase on the coreset construction time as data sizes grow (see Tables 5 and 6 in our appendix). Moreover, coreset-based methods often produce worse clustering quality than algorithms that directly optimize the original objective.\n\nIn contrast, our method bypasses the data compression step required by coreset construction and achieves comparable time complexity with much better empirical performance. As demonstrated in our experiments (see our Appendix), it delivers up to 80× speedup over coreset-based methods on large-scale datasets such as SIFT (100M points), while achieving better clustering quality. Therefore, even in the presence of coreset methods with similar time complexity, our approach contributes a more practical and scalable solution with theoretical guarantees.\n\n**Weakness 3 and Question 1: Could the authors provide a comparison between their new algorithm and a faster coreset-based approach over one million data points? For instance, a possible method is to apply uniform sampling to obtain a weak coreset and run -Means++ on the coreset. If additional experimental results could demonstrate that the proposed algorithm still outperforms such coreset-based methods, I think it will increase the contribution.**\n\nResponse: We thank the reviewer for the constructive feedbacks. We agree with the reviewer that an unweighted weak coreset can be constructed efficiently via uniform sampling when the dataset exhibits the balance property. While this approach offers excellent runtime performance, it often results in worse clustering qualities. As it targets a different trade-off between accuracy and efficiency, we did not include it in our original comparison.\n\nTo address this concern, we conducted additional experiments on large-scale datasets such as USC1990 (1M), HIGGS (11M), and SIFT (100M), where we apply uniform sampling (with sizes ranging from 1,000 to 1,000,000) to construct the coreset followed by $k$-means++ to obtain the clustering solutions. The results (Tables 3-5) show that even with 1 million samples, our method consistently achieves much better clustering quality than uniform coresets while maintaining comparable runtime. On average, our method reduces the clustering cost by 5% compared to the uniform coreset baseline. \n\nWe will include these results in the revised version to provide a more complete comparison.\n\n*Table 3: Comparisons with weak coreset on USC1990*\n\n|**Dataset**|**k**|**Method**|**Cost**|**Time(s)**|\n|---|---|---|---|---|\n|USC1990|10|Uniform+KM(size=1000)|2.68E+08$\\pm$1.1E+07|0.17|\n|||Uniform+KM(size=10000)|2.57E+08$\\pm$1.2E+07|**0.08**|\n|||Uniform+KM(size=100000)|2.59E+08$\\pm$1.2E+07|0.31|\n|||Uniform+KM(size=500000)|2.58E+08$\\pm$1.3E+07|1.29|\n|||Uniform+KM(size=1000000)|2.54E+08$\\pm$4.9E+06|7.52|\n|||Ours|**2.40E+08$\\pm$4.1E+06**|1.12|\n||50|Uniform+KM(size=1000)|1.46E+08$\\pm$1.8E+06|**0.08**|\n|||Uniform+KM(size=10000)|1.34E+08$\\pm$1.3E+06|0.13|\n|||Uniform+KM(size=100000)|1.33E+08$\\pm$6.0E+05|0.34|\n|||Uniform+KM(size=500000)|1.34E+08$\\pm$4.4E+05|1.96|\n|||Uniform+KM(size=1000000)|1.36E+08$\\pm$4.2E+05|11.76|\n|||Ours|**1.29E+08$\\pm$1.6E+06**|2.10|\n||100|Uniform+KM(size=1000)|1.17E+08$\\pm$9.8E+05|**0.06**||\n|||Uniform+KM(size=10000)|1.05E+08$\\pm$4.8E+05|0.25|\n|||Uniform+KM(size=100000)|1.04E+08$\\pm$1.3E+05|0.52|\n|||Uniform+KM(size=500000)|1.03E+08$\\pm$1.6E+05|3.25|\n|||Uniform+KM(size=1000000)|1.03E+08$\\pm$5.0E+04|18.14|\n|||Ours|**9.99E+07$\\pm$1.5E+06**|3.58|\n\n*Table 4: Comparisons with weak coreset on HIGGS*\n\n|**Dataset**|**k**|**Method**|**Cost**|**Time**|\n|---|---|---|---|---|\n|HIGGS|10|Uniform+KM(size=1000)|1.67E+08$\\pm$1.2E+06|0.71|\n|||Uniform+KM(size=10000)|1.69E+08$\\pm$2.2E+06|**0.26**|\n|||Uniform+KM(size=100000)|1.66E+08$\\pm$4.5E+05|3.28|\n|||Uniform+KM(size=500000)|1.67E+08$\\pm$1.1E+05|7.33|\n|||Uniform+KM(size=1000000)|1.67E+08$\\pm$7.7E+05|12.03|\n|||Ours|**1.56E+08$\\pm$5.4E+05**|5.73|\n||50|Uniform+KM(size=1000)|1.36E+08$\\pm$4.6E+05|**0.77**|\n|||Uniform+KM(size=10000)|1.33E+08$\\pm$3.4E+05|0.79|\n|||Uniform+KM(size=100000)|1.32E+08$\\pm$4.1E+05|8.06|\n|||Uniform+KM(size=500000)|1.32E+08$\\pm$1.5E+05|10.44|\n|||Uniform+KM(size=1000000)|1.32E+08$\\pm$2.6E+05|19.75|\n|||Ours|**1.26E+08$\\pm$3.3E+05**|7.88|\n||100|Uniform+KM(size=1000)|1.27E+08$\\pm$5.3E+05|**0.45**||\n|||Uniform+KM(size=10000)|1.20E+08$\\pm$2.0E+05|0.84|\n|||Uniform+KM(size=100000)|1.20E+08$\\pm$1.8E+05|3.48|\n|||Uniform+KM(size=500000)|1.20E+08$\\pm$2.0E+05|4.09|\n|||Uniform+KM(size=1000000)|1.19e+08$\\pm$2.6E+05|14.68|\n|||Ours|**1.14E+08$\\pm$1.3E+05**|9.87|\n\n*Table 4: Comparisons with weak coreset on SIFT*\n\n|**Dataset**|**k**|**Method**|**Cost**|**Time**|\n|---|---|---|---|---|\n|SIFT|10|Uniform+KM(size=1000)|1.10E+13$\\pm$8.6E+10|1.32|\n|||Uniform+KM(size=10000)|1.10E+13$\\pm$1.7E+10|**1.10**|\n|||Uniform+KM(size=100000)|1.10E+13$\\pm$9.2E+09|8.25|\n|||Uniform+KM(size=500000)|1.10E+13$\\pm$2.0E+09|14.76|\n|||Uniform+KM(size=1000000)|1.09E+13$\\pm$1.3E+08|33.24|\n|||Ours|**1.05E+13$\\pm$3.9E+10**|99.81|\n||50|Uniform+KM(size=1000)|9.51E+12$\\pm$2.7E+10|**0.39**|\n|||Uniform+KM(size=10000)|9.21E+12$\\pm$6.6E+09|0.58|\n|||Uniform+KM(size=100000)|9.14E+12$\\pm$1.4E+10|1.96|\n|||Uniform+KM(size=500000)|9.17E+12$\\pm$7.9E+08|7.50|\n|||Uniform+KM(size=1000000)|9.15E+12$\\pm$1.5E+09|72.05|\n|||Ours|**8.74E+12$\\pm$1.1E+10**|125.89|\n||100|Uniform+KM(size=1000)|9.17E+12$\\pm$5.4E+09|**0.11**|\n|||Uniform+KM(size=10000)|8.58E+12$\\pm$1.2E+10|0.58|\n|||Uniform+KM(size=100000)|8.49E+12$\\pm$1.2E+10|2.11|\n|||Uniform+KM(size=500000)|8.49E+12$\\pm$1.2E+10|7.02|\n|||Uniform+KM(size=1000000)|8.46E+12$\\pm$3.7E+09|96.31|\n|||Ours|**8.13E+12$\\pm$5.0E+09**|168.32|"}, {"review_id": "oRP2DirD7O", "reviewer": "Reviewer_MdzR", "summary": "This paper introduces two fast local-search algorithms for Euclidean $k$-means and $k$-median that avoid the traditional $\\Omega(n d k)$ running time for local search algorithms.\n\n1. BanditLS combines an adaptive $D^2$-sampling surrogate—based on a binary tree that stores distance weights and is updated per swap—with a multi-armed-bandit selection rule that groups swap pairs and prunes low-impact arms. Under a sub-Gaussian prior on swap-pair gains, it achieves a constant-factor approximation for k-means in expected $\\tilde O(nd + \\sqrt{n} k^4 d)$ time.\n1. A second variant replaces the tree walk by a Metropolis–Hastings approximation of D²-sampling, removing the $\\sqrt{n}$ factor under some additional assumptions on optimal clusters, yielding $\\tilde O(nd + k^4 d)$ time.\n1. The same adaptive-sampling idea, plus weighted estimators, is extended to $k$-median, giving expected $\\tilde O(nd + \\sqrt{n} k^3 d)$ runtime without distributional assumptions\n\nLarge-scale experiments show running-time improvements over prior local-search methods and coreset-based techniques while outputting high-quality cluster centers.", "questions": "1. How reasonable are the sub-Gaussian assumptions for Theorems 1 and 2?\n1. Does the algorithm scale to high-dimensional datasets? The experiments seems to vary $n, k$ but not $d$.", "limitations": "yes", "rating": 5, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 3, "significance": 3, "originality": 3, "strengths_and_weaknesses": "### Strengths\n1. To the best of my knowledge, this is the first work that breaks the $O(ndk)$ running-time barrier for local-search-based $k$-means and $k$-medians clustering.\n1. Empirical results are pretty convincing: the running time beats $k$-means++ while outputting better quality clusters on reasonably large datasets.\n1. The dynamic tree-based importance sampling and best-arm identification ideas may extend to other tasks.\n\n\n---\n\n### Weaknesses\n1. The theoretical approximation ratios are quite large (2000+)\n2. The running time guarantees depend on some non-standard assumptions, and it is unclear how reasonable these assumptions are.", "was_revised": true, "final_justification": "The rebuttal more than addressed my concerns. I believe this is a nice contribution and maintain my current score.", "rebuttal": "We thank the reviewer for the positive rating and thoughtful comments. Below, we address the concerns.\n\n**Question 1: How reasonable are the sub-Gaussian assumptions for Theorems 1 and 2?**\n\nResponse: We thank the reviewer for this important question. For the $k$-means objective, it is based on squared distances. This makes our assumption-free weighted sampling strategy (originally designed for $k$-median) inapplicable, as it can lead to significant approximation errors in estimating swap gains. To address this issue, we formulate the best swap pair selection as a bandit problem. In the bandit setting, the sub-Gaussian assumption on arm gains is a standard and widely adopted condition for establishing theoretical guarantees (e.g., [1], [2] and [6], [30], [31] in our ref list). Thus, the sub-Gaussian assumption used in Theorems 1 and 2 is both reasonable and well supported by prior work.\n\nMoreover, previous assumptions for bandit-based methods are much stronger than ours. In [6], [30], and [31] (see our ref list), it is assumed that each arm’s gain follows a fixed distribution with mean $\\gamma$ and variance 1, or that the variance-to-maximum-gain ratio is bounded by a constant. In contrast, we only require the variance to be bounded by a function of the average clustering cost, making our assumption much weaker than previous ones.\n\nTo empirically verify the sub-Gaussian assumption, we use the fact that any variable with range $2\\sigma$ is $\\sigma$-sub-Gaussian. Thus, we assess sub-Gaussianity by computing the ratio between the gain range and the average clustering cost $\\delta^2(P, C)/n$, where the gain of a point $p$ under swap pair $(u, v)$ and a center set $C$ is defined as $\\delta^2(p, C \\setminus u \\cup v) - \\delta^2(p, C)$. We run our algorithm 5 times with 400 local search steps and calculate the ratios. In Tables 1 and 2, we give the mean, variance, and maximum of these ratios. The results show that the ratio remains below 50 on most datasets except for dataset SYN and HCV_L. This indicates that the sub-Gaussian assumption is both reasonable and empirically justified in practice.\n\n*Table 1: Ratio of gain range to average cost on large datasets*\n\n|**Dataset**|**$k$**|**Mean**|**Variance**|**Max**|\n|---|---|---|---|---|\n|SYN|10|9.75|11.89|77.92|\n||30|7.97|18.98|211.73|\n||50|6.54|19.50|185.31|\n||100|6.23|25.89|370.32|\n|USC1990|10|3.77|7.96|41.67|\n||30|1.71|2.81|48.99|\n||50|1.69|1.33|12.31|\n||100|1.62|1.17|11.00|\n|SUSY|10|4.05|4.83|34.08|\n||30|2.99|4.51|33.96|\n||50|2.87|3.74|28.71|\n||100|2.41|3.74|42.28|\n|HIGGS|10|1.66|1.14|12.82|\n||30|1.56|1.25|14.94|\n||50|1.41|1.08|22.56|\n||100|1.41|2.04|41.22|\n\n*Table 2: Ratio of gain range to average cost on small datasets*\n\n|**Dataset**|**$k$**|**Mean**|**Variance**|**Max**|\n|---|---|---|---|---|\n|iris|3|4.98|3.32|13.61|\n||5|3.54|3.16|15.42|\n||10|2.36|1.21|7.02|\n|seeds|3|3.84|2.39|11.40|\n||5|2.99|1.46|8.41|\n||10|2.66|1.45|9.81|\n|glass|3|4.31|3.15|16.21|\n||5|3.65|3.56|14.74|\n||10|3.01|3.12|15.88|\n|Who_FL|3|5.37|4.74|32.04|\n||5|6.19|7.74|37.71|\n||10|4.67|8.42|59.39|\n|HCV_L|3|7.21|13.31|62.19|\n||5|8.34|19.14|96.74|\n||10|3.68|8.22|78.61|\n|TRR_FL|3|1.31|0.34|2.53|\n||5|1.36|0.41|2.75|\n||10|1.24|0.38|2.73|\n|urbanGB_L|3|4.89|3.11|14.68|\n||5|5.47|4.22|22.74|\n||10|5.37|5.16|28.84|\n\n[1] Baharav T Z, Kang R, Sullivan C, et al. Adaptive sampling for efficient softmax approximation[C]//Proceedings of the 38th International Conference on Neural Information Processing Systems. 2024: 117580-117613.\n\n[2] Tiwari M, Kang R, Lee J Y, et al. Faster maximum inner product search in high dimensions[C]//Proceedings of the 41st International Conference on Machine Learning. 2024: 48344-48361.\n\n**Question 2: Does the algorithm scale to high-dimensional datasets? The experiments seems to vary n, k but not d.**\n\nResponse: We thank the reviewer for this insightful question. Our algorithm consists of two main components: (1) sampling via the Metropolis-Hastings process, and (2) bandit-based swap pair selection. In the Metropolis-Hastings step, samples are drawn uniformly, and only their distances to the current centers are computed, requiring time linear in the dimension $d$. Similarly, the bandit-based selection evaluates swap gains on a small sample of size $\\tilde{O}(k^2)$, where each step also involves distance computations with runtime linear in $d$. Therefore, the overall runtime scales linearly with $d$, making the algorithm suitable for high-dimensional datasets.\n\nTo empirically verify this, we conducted additional experiments on several high-dimensional datasets from diverse domains. These include: (1) the Twin Gas Sensor Arrays (TGas) dataset with 40,000 points and 2,000 dimensions; (2) the CIFAR-10 image dataset represented as 10,000 flattened 3,072-dimensional vectors; and (3) the Electricity Load Diagrams (ELD) dataset with 370 instances and 140,256 dimensions. These datasets were selected to cover a wide range of $d$ values and application domains. \n\nTables 3-5 give the comparison results on these datasets. It can be seen that our algorithm consistently scales well in all cases, achieving an average 10x speedups over baseline local search methods. Moreover, it maintains comparable or even better clustering quality despite the increase in dimensionality. These findings demonstrate that our method is robust not only to data size $n$ and number of clusters $k$, but also to the dimensionality $d$.\n\nWe will include these additional experimental results and corresponding runtime analysis in the appendix of the revised manuscript.\n\n*Table 3: Comparison results on dataset TGas*\n\n|**Dataset**|**k**|**Method**|**Cost**|**Time(s)**|\n|---|---|---|---|---|\n|TGas(400000, 2000)|10|k-means++|2.69E+05$\\pm$5.0E+02|**26.77**|\n|||MLS|-|>6h|\n|||LSDS++|2.68E+05$\\pm$1.1E+03|984.67|\n|||LS++|2.63E+05$\\pm$3.1E+02|482.69|\n|||Ours|**2.62E+05$\\pm$1.2E+03**|31.13|\n||50|k-means++|2.61E+05$\\pm$5.9E+02|111.96|\n|||MLS|-|>6h|\n|||LSDS++|2.60E+05$\\pm$1.1E+03|1047.36|\n|||LS++|**2.57E+05$\\pm$1.6E+02**|624.44|\n|||Ours|2.58E+05$\\pm$1.4E+02|**32.53**|\n||100|k-means++|2.58E+05$\\pm$4.0E+02|228.96|\n|||MLS|-|>6h|\n|||LSDS++|2.58E+05$\\pm$1.0E+03|1054.94|\n|||LS++|**2.57E+05$\\pm$2.4E+02**|881.86|\n|||Ours|2.57E+05$\\pm$3.5E+02|**42.04**|\n\n*Table 3: Comparison results on dataset CIFAR10*\n\n|**Dataset**|**k**|**Method**|**Cost**|**Time(s)**|\n|---|---|---|---|---|\n|CIFAR10(10000, 3072)|10|k-means++|8.97e+10$\\pm$8.2E+09|**0.43**|\n|||MLS|8.58E+10$\\pm$1.5E+09|35.34|\n|||LSDS++|8.86E+10$\\pm$1.1E+09|38.93|\n|||LS++|8.63E+10$\\pm$4.1E+08|5.52|\n|||Ours|**8.58E+10$\\pm$2.0E+08**|2.50|\n||50|k-means++|7.51E+10$\\pm$2.4E+08|**1.18**|\n|||MLS|7.28E+10$\\pm$6.1E+08|375.17|\n|||LSDS++|7.47E+10$\\pm$9.7E+08|43.45|\n|||LS++|7.28E+10$\\pm$2.5E+08|6.64|\n|||Ours|**7.28E+10$\\pm$2.3E+08**|6.71|\n||100|k-means++|7.38E+10$\\pm$2.5E+08|**2.19**|\n|||MLS|7.00E+10$\\pm$3.7E+08|2855.34|\n|||LSDS++|7.14E+10$\\pm$5.7E+08|47.68|\n|||LS++|6.93E+10$\\pm$2.2E+08|7.24|\n|||Ours|**6.79E+10$\\pm$2.3E+08**|12.24|\n\n*Table 3: Comparison results on dataset ELD*\n\n|**Dataset**|**k**|**Method**|**Cost**|**Time**|\n|---|---|---|---|---|\n|ELD(370, 140256)|10|k-means++|3.62E+06$\\pm$8.6E+04|**0.20**|\n|||MLS|3.43E+06$\\pm$1.6E+04|11.00|\n|||LSDS++|3.53E+06$\\pm$7.6E+05|9.33|\n|||LS++|**3.43E+06$\\pm$6.5E+03**|16.69|\n|||Ours|3.43E+06$\\pm$7.3E+03|0.59|\n||50|k-means++|2.12E+06$\\pm$2.2E+04|**0.58**|\n|||MLS|2.04E+06$\\pm$5.9E+03|76.27|\n|||LSDS++|2.09E+06$\\pm$1.9E+05|14.00|\n|||LS++|**1.98E+06$\\pm$2.8E+03**|22.77|\n|||Ours|2.02E+06$\\pm$4.7E+03|1.06|\n||100|k-means++|1.69E+06$\\pm$3.8E+03|**1.09**|\n|||MLS|1.67E+06$\\pm$1.8E+03|419.52|\n|||LSDS++|1.66E+06$\\pm$8.6E+03|22.13|\n|||LS++|**1.61E+06$\\pm$2.7E+03**|30.86|\n|||Ours|1.65E+06$\\pm$9.5E+03|1.60|"}, {"review_id": "lc2n5wAuSh", "reviewer": "Reviewer_BsVn", "summary": "The paper designs faster algorithms for k-means in $R^d$, tackling the challenge of achieving running time better than $O(ndk)$, which is obviously important when all three parameters $n,d,k$ are large. In essence, the goal is to decouple the product $nk$, and obtain running time like $O(nd + d k^2)$. The paper designs two algorithms, that achieve O(1) approximation for k-means in time $\\tilde{O}(nd + \\sqrt{n} k^4 d)$ under some assumption, and better time $\\tilde{O}(nd + k^4 d)$ under more assumptions. The second result can be extended to k-median without assumptions. \n\nThe said assumptions about the data are somewhat indirect, and roughly assert that for the algorithm can evaluate a potential swap operation on the current set of centers using a random point from the data. The general approach is to speed up the well-known local search method (start with random centers and repeatedly swap one center), using two ideas. One idea is to represent the distances approximately using a randomized tree embedding based on [12], because the tree structure can quickly sample proportional to the tree distances, and also quickly processing an update that changes one of the centers. A second idea is to model a swap of one center as a bandit problem, and estimate the benefit of each swap using a sample of $\\tilde{O}(k^2)$ points, and here the assumption is crucial. \n\nThe empirical evaluation compares the new algorithm to 4 previous algorithms, with real-life datasets of different sizes, and with k ranging between 10 and 100. he new algorithm is faster than all other algorithms, especially as k gets larger, and the accuracy is comparable to the best among the 4 previous ones.", "questions": "Your work is important when k is large, but in that case why is it important to have exactly k clusters, and not employ bicriteria approximation, which allows $O(k)$ clusters, or even uniform facility location where the number of clusters is controlled indirectly via a penalty term? \n\nIn Lemma 1, can you explain how this bound on the expectation of the tree distance $\\delta_T$, but not on the expectation of its square, is useful for k-means, which corresponds to squared distances? In addition, why do you use [12] and not the well-known tree embedding that corresponds to a randomly shifted quadtree, and has stretch bound $O(d\\log\\Delta)$ instead of $O(d)$ in Lemma 1? \n\nIn Lemma 4, what is the intuition for the term $\\sqrt n$, I mean where it is coming from? In addition, this lemma shows that each point x is sampled almost with the right probability, up to factor $\\sqrt n$; can we use $O(\\sqrt n)$ such samples and rejection sampling to obtain one sample with a precisely correct probability? Doesn't this require computing many distances in $R^d$ (not in the tree) which would take too much time? \n\nIn assumption 2, what does it mean that \"P is average\" ?  In addition, this assumption sounds very strong, in the sense of making the problem much easier. For example, sampling the points uniformly at rate $(\\log n)/k^2$ will reduce the input size by factor almost $k^2$ and the optimal clusters/centers would probably (I have not verified it) change only a little bit. Why do you think this assumption is reasonable? Has it been used before?", "limitations": "yes", "rating": 4, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": 3, "clarity": 2, "significance": 4, "originality": 3, "strengths_and_weaknesses": "The main strength in my opinion is tackling the challenge of breaking below running time $O(ndk)$, particularly for k-means where a randomized tree embedding is usually not effective (it does not provide good guarantees for squared distances). Another advantage is that the new algorithm is indeed faster, both in theory and in practice, while still achieving high accuracy. The approach seems new, even though it is based on known components, like a randomized tree embedding, and formulating swap pair selection as a bandit problem. \n\nThe main weakness of the theoretical results is that the assumptions are rather indirect and not intuitive, and in fact refer to the algorithm itself (e.g. the current centers and the proposed swap). \n\nA possible weakness of the empirical evaluation is that it compares the different algorithms' performance on data that essentially does not cluster well. Table 2 shows two datasets, and in both of them, when the number of clusters k doubles (from 20-30 to 50 and then to 100), the clustering cost goes down by about 10%. I suspect that the good accuracy reported here simply misses the intended application, although I am not sure what is a good or standard way to control for this gap.  \n\nBottom line, I really like the direction and the results sound promising, however I have to be skeptic because the theoretical assumptions seem strong and some pieces of the proof do not fit, as mentioned in my questions for the authors.", "was_revised": true, "final_justification": "The rebuttal does not change my overall evaluation, but does clarify some technical issues.", "rebuttal": "We thank the reviewer for the positive rating and thoughtful comments. Below, we address the concerns.\n\n**Question 1: Your work is important when k is large, but in that case why is it important to have exactly k clusters, and not employ bicriteria approximation, which allows $O(k)$ clusters, or even uniform facility location where the number of clusters is controlled indirectly via a penalty term?**\n\nResponse: We thank the reviewer for raising this insightful question. In the standard $k$-means objective, enforcing exactly $k$ centers is crucial for fair comparisons against prior methods. By fixing $k$, it allows for more meaningful comparisons of theoretical results and clearer assessments of algorithmic effectiveness.\n\nWe agree with the reviewer that allowing more than $k$ centers to be opened can reduce the clustering cost and improve approximation ratios. Bicriteria approximation is therefore a practical and effective alternative in real-world scenarios. Extending our  framework to obtain bicriteria solutions and achieve better practical performance is an interesting direction for further studies.\n\nAs for the uniform facility location problem, to the best of our knowledge, existing local search algorithms do not achieve linear running time. This is mainly because the opening costs may affect the sampling process and complicate the design of sampling distributions. Developing  local search methods with linear runtime for facility location remains an interesting open problem for future work.\n\n**Question 2: In Lemma 1, can you explain how this bound on the expectation of the tree distance $\\delta_T$, but not on the expectation of its square, is useful for k-means, which corresponds to squared distances? In addition, why do you use [12] and not the well-known tree embedding that corresponds to a randomly shifted quadtree, and has stretch bound $O(d\\log\\Delta)$ instead of $O(d)$ in Lemma 1?**\n\nResponse: We appreciate the reviewer’s insightful question. \n\nFor the $k$-means problem, our rejection sampling relies on a tree embedding to approximate the squared-distance-based sampling distribution. A key challenge is that directly bounding the expectation of squared tree distances through tree metric leads to large distortion. Prior work (e.g., [1]) shows that the expected squared distance distortion under tree embeddings can be as large as $\\Omega(n)$, which adds an $O(n)$ factor loss in sampling probabilities and runtime.\n\nTo address this, we instead bound the expected tree distance under the $\\ell_1$ objective, as stated in Lemma 1. Then, by carefully relating $\\ell_1$ objective to squared distances through the Cauchy–Schwarz Inequality, we reduce the distortion in sampling probabilities from $O(n)$ to $\\tilde{O}(\\sqrt{n})$.\n\nRegarding the quadtree embedding, the method in [12] is indeed based on the well-known quadtree embedding technique for $\\ell_1$ distances, which has a distortion of $O(d \\log \\Delta)$. We apologize for the typo in Lemma 1, where the distortion should be $\\tilde{O}(d)$ instead of $O(d)$ to reflect the hidden logarithmic dependence.\n\nTo ensure that the $\\log \\Delta$ factor does not affect the overall runtime, we include a preprocessing step (lines 198–199 in our paper) that rescales the dataset to reduce the aspect ratio $\\Delta$ to $\\mathrm{poly}(n, d)$. This step runs in $\\tilde{O}(nd)$ time and guarantees that $\\log \\Delta = \\tilde{O}(1)$. As a result, the distortion becomes $\\tilde{O}(d)$, which only contributes an $\\tilde{O}(1)$ overhead and does not change the overall time complexity of our algorithm. Thus, the quadtree embedding remains compatible with our rejection sampling framework.\n\nWe will clarify this in the revised version.\n\n[1] Cohen-Addad et al. Streaming euclidean k-median and k-means with o(log n) space. Proc. FOCS 2023, pp. 883-908.\n\n**Question 3：In Lemma 4, what is the intuition for the term $O(\\sqrt{n})$, I mean where it is coming from? In addition, this lemma shows that each point $x$ is sampled almost with the right probability, up to factor $O(\\sqrt{n})$; can we use such samples and rejection sampling to obtain one sample with a precisely correct probability? Doesn't this require computing many distances in $R^d$ (not in the tree) which would take too much time?**\n\nResponse: We thank the reviewer for the thoughtful question. \n\nThe $O(\\sqrt{n})$ factor loss arises from approximating squared distances using $\\ell_1$ objective in the sampling probability. For our rejection sampling process, to avoid the $O(n)$ factor loss caused by squared distance embedding, the $\\ell_1$ distances are used. However, we still need to relate the $\\ell_1$-based sampling distribution to the one suitable for squared distances. This introduces an $O(\\sqrt{n})$ factor loss through Cauchy–Schwarz inequality.\n\nFor the sampling process, the samples drawn from rejection sampling cannot recover the precise sampling probabilities. However, our goal is to ensure that we can sample a “good” point for local search swaps with constant probability. For these points, it is observed that the clustering cost of them account for a certain fraction of the total clustering cost (as shown in Lemma B.3). Since our sampling probabilities are within an $O(1/\\sqrt{n})$ factor, by repeating the rejection sampling process for $\\tilde{O}(\\sqrt{n})$ rounds, it suffices to select a good point with constant probability.\n\nAs for computation cost, each round of rejection sampling only requires computing the true squared distance from a single candidate point to its nearest center in $\\mathbb{R}^d$. Thus, the total cost over $\\tilde{O}(\\sqrt{n})$ rounds is $\\tilde{O}(\\sqrt{n}kd)$, which does not significantly impact the overall runtime.\n\n**Question 4: In assumption 2, what does it mean that $P$ is average\"? In addition, this assumption sounds very strong, in the sense of making the problem much easier. For example, sampling the points uniformly at rate $\\log(n)/k^2$ will reduce the input size by factor almost $k^2$ and the optimal clusters/centers would probably (I have not verified it) change only a little bit. Why do you think this assumption is reasonable? Has it been used before?**\n\nResponse: We thank the reviewer for this thoughtful question. In Assumption 2, $P$ is “average” means that each optimal cluster has size $\\Omega(k^2)$, which implies the clustering instance is not heavily imbalanced. \n\nThis type of assumption is commonly used in the clustering literature. Prior works on clustering algorithm design (e.g., [2]–[4]) often assume that each optimal cluster has size $\\Omega(n/k)$ to ensure meaningful sampling guarantees. Our assumption, which requires each cluster to have size $\\Omega(k^2)$, is consistent with these when $k^3 = O(n)$. This is a reasonable setting in practice, as it typically holds that $n$ is much larger than $k$. \n\nRegarding the reviewer’s observation about uniform sampling at rate $\\log(n)/k^2$, we agree with the reviewer that this helps to reduce the input size to $\\tilde{O}(n/k^2)$ while guaranteeing the theoretical loss (e.g., [2]). However, this approach often leads to worse clustering quality in practice. \n\nTo evaluate this strategy, we conducted additional experiments comparing our method with a baseline that uniformly samples $\\beta$ points (ranging from 1,000 to 1,000,000) and applies a fast local search algorithm (LS++) to the subsample. Tables 1–3 give the comparison on large-scale datasets including SUSY (5M) and SIFT (100M). The results show that while using small samples (e.g., 1,000–10,000) yields faster runtimes than ours, it leads to at least a 5% increase in clustering cost. To match the clustering quality of ours, it requires to sample over 500,000 points. in which case our algorithm runs much faster than the uniform sampling baseline.\n\n[2] Huang et al. The power of uniform sampling for k-median. Proc. ICML 2023, pp. 13933-13956.\n\n[3] Meyerson et al. A k-median algorithm with running time independent of data size. Machine Learning, 2004, 56(1): 61-87.\n\n[4] Huang et al. Bi-criteria Sublinear Time Algorithms for Clustering with Outliers in High Dimensions. Proc. COCOON 2024, pp. 91-103. \n\n*Table 1: Comparisons with uniform sampling on SUSY*\n\n|**Dataset**|**k**|**Method**|**Cost**|**Time(s)**|\n|---|---|---|---|---|\n|SUSY|10|size=1000|3.13E+07$\\pm$4.3E+05|0.49|\n|||size=10000|3.05E+07$\\pm$1.6E+05|**0.43**|\n|||size=100000|3.05E+07$\\pm$2.1E+05|9.32|\n|||size=500000|3.02E+07$\\pm$1.1E+05|47.03|\n|||size=1000000|3.04E+07$\\pm$3.0E+05|558.41|\n|||Ours|**3.01E+07$\\pm$2.3E+05**|3.67|\n||50|size=1000|2.18E+07$\\pm$1.3E+05|0.89|\n|||size=10000|2.01E+07$\\pm$1.4E+05|**0.83**|\n|||size=100000|1.96E+07$\\pm$4.1E+04|20.64|\n|||size=500000|1.96E+07$\\pm$1.5E+05|83.91|\n|||size=1000000|1.96E+07$\\pm$3.1E+05|1052.03|\n|||Ours|**1.95E+07$\\pm$2.2E+05**|4.22|\n||100|size=1000|1.89E+07$\\pm$1.2E+05|0.94|\n|||size=10000|1.70E+07$\\pm$6.7E+04|**0.81**|\n|||size=100000|1.65E+07$\\pm$3.7E+04|20.47|\n|||size=500000|1.64E+07$\\pm$4.8E+04|99.38|\n|||size=1000000|1.64E+07$\\pm$2.5E+05|1647.13|\n|||Ours|**1.63E+07$\\pm$2.9E+04**|6.54|\n\n*Table 2: Comparisons with uniform sampling on SIFT*\n\n|**Dataset**|**k**|**Method**|**Cost**|**Time**|\n|---|---|---|---|---|\n|SIFT|10|size=1000|1.08E+13$\\pm$6.8E+10|3.46|\n|||size=10000|1.07E+13$\\pm$7.6E+10|**2.89**|\n|||size=100000|1.06E+13$\\pm$2.1E+10|31.72|\n|||size=500000|1.06E+13$\\pm$1.9E+10|139.13|\n|||size=1000000|1.06E+13$\\pm$7.5E+09|2814.43|\n|||Ours|**1.05E+13$\\pm$3.9E+10**|99.81|\n||50|size=1000|9.41E+12$\\pm$1.4E+10|**2.74**|\n|||size=10000|8.90E+12$\\pm$4.7E+09|3.69|\n|||size=100000|8.89E+12$\\pm$4.8E+10|42.26|\n|||size=500000|8.85E+12$\\pm$1.3E+10|206.24|\n|||size=1000000|8.80E+12$\\pm$1.5E+10|4693.98|\n|||Ours|**8.74E+12$\\pm$1.1E+10**|125.89|\n||100|size=1000|9.17E+12$\\pm$1.6E+10|**2.73**|\n|||size=10000|8.35E+12$\\pm$2.1E+10|3.59|\n|||size=100000|8.21E+12$\\pm$8.0E+09|44.71|\n|||size=500000|8.20E+12$\\pm$8.1E+09|232.58|\n|||size=1000000|8.15E+12$\\pm$1.9E+09|5537.51|\n|||Ours|**8.13E+12$\\pm$5.0E+09**|168.32|"}]

Author rebuttal:

Meta-review: The overall sentiment from reviewers is positive. In particular, they appreciated the additional experiments provided in the authors’ response, which addressed most concerns regarding the reasonableness of the assumptions and the comparison with coreset methods. Therefore, I believe this paper is suitable for acceptance.
Decision: Accept (poster)