BanditPAM++: Faster $k$-medoids Clustering
Authors: Mo Tiwari
Ryan Kang
Donghyun Lee
Sebastian Thrun
Ilan Shomorony
Martin Jinye Zhang
Venue: NeurIPS 2023 poster (2023)
Keywords: multi-armed bandits
clustering
k-medoids
best-arm identification
TL;DR: We propose a highly efficient k-medoids clustering algorithm with $O(n \log n)$ complexity per iteration and up to 10x faster than state-of-the-art BanditPAM algorithm.

Abstract: Clustering is a fundamental task in data science with wide-ranging applications. In $k$-medoids clustering, cluster centers must be actual datapoints and arbitrary distance metrics may be used; these features allow for greater interpretability of the cluster centers and the clustering of exotic objects in $k$-medoids clustering, respectively. $k$-medoids clustering has recently grown in popularity due to the discovery of more efficient $k$-medoids algorithms. In particular, recent research has proposed BanditPAM, a randomized $k$-medoids algorithm with state-of-the-art complexity and clustering accuracy. In this paper, we present BanditPAM++, which accelerates BanditPAM via two algorithmic improvements, and is $O(k)$ faster than BanditPAM in complexity and substantially faster than BanditPAM in wall-clock runtime. First, we demonstrate that BanditPAM has a special structure that allows the reuse of clustering information $\textit{within}$ each iteration. Second, we demonstrate that BanditPAM has additional structure that permits the reuse of information $\textit{across}$ different iterations. These observations inspire our proposed algorithm, BanditPAM++, which returns the same clustering solutions as BanditPAM but often several times faster. For example, on the CIFAR10 dataset, BanditPAM++ returns the same results as BanditPAM but runs over 10$\times$ faster. Finally, we provide a high-performance C++ implementation of BanditPAM++, callable from Python and R, that may be of interest to practitioners at https://github.com/motiwari/BanditPAM. Auxiliary code to reproduce all of our experiments via a one-line script is available at https://github.com/ThrunGroup/BanditPAM_plusplus_experiments.

BanditPAM++: Faster k-medoids Clustering
Mo Tiwari Stanford University motivari@stanforr edu
Ryan Kang* Stanford University
Donghyun Lee*
University College London
txryank@stanford.co edu
donghyun. lee. 2 21@ucllc....mmm ac. uk
Sebastian Thrun
Chris Piech Stanford University piechêcs siech@cs.stanforee stanford...rror
Stanford University
thrundstanford edu
edu
Ilan Shomoronyt
Martin Jinye Zhangi Camegii Mcllon University martinzh@andrew cmu, edu
University of Illinois at Stbana-Champaign ilams@illinois edu
Abstract
Clustering a fundamental task in data science with wide-ranging applications. In k-medoids clustering, cluster centers must be actual datapoints and arbitrary distance metrics may be used; these features allow for greater interpretability the cluster centers and the clustering of exotic objects in k-medoids clustering, respectively. k-medoids clustering has recently grown in popularity due to the discovery of more efficient k-medoids algorithms. particular, recent research has proposed BanditPAM, a randomized k-medoids algorithm with state-of-the-art complexity and clustering accuracy. In this paper, we present BanditPAMM++ which accelerates BanditPAM via two algorithmic improvements, and is ((k) faster than BanditPAM in complexity and substantially faster than BanditPAM in wall-clock runtime. First, we demonstrate that BanditPAM has special structure that allows the reuse of clustering information within each iteration. Second, we demonstrate that BanditPAM has additional structure that permits the reuse information across different iterations. These observations inspire our proposed algorithm, BandidPAM+++ which returns the same clustering solutions as BanditPAM but often several times faster. For example, on the CIFARI0 dataset, BanditPAM++ returns the same results as BanditPAM but runs over 10x faster. Finally, we provide a high high-performance C++ implemen- tation of BanditPAM++. callable from Python and R, that may be of interest to practitioners at https:/ http///gttttt com/ mmottiarri/Baddddae code to reproduce all of experiments via one-line script is available https: //gi:t higg.th.b....... com www.hhhnnrru/Baaniiiiiroa PhrrunGroop/Bmmiii..A.Pllpggug plusplus, experiments/
Auxiliary
1 Introduction
Clustering is critical and ubiquitous task in data science and machine learning. Clustering aims to separate dataset of n datapoints into & disjoint sets that form partition of the original dataset. Intuitively, datapoints within cluster are similar and datapoints in different clusters are dissimilar. Clustering problems and algorithms have found numerous applications textual data 200 social
network analysis [18] biology (28). and education (28).
37th Conference on Neural Information Processing Systems (NcurIPS 2023).

---

A common objective function used in clustering is Equation |
L(C) = Σ min d(c, zy).
aEC
Under this loss function, the goal becomes to minimize the distance, defined by the distance function d, from each datapoint its nearest cluster center e among the set of cluster centers Note that this formulation is general and does not require the datapoints to be vectors or assume specific functional form of the distance function d.
Specific choices of C, dataset X' and distance function give rise to different clustering problems. Perhaps of the most commonly used clustering methods is k-means clustering [7 [6]. In k-means clustering, each datapoint is a vector in IR and the distance function d is usually taken to be squared L2 distance; there are no constraints on other than that it subset of R². Mathematicallyy the k-means objective function is
L(C) = Σ j=1
min J cl12 aEC
subject to |C| k.
The most common algorithm for the k-means problem Lloyd iteration [[6], which has been improved by other algorithms such as k-means++ (2). These algorithms are widely used in practice due their simplicity and computational efficiency. Despite widespread use in practice, k-means clustering suffers from several limitations. Perhaps its most significant restriction is the choice of d as the squared L2 distance. This choice of d for computational efficiency - the mean of many points can be efficiently computed under squared L2 distance but prevents clustering with other distance metrics that may be preferred in other contexts [22 g For example, k-means difficult to use on textual data that necessitates string edit distance 20 social network analyses using graph metrics [8]. or sparse datasets (such as those found in recommendation systems [5]) that lend themselves to other distance functions. While k-means algorithms have been adapted to specific other metrics, e.g., cosine distance [23], these methods are bespoke to the metric and not readily generalizable. Another limitation of k-means clustering is that the set of cluster centers C may often be uninterpretable, as each cluster center is (generally) a linear combination datapoints. This limitation can be especially problematic when dealing with structured data, such parse trees in context-free grammars, where the mean of trees not necessarily well-defined, or images in computer vision, where the mean image appear as random noise 28 [5].
requires the cluster centers to In contrast with k-mmaan clustering, k-medoids clustering I0 I1
be actual datapoints, i.e., requires C X. More formally, the objective is to find a set of medoids M X (versus TRP in k-means) that minimizes
L(M)
Σ min d(m, 27,)
MEM
subject to [M] = A: Note that there is no restriction on the distance function d.
k-medoids clustering has several advantages over k-means. Crucially, the requirement that each cluster center is a datapoint leads to greater interpretability of the cluster centers because each cluster center can be inspected. Furthermore, k-medoids supports arbitrary dissimilarity measures; the Equation (3) need not be a proper metric, i.c., may be negative, asymmetric,
distance function or violate the triangle inequality. Because k-medoids supports arbitrary dissimilarity measures, can also be used to cluster "exotic" objects that are not vectors in R", such as trees and strings (8 without embedding them in RP.
The k-medoids clustering problem in Equation (3) is combinatorial optimization algorithm that is ha-hard in general [25]: as such, algorithms for k-medoids clustering are restricted to heuristic solutions. A popular early heuristic solution for the k-medoids problem was the Partitioning Around Medoids (PAM) algorithm [[: however, PAM is quadratic in dataset size in each iteration, which

---

is prohibitively expensive on large dataset. Improving the computational efficiency of these heuristic solutions is an active area of research. Recently, [28] proposed BanditPAM, the first subquadratic algorithm for the k-medoids problem that matched prior state-of-the-ar solutions in clustering quality. this work, we propose BanditPAM++, which improves the computational efficiency of BanditPAM while maintaining the same results. We anticipate these computational improvements will be important in the era of big data, when k-medoids clustering is used on huge datasets.
Contributions: We propose a new algorithm, BandiiPAM++, that is significantly more computation- ally efficient than PAM and BanditPAM but returns the same clustering results with high probability. BanditPAM++ is (Kk) faster than BanditPAM in complexity and substantially faster than BanditPAM in actual runtime wall-clock runtime. Consequently, BanditPAMM+ is faster than prior state-of-the-as k-medoids algorithms while maintaining the same clustering quality. BanditPAM++ is based on two observations about the structure of BanditPAM and the k-medoids problem, described in Section 4 The first observation leads to technique that we call Virtual Arms (VA). The second observation leads to technique that we refer to as Permutation Invariant Caching (PIC). We combine these techniques in BandidPAM++ and prove (in Section (5) and experimentally validate (in Section 6 that BanditPAM++ returns the same solution to the k-medoids problem as PAM and BanditPAM with high probability. but is more computationally efficient. In some instances, BanditPAM++ is over 10x faster than BanditPAM.. Additionally, we provide a highly optimized implementation of BanditPAM++ in C++ that is callable from Python and R and may be of interest practitioners.
2 Related Work
As discussed in Section [ global optimization of the k-medoids problem (Equation (3) is NP-hard in general [25). Recent work attempts to perform attempts global optimization and able to achieve an optimality gap of 0.1 % on one million datapoints, but is restricted to L2 distance and takes several hours run on commodity hardware (24
Because of the difficulty of global optimization of the k-medoids problem, many heuristic algorithms have been developed for the k-medoids problem that scale polynomially with the dataset size and number of clusters. The complexity of these algorithms is measured by their sample complexity. i.e., the number of pairwise distance computations that are computed; these computations have been observed to dominate runtime costs and, as such, sample complexity translates to wall-clock time via an approximately constant factor 28 (this is also consistent with our experiments in Section and Appendix
Among the heuristic solutions for the k-medoids problem, the algorithm with the best clustering loss Partitioning Around Medoids (PAM) I0 IT. which consists of two phases: the BUILD phase and the SWAP phase. However, the BUILD phase and each SWAP iteration of PAM perform O(kn²) distance computations, which can be impractical for large datasets or when distance computations expensive. We provide greater details about the PAM algorithm in Section 3 because it is an important baseline against which we assess the clustering quality of new algorithms.
Though PAM achieves the best clustering loss among heuristic algorithms, the era of huge data has necessitated the development of faster k-medoids algorithms in recent years. These algorithms have typically been divided into two categories: those that agree with PAM and recover the same solution to the k-medoids problem but scale quadratically in n. and those that sacrifice clustering quality for runtime improvements. In the former category, (25] proposed deterministic algorithm called FastPAMI, which maintains the same output as PAM but reduces the computational complexity of each SWAP iteration from (kn)²) to O(n²). However, this algorithm still scales quadratically in in every iteration, which is prohibitively expensive on large datasets.
Faster heuristic algorithms have been proposed but these usually sacrifice clustering quality; such algorithms include CLARA [II, CLARANS 21), and FastPAM (25) While these algorithms scale subquadratically in they return substantially worse solutions than PAM 28). Other algorithms with better sample complexity, such as optimizations for Euclidean space and those based on tabu search heuristics
G also return worse solutions. Finally, D attempts to minimize the number of unique pairwise distances or adaptively estimate these distances or coordinate- wise distances in specific settings [4 3], but all these approaches sacrifice clustering quality for runtime.
3

---

Recently, (28] proposed BanditPAM, state-of-the-and art k- k-medoiss algorithm that arrives at the same solution as PAM with high probability in Okn log n) time. BanditPAM borrows techniques from the multi-armed bandit literature to sample pairwise distance computations rather than compute all (n²)). In this work, we show that BanditPAM can be made more efficient by reusing distance computations both within iterations and across iterations.
We note that the use of adaptive sampling techniques and multi-armed bandits to accelerate algorithms has also had recent successes in other work, e.g., to accelerate the training of Random Forests 227 solve the Maximum Inner Product Search problem [27]. and more (26).
3 Preliminaries and Background
Notation: We consider a dataset X' of size n (that may contain vectors in IP or other objects). Our goal is to find a solution to the k-medoids problem, Equation (3) We are also given a dissimilarity function d that measures the dissimilarity between two objects in X. Note that we do not assume a
specific functional form of d. We use [r]] denote the set {I,..., ni,, and A b (respectively, V b) to denote the minimum (respectively, maximum) of and b.
Partitioning Around Medoids (PAM): The original Partitioning Around Medoids (PAM) algorithm IO I consists of two main phases: BUILD and SWAP. In the BUILD phase, PAM iteratively initializes each medoid in a greedy, one-by-one fashion: in each iteration, it selects the next medoid that would reduce the k-medoids clustering loss (Equation (3) the most, given the prior choices of medoids. More precisely, given the current set of medoids Mi {m1,' '. , ), the next point to add as a medoid is: m" arg min ZEX
d(x,, xz) A min d(m', 25)) mn E.M,
The output of the BUILD step is an initial set of the k medoids, around which local search is performed by the SWAP phase. The SWAP phase involves iteratively examining all k(n2- k) medoid- nonmedoid pairs and performs the swap that would lower the total loss the most. More preciscly, with M the current set of & medoids, PAM finds the best nondeddd pair to swap:
a
(m" xx) = arg min (M..))MM(XMMMM
Σ
dla m5))
A min d(mml, 25.) ml (m)
(5)
PAM requires Okk²²) distance computations for the A greedy searches in the BUILD step and Okk?? distance computations for each SWAP iteration [28]. The quadratic complexity of PAM makes it prohibitively expensive to run on large datasets. Nonetheless, we describe the PAM algorithm because it has been observed to have the best clustering loss among heuristic solutions to the k-medoids problem. More recent algorithms, such as BanditPAM (28,, achieve the same clustering loss but have significantly improved complexity of Okk log n.) in each step. Our proposed algorithm, BanditPAMM++ improves upon the computational complexity of BanditPAM by factor of (k)).
Sequential Multi-A Armed Bandits: BanditPAM 28 improves the computational complexity of PAM by converting each step of PAM to multi-armed bandit problem. A multi-armed bandit problem (MAB) is defined as a collection of random variables { {RI.... Rn J. called actions or arms. We are commonly interested in the best-arm identification problem, which the highest mean, i.e., arg max; E R.J. with a given probability of possible error d. Many algorithms for this problem exist, each of which make distributional assumptions about the random variables {R..... , R, }: popular ones include the upper confidence bound (UCB) algorithm and successive elimination. For an overview of common algorithms, we refer the reader to [..
identify the arm with
We define a sequential multi-armed bandit problem to be an ordered sequence of multi-armed bandit problems Q
BB,.... Br) where each individual multi-armed bandit problem B, { [R(...., RC] } has the same number of arms n, with respective, timestep-dependent means j pi.....,44, At each timestep f, our goal is to determine (and take) the best action a; arg max, B[ Crucially, our choices of o, will affect the rewards at future timesteps, i.e., the Ry for L' > L. Our definition of a sequential multi-armed bandit problem is similar to non-stationary multi-armed bandit problems. with the added restriction that the only non-stationarity in the problem comes from our previous actions.

---

We now make few assumptions for tractability. We assume that each R is observed by sampling an element from a set S with S possible values. We refer the values of S as the reference points. where each possible reference point is sampled with equal probability and determines the observed reward. With some abuse of notation, we write (x,))) for the reward observed from arm RÚ when the latent variable from S is observed to be 2'g. We refer to a sequential multi-armed bandit as permutation invariint, or as a SPIMAB (for Sequential, Permutation Invariant Multi-Armed Bandit), if the following conditions hold:
1. For every arm and timestep 1, R [D,, fao, 01,........................ ...........11.... ()-1)) for some known function / and some random variable D; with mean 14; == E(D.] E (xa,),
2. There exists a common set of reference points, S, shared amongst each D;.
3. It is possible to sample each D; in O(1) time by drawing from the points in S without replacement, and
4. f is computable in (1) time given its inputs.
Intuitively, the conditions above require that at each timestep, each random variable R2 is expressible as known function of another random variable D; and the prior actions taken in the sequential multi- armed bandit problem. Crucially, Dj does not depend on the limestep; R is only permitted to depend on the timestep through the agent's ' previously taken actions {ag. .i..... at-1). The SPIMAB conditions imply that if E[D,] lA; is known for each i, then μt = E R is also computable in ((1) time for each i and f. i.e. for each arm and timestep.
BanditPAM: BanditPAM 28 reduces the scaling with of each step of the PAM algorithm by reformulating each step as a best-arm identification problem. In PAM, each of the k BUILD steps has complexity O(n²) and each SWAP iteration has complexity Okn?) In contrast, the complexity of BanditPAM is O(n log n)) for each of the k BUILD steps and Okk log n) for each SWAP iteration. Fundamentally, BanditPA. achieves this reduction in complexity by sampling distance computations instead of using all O(n²) pairwise distances in each iteration. We note that all & BUILD steps of BanditPAM (respectively, PAM) have the same complexity as each SWAP iteration of BanditPAM (respectively. PAM). Since the number of SWAP iterations is usually O(k) ( ([28); see also our experiments in the Appendix (3). most of BandidPAM's 's runtime is spent in the SWAP iterations; this suggests improvements to BanditPAM should focus on expediting its SWAP phase.
4 BanditPAM+++ Algorithmic Improvements to BanditPAM
In this section, we discuss two improvements to the BanditPAM algorithm. We first show how each SWAP iteration of BanditPAM can be improved via a technique we call Virtual Arms (VA). With this improvement, the modified algorithm can be cast as SPIMAB. The conversion a SPIMAB permits second improvement via a technique we call the Permutation-IIvanl VA technique improves only the SWAP phase, the PIC technique improves both the BUILD and SWAP phases. The VA technique improves the complexity of each SWAP iteration by factor of O(k), whereas the PIC improves the wall-clock runtime of both the BUILD and SWAP phases.
Cache (PIC). Whereas the
4.1 Virtual Arms (VA)
As discussed in Section 3 most of the runtime of BanditPAM is spent in the SWAP iterations. When evaluating a medoid-nonmedde pair (m, x.) to potentially swap, BanditPAM estimates the quantity:
1
Amm,2
n
for each medoid m and candidate nonmedoid Ti, where
Al Alm,,, (Ix) = [d(x; It)) min d(m', sr,) 0O
(7)
is the change in clustering loss (Fquation (3) induced on point ær for swapping medoid m with nonmedoid x; the set of medoids M. Crucially, we will find that for a given F., each Alm,a, (x.)
for 1, I,...,h, k, except possibly one, is equal. We state this observation formally in Theorem
5

---

General SPIMAB Term
BanditPAM++ BUILD step
BanditPAM++ SWAP step
Arms, { [
Candidate points for medoids
Points to swap in as medoids
Reference points, S Timestep. (x(),)
Points of dataset X' t-th medoid to be added (t - k--th swap to be performed Distance between 2; and 3s Equation
Points of dataset X'
Distance between and Is Equation
[(D;, (ao, 01,....0--1))
Table 1: BanditPAM++'s 's two phases can each be cast in the SPIMAB framework.
Theorem 1. Let Amm,z, (Ix) be the change in clustering loss induced on point Dg by swapping medoid m with nonmedoid Ti, given in Equation ( with 2g and v, fixed. Then the values Alma, (Ix) for m I,...,k k are equal, except possibly where m is the medoid for reference point Zg.
Theorem is proven in Appendix 2 Crucially, Theorem tells that when estimating Alm,, in Equation (6 for fixed n; and various values of m, we can reuse a significant number of the summands across different indices m (across l medoids). We note that Theorem [ has been observed in alternate forms, e.g., as the "FastPAM1" trick, in prior work [28]. However, to the best of our knowledge, we are the first to provide a formal statement and proof of Theorem and demonstrate its use in an adaptive sampling scheme inspired by multi-armed bandits.
Motivated by Theorem I we present the SWAP step of our algorithm, BanditPAM++, in Algorithm BandidPAM++ uses the VA technique improve the complexity of each SWAP iteration by factor of (k). We call the technique "virtual arms" because it uses only a constant number of distance computations to update each of the le "virtual" arms for each of the "real" arms, where "real" arm corresponds to a datapoint.
4.2 crmmutaatio-Invariinn Caching (PIC):
The original BanditPAM algorithms considers each of the k(n k))
pairs as arms. With the VA technique described in Section 4.1 BanditPA.M++ instead considers each of the n datapoints (including existing medoids) arms in each SWAP iteration. Crucially, this implies that the BUILD phase and each SWAP iteration of BanditPAM++ consider the same set of arms. is this observation, induced by the VA technique, that allows us cast BanditPAM++ as SPIMAB and implement a second improvement. We call this second improvement the Permutation- Invariant Cache (PIC).
We formalize the reduction of BanditPAM++ to a SPIMAB in Table 1 In the SPIMAB formulation of BanditPAM++, the set of reference points S is the same as the set of datapoints X. Each D; is a random variable representing the distance from point x; to one of the sampled reference points and can be sampled in (1) time without replacement. Each μ4; E [D] corresponds to the average distance from point v; to all the points in the dataset X. Each arm R² corresponds to the point that would add to the set of medoids (for / < k) or swap in to the of medoids (for &. k), of which there are n at each possible timestep. Similarly, the actions {ao.... ap...., ar} correspond to points added the set of medoids (for S k) or swaps performed (for > k). Equation (8) provides the functional forms of f for each of the BUILD and SWAP steps.
f(Di(x)).A) (d(z1)2)) min d(m', d(m',zs) A A. m'EM
(8)
where M is set of medoids. For the BUILD step. A is sequence of actions that results in a set of medoids of size < k, and, for the SWAP step, A is a set of actions that results in a set of medoids M of size k.
The observation that BanditPAM++ is SPIMAB allows us to develop an intelligent cache design, which we call permutation-imvariann cache (PIC). We may choose a permutation π of the reference points S = X' and sample distance computations to these reference points the order of the permutation. Since we only need to sample some of the reference points, and not all of them, we do not need to compute all (n²) pairwise distances. Crucially, we can also reuse distance computations across different steps of BanditPAM++ to save computational cost and runtime.

---

Algorithm I Bandi SWAP Step ( (fz(D, {ai,.... 0r ), 6, ox, permutation of [ni]
1: Ssalation f- [n]
> Set of potential solutions to MAB
2: t' -0
Number of reference points evaluated
3: For all (1,5) E [rd] [k), set fr.y +- 0, C2,3 t- oo
D Initial means and CIs for all swaps
4: while t' < n and [Sootution] > 1 do
5:
+<)((()
Uses PIC
6: for i E do
7:
Let c(s) and (s) be the indices of v, 's closest and second closest medoids > Cached
8:
Compute distance to 's closest medoid d1 := d(mmc(s).2))
> Cached
9:
Compute distance to Tzy's second closest medoid d2 == d( (s): T.a)
> Cached
10:
Compute d, :- d v.)
Reusing x, 's across calls leads to more cache hits
Pi.c(s) - t.Pe.c di min(da.d.)
11:
>> Update running mean for v. 's medoid
12:
Ci,c(s) + 0.V
> Update confidence interval for 2 's medoid
F+1
13:
for all E {l,..., k)
{c(s)) do
14:
Pig f * (1D, h. )......)))
Update running means; does not depend on
P'11
1o
15:
Ols t -
> Update confidence intervals; does not depend on
241
16:
Ssolution + (i: Aj s.t. jag Cr. I mine, Ct.))
Filter suboptimal arms
17:
-t'- +1
18:
if
|Solution|
then
19:
return i* E Soolution and j arg miny i
20:
else
21:
Compute Pay exactly for all Ssolution
D At most 3n distance computations
22:
return (i" (²)) arg min (..) H4.
The full BanditPAM++ algorithm is given in Algorithm 1 Crucially, for each candidate point to swap into the set of medoids Line we only perform 3 distance computations (not ke) to update all k arms, each of which has mean and confidence interval (CI), on Lines 13-15. This is permitted by Theorem M and the VA technique which says that k "virttall" arms for a fixed will get the same update. The PIC technique allows to choose permutation of reference points (the x, 's) and reuse those n, 's across the BUILD and SWAP steps; as such, many values of d(2,, x.) can be cached.
We emphasize that BanditPAM++ uses the same BUILD step as the original BanditPAM algorithm, but with the PIC. The PIC also used during the SWAP step of BandiiPAM++, is the VA technique. We prove that the full BanditPAM+ algorithm returns the same results as BanditPAM and PAM in Section and demonstrate the empirical benefits of both the PIC and VA techniques in Section 6
Analysis of the Algorithm
In this section, we demonstrate that, with high probability, BanditPAM++ returns the same answer to the m-medoids clustering problem as PAM and BanditPAM while improving the SWAP complexity of BanditPAM by O(k) and substantially decreasing its runtime. Since the BUILD step of BanditPAM is the same as the BUILD step of BanditPAM++, it is sufficient to show that each SWAP step of BanditPAM M++ returns the same swap as the corresponding step of BanditPAM (and PAM). All of the following theorems are proven in Appendix 2
First, we demonstrate that PIC does not affect the results of BanditPAM++ in Theorem S
Theorem 2. Let X (ZI...., J be the reference points of D., and let be a random p (1,..., S}. Then for any S S. ΣG=ı (.((x(p))) has the same distribution as Σ Dr(Ip,). where each PG is drawn uniformly without replacement from 5).
permutation
Intuitively, Theorem 2 says that instead of randomly sampling new reference points at each iteration of BanditPAM++, we may choose a fixed permutation in advance and sample in permutation order at each step of the algorithm. This allows us to reuse computation across different steps of the algorithm.
7

---

We now show that BandidPAMM+ returns the same result as BanditPAM (and PAM) in every SWAP iteration and has the same complexity in as BanditPAM. First, we consider single call to Algorithm| T Let μ4; == min [ii. and let i" == arg min B4 be the optimal point to swap in to the set of medoids, so that the medoid to swap out is j == arg min jE [ Mi .j. For another candidate
point i E [n] with i + i", let A; := p4 py, and for i = i", let A, == minj (2) Ma.j minj Miji
where minj (2) denotes the second smallest value over the indices j. To state the following results, we will assume that, for fixed candidate point and a randomly sampled reference point 2s, the random is o;-sub--aaa aub-Gaussian for some known parameter oj (which, in practice, can be estimated from the data (28)):
variable (z,) A)
Theorem 3. For 5 = 1 1/kr². with probability at least 1 - 2 Algorithm ! to perform using total of M distance computations, where
returns the optimal swap
n'
E|M] < 6n +
Σ min RE(M]
(a, + oo, log kn + B,3n A
12
Intuitively, Theorem 3 states that with high probability, each SWAP iteration of BanditPAM+ returns the same result as BanditPAM and PAM. Since the BUILD step of BanditPAM++ is the same as the BUILD step of BanditPA.M, this implies that BanditPAM++ follows the exact same optimization trajectories as BanditPA.M and PAM over the entire course of the algorithm with high probability. We formalize this observation in Theorem
Theorem 4. If BanditPAM++ ++ is run on a dataset X' with 5 = 1/kn², then it returns the same set of medoids as PAM with probability 1- o(1). Furthermore, the total number of distance computations Mootal required satisfies
E|Mrotal] = O (n log kn).
Note on assumptions: For Theorem B] we assumed that the data is generated in way such that the observations f (Di( A) follow sub-Gaussian distribution. Furthermore, for Theorem { we assume that the A 's are not all close to 0, i.e., that we are not in the degenerate arm setting where many of the swaps are equally optimal, and assume that the a, 's are bounded (we formalize these assumptions in Appendix(2) 2). These assumptions have been found to hold in many real-world datasets 28); see Section and Appendices 1.1 and 2 for more formal discussions.
Additionally, we assume that both BanditPAM and BanditPAM++ place hard constraint on the maximum number of SWAP iterations that are allowed. While the limit on the maximum number of swap steps T may seem restrictive, it is not uncommon to place a maximum number of iterations on iterative algorithms. Furthermore, T has been observed empirically to be ((k) (25). consistent with our experiments in Section Dection[Band Appendix 3
We note that statements similar to Theorems 3 and can be proven for other values of 5. We provide additional experiments to understand the effects of the hyperparameters T and 8 in Appendix 3
Complexity in k: The original BanditPAM algorithm scales as (kc(k)) log kn). where c(k) is deeleetdependeet function of k. Intuitively, c(k) governs the "hardness" of the problem as a function of k; as more medoids are added, the average distance from each point to its closest medoid will decrease and the arm gaps (the A'ss) will decrease, increasing the sample complexity in Theorem 3 With the VA technique, BanditPAM++ removes the explicit factor of k: each SWAP iteration has complexity (cc(k)n log kn). The implicit dependence k may still enter through the term e(k) for a fixed dataset, which we observe in our experiments in Section
Empirical Results
Setup: BandidPAM++ consists of two improvements upon the original BanditPAM algorithm: the VA and PIC techniques. We measure the gains of each technique by presenting an ablation study in which we compare the original BanditPAM algorithm (BP), BanditPAM with only the VA technique (BP+VA), BanditPAM with only the PIC (BP+PIC), and BanditPAM with both the VA and PIC techniques (BP++, the final BanditPAM++ algorithm).
First, we demonstrate that all algorithms achieve the same clustering solution and loss as BanditPAM across a variety of datasets and dataset sizes. In particular, this implies that BanditPAM+ matches
8

---

Dataset Size (n):
10,000
15,000
20,000
25,000
30,000
MNIST (L2, k 10)
1.00
1.00
1.00
1.00
1,00
CIFARIO (L1, k 10)
1.00
1.00
1.00
1.00
1.00
20 Newsgroups (cosine, A = 5)
1.00
1.00
1.00
1.00
1.00
Table 2: Clustering loss BanditPAM++, normalized to clustering loss of BanditPAM, across
variety of datasets, metrics, and dataset sizes. In all scenarios, BanditPAM++ matches the loss of BanditPAM (in fact, returns the exact same solution).
prior state-of -aate-of-hhe-att in clustering quality. Next, we investigate the scaling of all four algorithms in both and across a variety of datasets and metrics. We present our results both sample-complexity and wall-clock runtime. BanditPAM+ M++ outperforms BanditPAM by up to 10x. Furthermore, our results demonstrate that each of the VA and PIC techniques improves the runtime of the original BanditPAM algorithm.
For an experiment on dataset of size n, sampled n datapoints from the original dataset with replacement. In all experiments using the PIC technique, we allowed the algorithm to store up to 1, 000 distance computations per point. For the wall-clock runtime and sample complexity metries, we divide the result of each experiment by the number of swap iterations +1, where the +1 accounts for the complexity of the BUILD step.
Datasets: We conduct experiments on several public, real-world datasets to evaluate BanditPA.M++'s performance: the MNIST dataset, the CIFAR1O dataset, and the 20 Newsgroups dataset. The MNIST dataset [33 contains 70,000 black-and-white images of handwritten digits. The CIFARIO dataset 12 comprises 60,000 images, where each image consists of 32 32 pixels and each pixel has colors. The 20 Newsgroups dataset [19] consist of approximately 18,000 posts on 20 topics split in two subsets: train and test. We used fixed subsample of 10,000 training posts and embedding them into 385-dimensional vectors using a sentence transformer from HuggingFace (7]. We use the L2, L1, and cosine distances across the MNIST, CIFAR1O, and 20 Newsgroups datasets, respectively.
6.1 Clustering/loss quality
First, we assess the solution quality all four algorithms across various datasets, metrics, and dataset sizes. Table 2 shows the relative losses of BanditPAM++ with respect to the loss of BanditPAM M; the results for BP+PIC and BP+VA are identical and omitted for elarity. All four algorithms return identical solutions; this demonstrates that neither the VA nor the PIC technique affect solution quality. particular, this implies that BanditPAM++ matches the prior state-of-the-a algorithms, BanditPAM and PAM, in clustering quality.
6.2 Scaling with
Figure D compares the wall-clock runtime scaling with & of BP, BP+PIC, BP+VA, and BP++ same datasets as Figure I Across all data subset sizes, metrics, and values of k. BP++ outperforms each of BP+VA and BP+PIC, both of which in turn outperform BP. As k increases, the performance gap between algorithms using the VA technique and the other algorithms increases. For example, on the CIFAR10 dataset with k 15, BanditPAM++ is over 10x faster than BanditPAMM. This provides empirical evidence for our claims in Section 5 that the VA technique improves the scaling of the BanditPAM algorithm with k.
We provide similar experiments that demonstrate the scaling with n of BanditPAM++ and each of the baseline algorithms in Appendix 3 The results qualitiatively similar to those shown here; in particular, BanditPAM++ outperforms BP+PIC and BP+VA, which both outperform the original BanditPAM algorithm.
Conclusions and Limitations
We proposed BanditPAM++, an improvement upon BanditPAM that produces state-of-the-art results for the k-medoids problem. BanditPAM++ improves upon BandiPPAM using the Virtual Arms

---

MNIST, La. n= 20000
CIFARIO, I.. o 20000
3.0
2.5
2.0
Timee
0.5
6
10 Number of modelta 01)
14
10
12
34
Number of modiaks (1)
(a)
(b)
102
20 Newsgroups, Cosine,
BP
2.75
BF-HIC
BP-Wh
1.52
7.00
025
Number medelds (8)
Figure 1: Average wall-clock runtime versus for various dataset sizes, metrics, and subsample sizes BP++ outperforms BP+PIC and BP+VA, both of which outperform BP. Negligible error bars omitted for clarity.
technique, which improves the complexity of each SWAP iteration by O(k), and the Permutation- Invariant Cache, which allows the reuse of computation across different phases of the algorithm. We prove that BanditPAM++ returns the same results BanditPAM (and therefore PAM) with high probability. Furthermore, our experimental evidence demonstrates the superiority of BanditPA.M++ over baseline algorithms; across a variety of datasets, BanditPAMM+ prior ttate-of-the-ete art while returning the same results. While the assumptions of BanditPAM and BanditPAM++ are likely to hold many practical scenarios, it is important to acknowledge that these assumptions can impose limitations on our approach. Specifically, when numerous arm gaps are narrow, BanditPAM++ might employ a naive and fall back to brute force computation. Similarly, if the distributional assumptions arm returns are violated, the complexity of BanditPA. may be no better than PAM. We discuss these settings greater detail Appendix [.]
up 10x faster than
Acknowledgements
We would like to thank the anonymous Reviewers, PCs, and ACs for their helpful feedback on our paper. Mo Tiwari was supported by a Stanford Interdisciplinary Graduate Fellowship (SIGF) and Standard Data Science Scholarship. The work of Ilan Shomorony was supported in part by the National Science Foundation (NSF) under grant CCF-2046991. Martin Zhang is supported by NIH RO1 MHI 15676.
10

---

References
[1] Amin Aghace, Mehrdad Ghadiri, and Mahdich Soleymani Baghshah. Active distance-based clustering using k-medoids. Pacific-Asia Conference Knowledge Discovery and Data Mining, 9651:253-22 264, 2016.
[2] David Arthur and Sergei Vassilvitskii. K-means++ the advantages of careful seeding. In Proceedings of the eighteenth annual ACM-SIAM symposium Discrete algorithms, pages 1027-1035, 2007.
[3] Vivek Bagaria, Tavor Z. Baharav, Govinda M. Kamath, and David Tse. Bandit-based monte carlo optimization for nearest neighbors. In Advances in Neural Information Processing Systems, pages 3650-3659, 2019.
[4] Vivek Bagaria, Govinda M. Kamath, Vasilis Ntranos, Martin J. Zhang. and David Tse. Medoids almost-linear time via multi-armed bandits. In International Conference on Artificial Intelli- gence and Statistics, pages 500-509, 2018.
[5] Paul S. Bradley, Olvi L. Mangasarian, and W. N. Street. Clustering via concave minimization. In Advances in Neural Information Processing Systems, pages 368-374, 1997.
[6] Vladimir Estivill-Castro and Michael E. Houle. Robust distance-based clustering with applica- tions to spatial data mining. Algorithmica, 30(2):216 2001.
[7] HuggingFace.
all-minilm-16--2 model
card.
attps://huggingee co/
sentence transformers/aaaeeoeooorr www.fformrrraaa-Miiillee.. -Min-MMi-MMM-------- L6-LLM--66--0 L6- v2 Accessed: 10-2-10-28.
[8] Anil K. Jain and Richard C. Dubes. Algorithms for clustering data. Prentice-Hae Hall, 1988.
[9] Kevin Jamieson and Robert Nowak. Best-arm identification algorithms for multi-armed bandits in the fixed confidence setting. In Annual Conference Information Sciences and Systems, pages 1-6, 2014.
[10] Leonard Kaufman and Peter J. Rousseeuw. Clustering by means of medoids. Statistical Data Analysis based on the LI Norm and Related Methods, pages 405-416, 1987.
[11] Leonard Kaufman and Peter J. Rousseeuw. Partitioning around medoids (program pam). Finding groups in data: an introduction to cluster analysis, pages 68 68-125, 1990.
[12] Alex Krizhevsky and Geoffrey Hinton. Learning multiple layers of features from tiny images. Technical report, Citeseer, 2009.
[13] Yann LeCun, Léon Bottou, Yoshua Bengio, and Patrick Haffner. Gradient-based learning applied to document recognition. Proceedings of the IEEE, 86(11):2278-2324, 1998.
[14] Daniel LeJeune, Richard G. Baraniuk, and Reinhard Heckel. Adaptive estimation for approxi- mate k-nearest-ner neighhor computations. In International Conference on Artificial Intelligence and Statistics, 2019.
[15] Jure Leskovec, Anand Rajaraman, and Jeffrey D. Ullman. Mining of massive data sets. Cam- bridge university press, 2020.
[16] Stuart Lloyd. Least squares quantization in pcm. IEEE transactions on information theory, 1892),1 1982.
[17] James MacQueen. Some methods for classification and analysis of multivariate observations. In Berkeley Symposium on Mathematical Statistics and Probability, volume 1, pages 281 281-297, 1967.
[18] Nina Mishra, Robert Schreiber, Isabelle Stanton, and Robert E. Tarjan. Clustering social networks. In International Workshop on Algorithms and Models for the Web-Graph pages 56-67. Springer, 2007.
[19] Tom Mitchell. Twenty Newsgroups. http://too.oog/003444442C3333333
UCI Machine Learning Repository, 1999. DOI:

---

[20] Gonzalo Navarro. A guided tour to approximate string matching. ACM computing surveys, 33(1))31-88, 2001.
[21] Raymond T. Ng and Jiawei Han. Clarans: method for clustering objects for spatial data mining, IEEE transactions knowledge and data engineering, 14(3):1003-1016, 2002.
[22] Michael L Overton. A quadratically convergent method for minimizing norms, Mathematical Programming, 27(1):34-65 63, 1983.
of euclidean
[23] Rameshwar Pratap, Anup Deshmukh, Pratheeksha Nair, and Tarun Dutt. A faster sampling algorithm for spherical k-means. Asian Conference on Machine Learning, pages 343-358.
PMLR, 2018.
[24] Jiayang Ren, Kaixun Hua, and Yankai Cao. Global optimal k-medoids clustering of one million samples. Advances in Neural Information Processing Systems, 35:982-994, 2022.
[25] Erich Schubert and Peter Rousseeuw. Faster k-medoids clustering: improving the pam, clara, and clarans algorithms. In International Conference Similarity Search and Applications, pages 171-187. Springer, 2019.
[26] Mo Tiwari. Accelerating machine learning algorithms with adaptive sampling. arXiv preprint arXiv:2309. 14221, 2023.
[27] Mo Tiwari, Ryan Kang, Jacyong Lee, Chris Piech, Ilan Shomorony, Sebastian Thrun, and Martin J Zhang. Mabsplit: Faster forest training using multi- bandits. Advances in Neural Information Processing Systems, 35: 1223 1237, 2022.
[28] Mo Tiwari, Martin Jinye Zhang, James Mayelin, Sebastian Thrun, Chris Piech, and Ilan Shomorony. Banditpam: Almost linear time k-medoids clustering via multi-armed bandits. In Advances Neural Information Processing Systems, 2020.
12


Reviews:
[{"review_id": "x7RkhVEQSQ", "reviewer": "Reviewer_oPsi", "summary": "The paper \n * presents a new algorithm for the k-medoids problem called BandidPAM++, which is an extension of BanditPAM\n * analyzes the algorithm theoretically showing that it produces essentially the same output as BanditPAM and uses O(n log n) distance computations\n * evaluates the algorithm empirically comparing it to BanditPAM", "questions": "1. Is it true that BanditPAM and BanditPAM++ come with the same running time bound on the BUILD iterations?\n1. Is the running time of each SWAP iteration O(k n log n) for BanditPAM vs O(n log n) for BanditPAM++? (Theorem 4 bounds the number of comparisons which is not the same as the running time)\n1. Do you use any results from the multi-arm bandit literature in any of the proofs?\n1. Since Theorem 1 is a reformulation of a known fact, what is the novelty in the virtual arms optimization?\n1. How similar is Theorem 3 to what was previously done for BanditPAM?\n1. What is the high-level idea behind Theorem 3?\n1. Could you also elaborate on the main idea behind permutation-invariant cache? I find the current description (lines 196-201) to be too compact to be understandable.\n1. Line 131 suggests that each $D_j$ is a random variable which represents a *uniform* distribution over $\\mathcal{S}$. Could you formally define $D_j$?\n1. Does Table 1 describe BanditPAM or BanditPAM++? The caption vs line 183 seem to contradict each other.", "limitations": "Yes", "rating": 7, "confidence": 2, "soundness": 3, "presentation": 2, "contribution": 3, "strengths": "* The paper new algorithm outperforms the baseline in experiments, delivering significantly better running times\n* While I'm not absolutely certain, it appears to me that some of the ideas involved are nontrivial", "weaknesses": "The main issue with the paper is the presentation. In the current form the description of the proposed algorithm consists of a short, high-level description + a nontrivial pseudocode. This is especially disappointing given that it would be easy to fit in more text in the first 9 pages (there are many subsection headers, figures could be made smaller, Algorithm 1 could potentially go to the appendix, etc). This, combined with the fact that I'm not familiar with the area of bandit algorithms, makes it difficult for me to make a confident decision.\n\nMy second (smaller) concern is about the empirical evaluation: while the motivation for using k-medoids is clustering points in non-Euclidean spaces, the experiments use two datasets, both of which define Euclidean space, which is non-ideal (but I guess acceptable for a paper whose main goal is to demonstrate speedups)\n\n[Edit: thank you for addressing this point in the rebuttal]", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer oPsi for their thoughtful comments and provide a point-by-point response below.\n\nWeakness 1: We have now added detailed explanations of the pseudocode at the end of Section 3 and made better use of space. BanditPAM++ improves upon BanditPAM via two techniques: the Virtual Arms (VA) technique, and Permutation-Invariant Caching (PIC). The Virtual Arms technique improves the sample complexity of BanditPAM++ relative to BanditPAM by $O(k)$ in each SWAP. The Permutation-Invariant Caching technique results in improvements to wall-clock time.\n\nThe Virtual Arms technique is implemented in lines 6-11 of Algorithm 2. Specifically, for a given reference point $x_r$ and a candidate point $x_c$, instead of evaluating all $k$ medoids as in BanditPAM (Line 6 of Algorithm 1), BanditPAM++ only performs 3 distance computations (Line 7 of Algorithm 2), which suffices to update the information of all $k$ medoids. This brings an $O(k)$ improvement of algorithm complexity and is permitted by Theorem 1.\n\nPermutation-Invariant Caching is implemented in line 5 of Algorithm 2. Instead of randomly sampling a reference point $x_r$ as in BanditPAM (line 5 Algorithm 1), BanditPAM++ prestores a permutation of reference points in cache and draws the reference point $x_r$ based on this prestored permutation. This procedure substantially increases the speed of BanditPAM++ in wall-clock time. Choosing a permutation of datapoints a priori and using this permutation across all phases of BanditPAM++ increases the proportion of cache hits when the cache is too small to hold all $n^2$ pairwise distances. Without a permutation order, many cache misses would occur because we are not reusing the same reference points (and therefore distances) across different steps of the algorithm.\n\nWeakness 2: The CIFAR experiment in the original submission used $L_1$ distance (lines 241-242 and captions of Figures 1, 3, and 4 of the original submission). \n\nWe have also run new experiments across two different datasets and metrics ($L_1$ distance on scRNA-seq data; cosine similarity on newsgroups data) to demonstrate the speedups that BanditPAM++ has over BanditPAM in non-Euclidean spaces. BanditPAM++ demonstrates significant speedups over BanditPAM (Figure 1 and Table 1 in the additional experiments).\n\nQ1: BanditPAM and BanditPAM++ have the same sample complexity for BUILD, because both use Algorithm 3 (in Appendix 1 of our original submission) for the BUILD phase. However, BanditPAM++ accelerates its BUILD steps using the Permutation-Invariant Caching technique described above to cache the distance computations and reuse them across iterations, so its wall-clock time is better (Table 1 of the additional rebuttal experiments). \n\nQ2: Yes, the running time of each SWAP iteration is $O(k n \\text{log} n)$ for BanditPAM vs $O(n \\text{log} n)$ for BanditPAM++.\n\nWe chose to focus on the number of distance computations in our theoretical analysis because the overall complexity is governed by the number of distance computations [19, 20]. We also evaluated the wall-clock time improvement using real-data experiments. As shown in Figure 3, the improvement of BanditPAM++ over BanditPAM increases with $k$. \n\nFor real-data, wall-clock time of BanditPAM++ also slightly increased with $k$, possibly because the MAB problem instances became harder (smaller $\\Delta_i$'s) as $k$ increased (assumed to be the same in theoretical analysis). Nonetheless, this will have the same impact on both algorithms. Therefore, the improvement of BanditPAM++ over BanditPAM is still $O(k)$. \n\nQ3: Yes, we borrowed some proof techniques from the multi-armed bandit literature. Algorithms 1 and 2 (as well as Algorithm 3 in Appendix 1) are similar to the successive elimination algorithm from multi-armed bandits. Our proof of Theorem 3, which proves correctness and the sample complexity of BanditPAM++, is similar to the proof of correctness and sample complexity for successive elimination. Additionally, the proof of Theorem 4 relies on proofs from the multi-armed bandit literature to convert a gap-dependent bound that depends on the $\\Delta_i$'s to a bound in terms of the dataset size $n$ [2, Appendix 2 of 3].\n\nQ4: Though Theorem 1 has been described elsewhere [19], our paper provided its first formal proof.  More importantly, it is non-trivial to incorporate the observation in Theorem 1 into a randomized algorithm, BanditPAM, to achieve a complexity improvement. Before our work, it was unclear that the observation in Theorem 1 could be applied to a randomized algorithm for algorithmic speedups.\n\nQ5: The proofs are different because BanditPAM considers datapoint-medoid pairs as separate “arms” whereas BanditPAM++ considers the $k$ medoids of a given datapoint share the same “arm” and are updated together. As a result, the proof from BanditPAM does not carry over to BanditPAM++. \n\nQ6: The  idea behind Theorem 3 is that BanditPAM++ is a) correct, and b) sample-efficient. Regarding a), Theorem 3 shows that with high probability, BanditPAM++ returns the same results as BanditPAM (and therefore PAM). a) is achieved by a union bound over the error probability of all confidence interval estimates. Regarding b), the bound is a sum of sample complexity over all datapoints (arms). Specifically, for a datapoint $x$, Theorem 3 states that its sample complexity depends on the difference between its quality and that of the best arm, i.e., the gap $\\Delta_x$'s. \"Bad\" arms with high $\\Delta_x$ will be thrown out with just a few samples, whereas \"good\" arms with low $\\Delta_x$ will be sampled more to determine the true best arm. \n\nQ7: Please see response to Weakness 1.\n\nQ8: We define $D_j$ as an r.v. which has a uniform distribution over the possible distances from point $j$ to any point in the dataset; in this way, $E[D_j] = \\mu_j$ is the average distance from point $j$ to all points in the dataset.\n\nQ9: Line 183 should read \"BanditPAM\" instead of \"BanditPAM++\"; we have fixed this."}, {"review_id": "l7AvIbxGfE", "reviewer": "Reviewer_zP71", "summary": "The authors developed an algorithm (called BanditPAM++) which is claimed to be more efficient than banditPAM, the state-of-the-art algorithm for k-medioid (which is the version of the k-means problem where the centers must be input points). They prove that with high probability the two algorithms produce the same results, with BanditPAM++ being more efficient. An experimental evaluation on two datasets complements their results. ", "questions": "Questions:\n- from a theoretical point of view, what is the advantage of BanditPAM++ over BanditPAM? (both have running time O(n logn) as far I could check). Update after rebuttal: I had missed that BanditPAM++ is faster than BanditPAM by a factor of k. This should be stressed in the text. Perhaps a short paragraph right after Theorem 4 might help. \n- discuss the role of delta and the max number of swaps T and how they might affect the results\n\nother comments:\n- Figure 4 does not offer many insights, all approaches are superimposed and we just see a straight line.\n- there is a version of k-means with cosine similarity, in contrast with what is claimed in the paper (references [4,6,18] are old and not up to date)\n\nA Faster Sampling Algorithm for Spherical k-means, Journal of machine learning 2018. \n\n- I recommend to evaluate the impact of delta and T in the experiments, for future work\n- in the experiments it should be immediately clear which curves correspond to BanditPAM++ (all of the variants? some?)", "limitations": "No potential negative societal impact, as far as I could check. ", "rating": 4, "confidence": 4, "soundness": 3, "presentation": 2, "contribution": 2, "strengths": "- the work focuses on improving the scalability of an effective algorithm for clustering\n- theoretical results support their claims, which are validated empirically on two datasets", "weaknesses": "- it is not clear how some of the parameters (in particular T=max number of swaps and delta) affect the results.\n- only two datasets are included in the experimental evaluation. I recommend a more extensive evaluation, given that it is not clear what is the advantage of their approach from a theoretical point of view\n- the analysis focuses on the case when the input points follow a subgaussian distribution, while for example k-means++ guarantees hold for any distribution\n- the presentation of the paper could be improved (e.g. the role of delta is not discussed, not immediately clear which curve in the experiments correspond to BanditPAM++, some claims are inaccurate, e.g. there is a version of k-means that can deal with cosine similarity)\n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer zP71 for their thoughtful comments and helpful feedback. We provide a point-by-point response below.\n\nResponse to W1 and Q1: From a theoretical point of view, the computational complexity of BanditPAM++ is same as that of BanditPAM with respect to the number of datapoints $n$ (both $O(n$ log $n)$), but it is smaller than that of BanditPAM with respect to the number of medoids $k$,  ($O(1)$ vs. $O(k)$) in each SWAP. Overall, BanditPAM++ is $O(k)$ faster than BanditPAM in each SWAP step.  \n\nFor real-data, wall-clock time of BanditPAM++ may slightly increase with $k$, possibly because the MAB problem instances became harder (smaller $\\Delta_i$'s) as $k$ increases (assumed not to change in theoretical analysis). Nonetheless, this will have the same impact on both algorithms. Therefore, the improvement of BanditPAM++ over BanditPAM is still $O(k)$.\n\nIn addition, on real-data experiments, we verified that BanditPAM++ is substantially faster than BanditPAM while producing same results and that this improvement increased with $k$ (e.g., $4.75\\times$ faster with $k=5$, $8.23\\times$ faster with $k=10$, and $10.77\\times$ faster with $k=15$ on the same dataset across various dataset sizes; Table 1 of the additional rebuttal experiments. Figure 3b of the original submission and Figure 1 of the additional rebuttal experiments show similar trends). \n\nResponse to W2 and Q2: We have run additional experiments to understand the effects of $T$ and $\\delta$ and provided them in the \"global\" response to all reviewers.\n\nWe note that the choice of hyperparameter $\\delta = n^{-3}$ is not very important and is only chosen to show that the correctness probability is $1 - o(1)$; an analogous claim can be made for arbitrary $\\delta$ (Remark A1 on Lines 407-40). In practice, even setting $\\delta$ very large does not affect the quality of our results. We present an additional experiment demonstrating this in the \"global\" rebuttal. Table 2 of our additional results demonstrates that on two different datasets and metrics, $\\delta$ can be varied significantly without affecting the quality of the clustering results.\n\nNote that increasing $T$ will decrease the overall loss of the clustering solution as more swaps can be performed; however, we find that beyond $T = k$, increasing $T$ further does not lower the overall loss very significantly. This is specified on Line 121 of our original submission and has been observed in prior work [19, 20]. Figure 2 of our additional experiments shows the loss vs. swap iteration ($T$) for the MNIST dataset with $L_2$ distance and the CIFAR10 dataset with $L_1$ distance, for $k=5, 10$. Increasing $T$ beyond $k$ only leads to marginal decreases in overall loss.\n\nResponse to W3: We emphasize that from a theoretical point of view, BanditPAM++ is $O(k)$ faster than BanditPAM in each SWAP. \n\nWe have also run new experiments across two different datasets and metrics ($L_1$ distance on scRNA-seq data; cosine similarity on newsgroups data) to demonstrate the speedups that BanditPAM++ has over BanditPAM in non-Euclidean spaces. On both new datasets and metrics, BanditPAM++ demonstrates significant speedups over BanditPAM. Please see Figure 1 and Table 1 in the additional experiments.\n\nResponse to W4: We followed the same assumptions of the original BanditPAM paper [20]. The subgaussian assumption is a reasonable assumption in many real-world datasets; for example, these assumptions hold in any bounded dataset. Furthermore, prior work [20, Appendix 1.3] has shown that these assumptions hold in most practical datasets, e.g., the MNIST dataset with either $L_1$ or $L_2$ distances and the scRNA dataset with $L_1$ distance. Similar observations have been made in [3, Figures 2 and 3]\n\nThese assumptions are also common in the multi-armed bandit literature [2, 3]. The scaling of our algorithms in Section 5 also demonstrates that the sub-Gaussianity assumptions hold in practice. \n\nHowever, we also acknowledge that more general assumptions may be more desirable, e.g., for unbounded or heavy-tailed distributions. We added a discussion of this limitation in the revised paper. \n\nResponse to W5 (curves): The algorithm we refer to as BanditPAM++ is the version of BanditPAM with both the Virtual Arms (VA) technique and the Permutation-Invariant Caching (PIC) technique, and corresponds to the red lines in all plots. This is our best-performing algorithm. \"BanditPAM Original without caching\" refers to the original BanditPAM algorithm without the VA or PIC techniques, \"BanditPAM Original with caching\" refers to the original BanditPAM algorithm with a cache but without the VA technique, and \"BanditPAM VA without caching\" refers to the original BanditPAM algorithm with the VA technique but without PIC.\n\nWe have clarified the legends in the additional experiments presented in the \"global\" rebuttal and will clarify the plots in the original submission in our final version.\n\nResponse to W5 and Q4 (k-means comments): Thank you for pointing this out. We will clarify this and cite \"A Faster Sampling Algorithm for Spherical k-means, Journal of machine learning 2018.\" in our final paper. \n\nWe also would like to note that the general point is still valid: $k$-means typically requires specific distance metrics for efficient solutions. In contrast, $k$-medoids is agnostic to the distance metric being used. In fact, the distance function being used in $k$-medoids need not even be a metric and can be an arbitrary dissimilarity function (i.e., may violate the triangle inequality, may be negative, and may not be symmetric). For this reason, $k$-medoids has been used to cluster \"exotic\" objects such as trees [20].\n\nResponse to Q3: Thank you for this feedback. We will convert Figure 4 to two much smaller tables that show the losses of each algorithm relative to BanditPAM, where the entries of each table will simply be \"1.00x\". This is also what we have done for Table 2 of our additional rebuttal experiments."}, {"review_id": "ynB0ojx7oh", "reviewer": "Reviewer_hYVb", "summary": "This paper studies the $k$-mediods clustering problem. Given an instance $P$, the goal of the $k$-mediods problem is to find a subset $C \\subseteq P$ of size $k$ such that the sum of the defined distances from data points to $P$ to the centers in $C$ is minimized. In previous heuristic solutions for the $k$-mediods problem, the state-of-the-art $k$-mediods algorithm is the BanditPAM method based on the multi-armed bandit strategy which can return the same solution as PAM (the best heuristic algorithm using exhaustive searching) with high probability and time complexity $O(nlogn)$ and $O(knlogn)$ for each BUILD and SWAP iteration, respectively. In this paper, an improved algorithm is proposed (called BanditPAM++) with faster running time and can return the same solution to BanditPAM with high probability. By designing Virtual Arms (VA) technique and using a permutation-invariant cache (PIC) structure, this paper shows that the total number of distance computations in expectation can be reduced to $O(nlogn)$, which is much faster than BanditPAM method. Experimental results show that the proposed BanditPAM++ method runs 10 times faster than BanditPAM method with the same clustering results on some of the commonly used datasets for clustering tasks.", "questions": "1. Can the authors provide detailed comparison between the running time of BanditPAM++ and BandiPAM for both BUILD iteration and SWAP iteration.\n\n2. Does the maximum iteration $T$ required for BanditPAM and BanditPAM++ to reach a good convergence the same?\n\n3. Why the authors claim that the assumptions of BanditPAM and BanditPAM++ are likely to hold in many practical scenarios?\n\n4. If the sub-gaussian distribution assumptions do not hold, what will the guarantee be like (maybe some worse guarantee)?", "limitations": "I don't think this paper has potential negative societal impact of their work. ", "rating": 4, "confidence": 2, "soundness": 2, "presentation": 2, "contribution": 2, "strengths": "1. This paper proposes a new heuristic method for the $k$-mediods clustering problem. The authors show that  the proposed method runs much faster than the current SOTA (BanditPAM method) and can give the same solution as the current SOTA with high probability. \n\n2. This paper proposed two novel techniques and structures, namely the Virtual Arms (VA) technique and the permutation-invariant cache (PIC), which is proved to be efficient for accelerating the BanditPAM method in both theory and practice.\n\n3. The experimental results show that the proposed BanditPAM++ method has significant advantages on the running time compared with BanditPAM method, while the clustering quality of BanditPAM++ matches that of BanditPAM", "weaknesses": "1. The theoretical proofs of this paper is mainly based on the assumptions that (1): during the SWAP process, the samples are drawn from a sub-Gaussian distribution. (2) the maximum number of swap steps is limited. However, there is no detailed discussion about the assumptions made (especially the sub-Gaussian distribution assumption). Although the authors mentioned in the conclusions and limitations part that \"while the assumptions of BanditPAM and BanditPAM++ are likely to hold in many practical scenarios, it is important to acknowledge that these assumptions can impose limitations on our approach\", there is no clear explanation about the assumptions made and it is unclear why the Gaussian distribution is likely to hold in many practical scenarios. Besides, in a clustering algorithm, the iteration rounds needed to reach the convergence is not usually a constant $T$. For example, in local search algorithms for clustering, the iteration rounds could be as large as $O(klog\\Delta)$ if the initialization is bad, where $\\Delta$ is the aspect ratio of the given instance. Thus, it is unlcear how the number of swap steps could affect the results of BanditPAM++\n\n2. For the proposed BanditPAM++ method, this paper claims that it runs much faster than BanditPAM. However, there is no direct comparison between each phase (i.e., the BUILD and SWAP iteration). It is unclear how the time complexity is improved compared with BanditPAM method where each BUILD iteration requires $O(nlogn)$ running time and each SWAP iteration requires $O(nklogn)$ running time.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer hYVb for their thoughtful comments and helpful feedback. We provide a point-by-point response below.\n\nResponse to W1, Q3, and Q4 (sub-Gaussianity): We agree that it is ideal to remove the sub-Gaussian assumption and discuss this limitation in the revised paper. However, we argue this is a reasonable assumption for many real-world datasets. First, this assumption holds in any bounded dataset, since a distribution with bounded support is always $\\sigma$-sub-Gaussian for some $\\sigma$. Furthermore, prior work [20, Appendix 1.3] has shown that these assumptions hold in most practical datasets, e.g., the MNIST dataset with either $L_1$ or $L_2$ distances and the scRNA dataset with $L_1$ distance. Similar observations have been made in [3, Figures 2 and 3]\n\nOur assumptions also follow that of the original BanditPAM paper [20], so are as general as in the original BanditPAM algorithm. These assumptions are also common in the multi-armed bandit literature [2, 3]. The scaling of our algorithms in Section 5 also demonstrates that the sub-Gaussianity assumptions hold in practice. We acknowledge, however, that these assumptions may introduce a limitation in pathological datasets, e.g., datasets that are unbounded.\n\nIf the sub-Gaussianity assumptions do not hold, it is possible that, in the worst case, BanditPAM++ will not identify the correct medoids. However, this is the same behavior as BanditPAM and in practice, we do not observe that this happens. Our results in Section 5 and prior work [20] demonstrate that BanditPAM++ and BanditPAM return the same results as PAM.\n\nResponse to W1 and Q2 ($T$): We have run additional experiments to understand the effect of $T$ and provided them in the \"global\" response to all reviewers.\n\nNote that increasing $T$ will decrease the overall loss of the clustering solution as more swaps can be performed; however, we find that beyond $T = k$, increasing $T$ further does not lower the overall loss very significantly. This is described on Line 121 of our original submission and has been observed in prior work [19, 20].\n\nWe have provided new empirical evidence for this claim in the additional experiments in the \"global\" rebuttal. Figure 2 of our additional experiments shows the loss vs. swap iteration ($T$) for the MNIST dataset with $L_2$ distance and the CIFAR10 dataset with $L_1$ distance, for $k=5, 10$. Increasing $T$ beyond $k$ only leads to marginal decreases in overall loss. BanditPAM++ and BanditPAM track the same optimization trajectory, so are superimposed.\n\nFinally, we note that setting a maximum number of iterations $T$ for iterative algorithms is common practice, e.g., in $k$-means. We are aware of convergence results for $k$-means, (e.g., [F]), however it is not immediately clear how these carry over to $k$-medoids. We acknowledge that we do not provide a strict theoretical bound on $T$ for the convergence of $k$-medoids algorithm, however, to the best of our knowledge no such bound exists for PAM or BanditPAM either (which would imply a bound for BanditPAM++, since they all follow the same optimization trajectory). Proving such a bound on $T$ has been observed to be a difficult problem [19, 20].\n\nWe also note that the number of iterations required for BanditPAM and BanditPAM++ to reach a good solution are the same; this is because BanditPAM++ is guaranteed to return the same solution as BanditPAM to every SWAP step. This is the key point of Theorem 3, which states that BanditPAM++ follows the same optimization trajectory as BanditPAM and PAM. However, we agree that this is currently unclear as written because Theorem 4 states that only the final set of $k$ medoids is the same. We will clarify that BanditPAM++ tracks the exact same optimization trajectory as BanditPAM and PAM -- in particular, for each BUILD and SWAP step -- in the final paper.\n\n[F: \"A theoretical analysis of Lloyd's algorithm for k-means clustering\" Bhowmick 2009]\n\nResponse to W2 and Q1: Thank you for the comment. We have run additional experiments to understand the speedups in the BUILD step only, the SWAP steps only, and in the BUILD + SWAP steps. We present our results in Table 1 in the additional experiments in a \"global\" response to all reviewers.\n\nWe observe that the gains for the BUILD step are slightly less than those for the SWAP steps; this is expected because the Virtual Arms (VA) technique improves only the SWAP step and not the BUILD step. However, we still observe significant gains in the BUILD step due to the Permutation-Invariant Caching (PIC) technique, which caches distance computations and reuses them across iterations of the BUILD step. For example, on the scRNA dataset with $L_1$ distance, the BUILD step of BanditPAM++ is still over $3\\times$ faster than BanditPAM with $k=15$. On the 20 Newsgroups dataset with cosine distance, the BUILD step of BanditPAM++ is still $3\\times$ faster than the BUILD step of BanditPAM when $k=15$.\n\nIn the same table, we also provide results analyzing the improvements to the SWAP step alone, as well as the overall algorithm (BUILD + SWAP).\n\nOn both new datasets and metrics, we analyze the scaling of BanditPAM++ with respect to both $n$ and $k$. In both additional settings, BanditPAM++ demonstrates significant speedups over BanditPAM. On the scRNA dataset, BanditPAM++ is approximately 6.4x faster than the original BanditPAM algorithm with $k=15$ and $L_1$ distance, even with $n = 10,000$ on this moderately-sized dataset. On the 20 newsgroups dataset with $k = 15$ and cosine distance (a highly non-Euclidean space), BanditPAM++ is up to 8.5x faster than the original BanditPAM algorithm with $n = 10,000$.\n\nWe will include a full description of datasets in the final paper.\n"}, {"review_id": "hHi6KNs5uh", "reviewer": "Reviewer_s8tF", "summary": "The paper introduces BanditPAM++, a faster version of the BanditPAM algorithm for k-medoids clustering. It incorporates two improvements: the Virtual Arms technique, which reuses loss changes within each SWAP iteration, and the Permutation-Invariant Caching method, which samples reference points and reuses distance computations across steps. The authors provide theoretical guarantees that BanditPAM++ produces the same solution as BanditPAM and showcase its significant speedup on multiple datasets.", "questions": "See the weakness", "limitations": "The exact speedup of this method cannot be theoretically derived; only empirical speedup can be obtained. Furthermore, the authors have highlighted cases in which the VA technique fails, as mentioned in section 6.", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "The originality of the paper lies in the combination of techniques from the multi-armed bandit literature with k-medoids clustering, resulting in a novel algorithm that achieves significant speedup. Most part of the paper is well written and organized.", "weaknesses": "1. The presentation of empirical results is somewhat confusing. Here are a few specific concerns:\n\n1.1 The names of methods are not consistent between the text and figures. For example, \"BanditPAM\" in Line 234 is referred to as \"BanditPAM Original with caching\" in the figures. It would be helpful to ensure consistency in naming.\n\n1.2 It seems that the results in Figure 1 and Figure 3 are intended to reflect the same experimental setting, but there are inconsistencies. For instance, in the MNIST case with n=70000, the runtime of each method differs between the two plots. Additionally, in the CIFAR10 case, the runtime for k=10 in Figure 3 is significantly smaller than that in Figure 1. It would be beneficial for the author to review and revise the results or provide an explanation for these discrepancies.\n\n1.3 Since the specific speedup of BanditPAM++ can only be evaluated based on the empirical results, it would be clearer to include a table presenting the specific values of speedup. For example, adding a table of average speedup over all dataset sizes on different k. This approach would make it easier for readers to grasp the key points compared to the current figures, which require a non-trivial comparison between several plots.\n    \n2. In addition to heuristic methods, there is currently a growing trend in solving clustering problems, including K-Medoids, using global optimization approaches. It would be valuable for the author to discuss these methods in the related work section and provide a comparison between global and heuristic methods. This discussion could highlight the advantages and limitations of each approach, shedding light on the current state-of-the-art in clustering techniques.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer s8tF for their thoughtful comments and helpful feedback. We provide a point-by-point response below.\n\nResponse to W1.1: We apologize for the confusion. The algorithm we refer to as BanditPAM++ is the version of BanditPAM with both the Virtual Arms (VA) technique and the Permutation-Invariant Caching technique, and corresponds to the red lines in all plots. This is our best-performing algorithm. \"BanditPAM Original without caching\" refers to the original BanditPAM algorithm without the Virtual Arms or Permutation-Invariant Cache techniques, \"BanditPAM Original with caching\" refers to the original BanditPAM algorithm with a cache but without the Virtual Arms technique, and \"BanditPAM VA without caching\" refers to the original BanditPAM algorithm with the Virtual Arms technique but without a cache.\n\nWe have clarified the legends in the additional experiments presented in the \"global\" response to reviewers and will clarify the plots in the original submission in our final version.\n\nResponse to W1.2: We apologize for the confusion. The reason for this discrepancy is that each set of experiments was run on a different machine, each with different hardware specifications. When run on the same machine, these discrepancies disappear; we will re-run all experiments on a single machine and update the plots for the final version.\n\nAdditionally, we have rerun experiments on a single machine for $n=10000$ and $k=5, 10, 15$ on the MNIST dataset with $L_2$ distance, the CIFAR10 dataset with $L_1$ distance, and for two new datasets: the scRNA dataset with $L_1$ distance, and the 20 Newsgroups dataset with cosine distance; our results are presented in Table 1 in the additional experiments in a \"global\" response to all reviewers. The general trend still persists and we observe that BanditPAM++ is significantly faster than the original BanditPAM algorithm, sometimes even an order of magnitude faster.\n\nOn both new datasets and metrics, we analyze the scaling of BanditPAM++ with respect to both $n$ and $k$. In both additional settings, BanditPAM++ demonstrates significant speedups over BanditPAM. On the scRNA dataset, BanditPAM++ is approximately 6.4x faster than the original BanditPAM algorithm with $k=15$ and $L_1$ distance, even with $n = 10000$ on this moderately-sized dataset. On the 20 newsgroups dataset with $k = 15$ and cosine distance (a highly non-Euclidean space), BanditPAM++ is up to 8.5x faster than the original BanditPAM algorithm with $n = 10000$. Please see Figure 1 of the additional experiments in the \"global\" rebuttal for more details.\n\nWe will include a full description of datasets in the final paper.\n\nResponse to W1.3: Thank you for this feedback. We have created a Table of the suggested format in Table 1 of our additional experiments and included it in the \"global\" response to all reviewers. We agree that this is clearer than the plots as currently shown and will include these tables in the final paper. We will also convert the relevant plots in the initial submission to tables when we have enough time to rerun all experiments for different values of $n$.\n\n\nResponse to W2: We agree with this feedback and will include global optimization approaches, (e.g., [D] and [E]) in our related work. Our method is crucially different from global optimization approaches such as [D] and [E] in that those methods scale at least quadratically in space and/or time with dataset size. In contrast, our method scales almost linearly in both space and time with dataset size. As such, our heuristic approach scales to even larger datasets, at the expense of a possibly worse optimality gap. We will include this discussion in the final paper.\n\n[D: Global Optimal K-Medoids Clustering of One Million Samples, Jiayang Ren, Kaixun Hua, by Yankai Cao]\n[E: K-Medoids Clustering Is Solvable in Polynomial Time for a 2d Pareto Front]\n"}, {"review_id": "QxXLKeFXNk", "reviewer": "Reviewer_mktW", "summary": "This paper studies a clustering problem (k-medoid) in which we have to select k centers among the datapoints in order to minimize the sum of distances of points to their closest center. For this problem, there is an intuitive algorithm which works well empirically (the PAM algorithm): (1) add greedily centers one by one by selected the next center to be the one which maximizes the decrease in cost (BUILD phase), (2) once we have selected k-centers, we try to swap one center for another datapoint again greedily (SWAP phase). We stop after some iterations (usually we stop after O(k) iterations). This algorithm works well in practice but is quite expensive in terms of computation. Therefore, several works have designed variants of the PAM algorithm to lower the computational cost without sacrificing too much of the quality of the solution. BanditPAM is one of these algorithm which returns the same solution as PAM (with high probability) but performs much less distance computations. \n\nThis paper further improve BanditPAM to propose BanditPAM++. They show that BanditPAM++ also returns the same solution as PAM with high probability, and uses less comparisons compared to BanditPAM. On some real-world datasets, BanditPAM++ can be up to 10x faster than BanditPAM.\n", "questions": "In the beginning of intro, I got confused by the remark that k-means does not generalize to other metrics. What do you mean by that? Isn’t equation (1) such a generalization?\nAbout the PAM algorithm in general, it seems that there is no worst-case guarantee for this method. Is there any practical algorithm that also gives some kind of worst-case guarantee in terms of clustering cost?", "limitations": "yes", "rating": 6, "confidence": 2, "soundness": 2, "presentation": 3, "contribution": 2, "strengths": "I think this is an interesting and useful result. The experimental improvements seem substantial.", "weaknesses": "The techniques used to obtain this algorithm did not seem very innovative (I am not an expert in bandit algorithms however).\n\ntypos:\nLine 118: fundamentally\nLine 161: missing Delta I think\nPseudocode page 6: I guess m_r^(2) the second closest medoid ? I am not sure this is defined in this page. \n", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "We thank Reviewer mktW for their thoughtful comments and helpful feedback. We provide a point-by-point response below.\n\nResponse to W1: We would like to clarify our contributions and their novelty.\n\nBanditPAM++ improves upon BanditPAM via two novel techniques: the Virtual Arms (VA) technique and Permutation-Invariant Caching (PIC). Both techniques are tailored for BanditPAM: VA relies on the observation that the updates of $k$ medoids for a given candidate point are correlated, and PIC relies on the observation that the sequence of MAB searching problems uses the same set of reference points. Moreover, both techniques bring substantial improvements: an $O(k)$ complexity improvement in the SWAP step by VA and wall-clock time speed up by PIC. Our resulting algorithm, BanditPAM++ can be an order of magnitude faster than BanditPAM in some settings (Table 1 of the additional experiments provided in the \"global\" response to all reviewers). Since these are not general techniques, and they substantially improve the state-of-the-art $k$-medoids clustering algorithm, we argue that they are non-trivial and significant. \n\nThe VA technique is applied to accelerate the SWAP steps in BanditPAM++. Whereas the original BanditPAM algorithm considers all $kn$ possible (medoid, datapoint) pairs as \"arms\" during each SWAP step, BanditPAM++ considers only each of the $n$ datapoints as \"arms\". This is emphasized in blue in Line 1 of Algorithm 2. This methodological change is permitted by Theorem 1, which codifies our observation that at least $k-1$ of the arms from the original BanditPAM algorithm are receiving the same update in Line 7 of Algorithm 1. We show that this methodological change in BanditPAM++ results in the same answers as PAM and BanditPAM in Theorems 3 and 4, which means that BanditPAM++ matches state-of-the-art algorithms in clustering quality.\n\nCrucially, the fact that BanditPAM++ now considers every datapoint as an \"arm\" during each of the BUILD and SWAP steps permits our second observation and the PIC technique. Because the set of objects that are considered \"arms\" is the same across each phase of BanditPAM++, we can now cache the quantities that are expensive to compute: the distances between datapoints. Choosing a permutation of datapoints a priori and using this permutation across all phases of BanditPAM++ increases the proportion of cache hits when the cache is too small to hold all $n^2$ pairwise distances, as is the case in the setting we are interested in: when $n$ is large. These observations are formalized in the SPIMAB framework of our paper. The Permutation-Invariant Caching technique improves both the BUILD and SWAP procedures of BanditPAM in BanditPAM++.\n\nResponse to W2: Thank you for pointing out these typos; we will correct them in the final paper. Yes, $m^{(2)}_r$ is the second-closest medoid.\n\nResponse to Q1 (k-means): We meant to emphasize what is written in text above Equation 1, on Line 24; namely, that in Equation 1, $k$-means must use a specific form of distance metric $d(u, x_j)$ for algorithmic efficiency. Usually, this is the squared $L_2$ distance, though other very particular metrics can also be used for $k$-means.\n\nThe general point is still that $k$-means typically requires specific distance metrics for efficient solutions. In contrast, $k$-medoids is agnostic to the distance metric being used. In fact, the distance function being used in $k$-medoids need not even be a metric and can be an arbitrary dissimilarity function (i.e., may violate the triangle inequality, may be negative, and may not be symmetric). For this reason, $k$-medoids has been used to cluster \"exotic\" objects such as trees [20].\n\n\nResponse to Q1 (worst-case): One may notice that the optimization objective in Equation 2 (when negated) is actually submodular with respect to the set $\\mathcal{M}$. It is a well-known result in submodular optimization theory that greedy optimization of submodular results in a solution that is at worst $1 - \\frac{1}{e} \\approx 0.63$ of optimal (e.g., [G]).\n\nFor this reason, after the BUILD steps of PAM (or BanditPAM or BanditPAM++, since they all follow the same optimization trajectory), the loss of the current set of $k$ medoids is no worse than $\\approx 0.63$ of optimal, when the loss is appropriately negated. Since the SWAP steps of every algorithm only improve the loss, this guarantee holds for the final set of returned medoids as well.\n\nTo the best of our knowledge, this observation is novel. We have not codified it in our submission and will add this discussion to an appendix.\n\n[G: Submodular Function Maximization, Krause and Golovin]\n\nAdditional Experiments: We have also run experiments across two new datasets and metrics to demonstrate the speedups that BanditPAM++ has over BanditPAM. We provide the results for these experiments in our \"global\" rebuttal in Figure 1.\n\nOn both datasets and metrics, we analyze the scaling of BanditPAM++ with respect to both $n$ and $k$. In both additional settings, BanditPAM++ demonstrates significant speedups over BanditPAM. On the scRNA dataset, BanditPAM++ is approximately 6.4x faster than the original BanditPAM algorithm with $k=15$ and $L_1$ distance. On the 20 newsgroups dataset with $k = 15$ and cosine distance (a highly non-Euclidean space), BanditPAM++ is up to 8.5x faster than the original BanditPAM algorithm. Please see Figure 1 and Table 1 of the \"global\" rebuttal for more details.\n\nWe have also run additional experiments to understand the effects of $T$ and $\\delta$. Table 2 of our additional results demonstrates that on two different datasets and metrics, $\\delta$ can be varied significantly without affecting the quality of the clustering results. Additionally, Figure 2 of our additional experiments shows the loss vs. swap iteration ($T$) for the MNIST dataset with $L_2$ distance and the CIFAR10 dataset with $L_1$ distance, for $k=5, 10$. Increasing $T$ beyond $k$ only leads to marginal decreases in overall loss.\n"}]

Author rebuttal:
Attached please find additional experimental results that we refer to in our rebuttals to reviewers.

Meta-review: The paper presents a new algorithm for the k-medoids clustering. The algorithm improves on previous work in the area, in particular on the  state-of-the-art algorithm, the BanditPAM algorithm.

The ideas introduced in the paper are interesting and the experimental results nice but only on Euclidean datasets. This is a bit surprising because the k-medoids problem is interesting mainly in metric spaces. Nevertheless the theoretical results and the improvement on the running time are interesting and worth publishing.

For the final version we would suggest to extend the paper with:
- an in-depth discussion on the fact that the results hold on under some assumptions
- some experiment on non-euclidean datasets
Decision: Accept (poster)